1. Sum of squared differences from observations for each gridpoint.
2. Kernel regression. 
3. Vectors and matrices in C. 
4. Calling C functions from C. 
5. Running C from terminal. 
6. Reading in from a file. 

1. Sum of squared differences from observations for each gridpoint.

In sumsgrid.c,

	#include <R.h>
	#include <Rmath.h>

	void ssg2 (double *x, int *n, double *g, int *m, double *y)
	/* x will be the vector of data of length n,
	   g will be the vector of m gridpoints,
	   and y will be a vector of length m, of squared differences 
	   from gridpoint i to the n observations.
	*/
	{
		int i;
		int j;
		double a;
		for(i = 0; i < *m; i++){
			a = 0.0;
			for(j=0; j < *n; j++){
				a += pow(g[i] - x[j], 2);
			}
			y[i] = a;
		}
	}

in R,
	system("R CMD SHLIB sumsgrid.c")
	dyn.load("sumsgrid.so")
	sum4 = function(data2, grid2){
		n = length(data2)
		m = length(grid2)
		a = .C("ssg2", as.double(data2), as.integer(n),
			as.double(grid2), as.integer(m),
			y=double(m))
		a$y  ## or equivalently a[[5]]
	}

	b = c(1,3,4)
	gr = c(0,1,2,3,4,5)
        sum4(b,gr)

This function ssg2 will be similar to your kernel density code, except 
that you will also have a bandwidth as an additional input, and you will 
replace the line a += pow(g[i] - x[j], 2); with the kernel density 
estimate.

2. Kernel regression.

Given vectors X and Y of length n, suppose you want to look at the 
relationship between X and Y, without assuming a linear model or other 
model. You could just plot(x,y), but perhaps you want to look at a smooth 
summary of the relationship between X and Y.

Suppose Y_i = f(X_i) + eps_i, for i = 1 to n,
where f is some smooth function, and where eps_i are iid with mean 0.

The Nadaraya-Watson or kernel regression estimate of f is 
fhat(x) = Sum [Y_i k(x-X_i; b)] / Sum [k(x-X_i; b)],
where k is some kernel function,
b is some bandwidth,
and the sums are from i = 1 to n.

Note that k(x-X_i; b) = k{(x-X_i)/b)}/b.
[for instance, compare dnorm(7.0/33.3)/33.3 with dnorm(7.0, sd=33.3).]

A commonly used kernel is the Epanechnikov kernel, 
k(u) = 3/4 (1-u^2), if |u| ≤ 1, and 0 otherwise.

This was proposed in 
  Epanechnikov, V.A. (1969). "Non-parametric estimation of a multivariate 
  probability density". Theory of Probability and its Applications 14: 153–158.
and shown to have minimal variance under certain conditions.

in epanech.c, 

#include <R.h>
#include <Rmath.h>

/* x and y will be the data vectors each of length n, 
b is the bandwidth,
g will be the vector of m gridpoints where we want to calculate the 
kernel regression estimates, est, which is thus also of length m. 
*/

void kernreg2 (double *x, double *y, int *n, double *b,
	       double *g, int *m, double *est)
{
    int i,j;
    double a1,a2,c;
    for(i = 0; i < *m; i++){
	a1 = 0.0;
	a2 = 0.0;
	for(j=0; j < *n; j++){
	    if(fabs((x[j] - g[i])/ *b) <= 1.0){
		c = 0.75 * (1.0-pow((x[j]-g[i])/ *b,2.0))/ *b;
		a1 += y[j] * c;	
		a2 += c;
	    }
	}
	if(a2 > 0.0) est[i] = a1/a2; 
	else est[i] = 0.0;
	Rprintf("I'm on gridpoint %f and index %d right now.\n", g[i], i);
    }
}

// Watch out for /* in your code. This is treated as a comment started.
// Use / *b instead of /*b.
// These double slashes can be used in C to comment an individual line, 
// like # in R.

In R,

	system("R CMD SHLIB epanech.c")
	dyn.load("epanech.so")
	n = 1000
	x = runif(n)*4
	y = 0.1*x^5 - 2*x^3 + 10*x + rnorm(n)
	bw = bw.nrd(x)
	m = 100
	g2 = seq(0,4,length=m)
## This time, I'm going to call the C function directly, 
## without making a separate R function to call it.
	a = .C("kernreg2", as.double(x), as.double(y), as.integer(n), 
	as.double(bw), as.double(g2), as.integer(m), est=double(m))
	plot(c(0,4),c(-3,17),type="n",xlab="x",ylab="y")
	points(x,y,pch=".",col="red")
	lines(g2,a$est,col="black")
	lines(g2,.1*g2^5 - 2*g2^3 + 10*g2,col="blue")

Here, bw ~ 0.3. Try it again with bw = .01 or bw = 2.
The Nadaraya-Watson estimator tends to be biased at the edges, 
especially when the bandwidth is large.

Play around with a[[1]], a[[2]], etc.
C functions called by R output back to R a list
containing all the arguments passed to them, though here the last argument, 
a$est, is changed of course after running the function kernreg2.

Note that, when calling your C function, the order of the arguments 
matters, and their type matters, but their names don't matter. For 
instance, here in R I called the bandwidth bw, but in C I called it b.
Same with g2 and g. 

3. Vectors and matrices in C.

You can just say double x[100] or double x[25][4], for example.
In mypi2.c,
  
  /* Calculate the first 100 terms in the approximation of pi. */
  #include <R.h>
  #include <Rmath.h>
  void pi100 (double *y){
    int i;
    double x[100];

    x[0] = 1.0;
    y[0] = sqrt(6.0);
    for(i = 1; i < 100; i++) {
      x[i] = x[i-1] + pow(i+1.0,-2.0);
      y[i] = sqrt(6.0 * x[i]);
    }
  }

In R,
    system("R CMD SHLIB mypi2.c")
    dyn.load("mypi2.so")
    a = .C("pi100",y=double(100))
    plot(1:100,a$y,type="l",xlab="n",ylab="Approx. of pi")

In mypi3.c,
  
  /* Calculate 1*2*3 + 4*5*6 + ... + 28*29*30. */
  #include <R.h>
  #include <Rmath.h>
  void weirdsum (double *y){
    int i,j;
    double x[10][3];
    *y = 0.0;
    for(i = 0; i < 10; i++) {
      for(j = 0; j < 3; j++){
        x[i][j] = 3*i + j + 1.0;
      }
      *y += x[i][0] * x[i][1] * x[i][2];
    }
  }

In R,
    system("R CMD SHLIB mypi3.c")
    dyn.load("mypi3.so")
    y = 1.0
    a = .C("weirdsum",as.double(y))
    a

If you want x to be an n by m matrix, where n and m are inputs, you can 
just do, for instance, in mypi4.c,

  /* Calculate 1*2*3*...*m + (m+1)(m+2)...(2m) + ... + [(n-1)m+1][(n-1)m+2]...(nm). */
  #include <R.h>
  #include <Rmath.h>
  void pi103 (double *y, int *n, int *m){
    int i,j;
    double x[*n][*m];
    double a;
    
    /* First define the matrix. */
    for(i = 0; i < *n; i++) {
	for(j = 0; j < *m; j++){
	    x[i][j] = *m*i + j + 1.0;
	}
    }
    
    /* Now multiply each row and add the product to y. */
    *y = 0.0;
    for(i = 0; i < *n; i++){
	a = 1.0;
	for(j = 0; j < *m; j++){
	    a *= x[i][j];
	}
	*y += a;
      }
    } 

In R,
    system("R CMD SHLIB mypi4.c")
    dyn.load("mypi4.so")
    y = 1.0
    a = .C("pi103",as.double(y),as.integer(10),as.integer(3))
    a[[1]]
    ## to check, 
    prod(1:3)+prod(4:6)+prod(7:9)+prod(10:12)+
    prod(13:15)+prod(16:18)+prod(19:21)+prod(22:24)+prod(25:27)+prod(28:30)
    a = .C("pi103",result75=as.double(y),as.integer(8),as.integer(4))
    a$result75
    prod(1:4)+prod(5:8)+prod(9:12)+prod(13:16)+
    prod(17:20)+prod(21:24)+prod(25:28)+prod(29:32)

Note that you have to say "y = as.double(y)" rather than just as.double(y) 
if you want to be able to call the result with a$y.

To do this in R, you could do
  n = 8
  m = 4
  x = matrix(1:(n*m),byrow=T,ncol=m)
  y = apply(x,1,prod)
  sum(y)

What if you want to define a vector of length n, or a matrix with 
a dimension n, and n is not an input or a known integer but instead 
something that must be calculated inside your function, or something that 
changes as you go through your function? Then you need to do dynamic 
memory allocation. We will hopefully get to this later.

4. Calling C functions from C.

In mysd1.c,
  
  /* Create a 10 x 4 matrix of integers and compute the sample SD of each row. */
  /* Row 1 will be (1, 4, 9, 16). Row 2 will be (25, 36, 49, 64), etc. */

  #include <R.h>
  #include <Rmath.h>
  
  double sd2(double *x, int n){
    int i;
    double xbar, b;
    xbar = 0.0;
    for(i = 0; i < n; i++){
      xbar += x[i];
    }
    xbar = xbar/n;
    b = 0.0; // b will be the sum of squared deviations
    for(i = 0; i < n; i++){
      b += pow(x[i] - xbar, 2);
    }
    return(pow(b/(n-1),0.5));
    }

  void rowsd (double *y){
    int i,j;
    double a;
    double x[10][4];
    for(i = 0; i < 10; i++) {
      for(j = 0; j < 4; j++){
        x[i][j] = pow(4*i + j + 1.0, 2);
      }
      y[i] = sd2(x[i], 4);
    }
  }
 
In R,
    system("R CMD SHLIB mysd1.c")
    dyn.load("mysd1.so")
    a = .C("rowsd",y=double(10))
    a$y
    ## to check it,
    sd(c(1,4,9,16))
    sd(c(25,36,49,64))

5. Running C from terminal.

In mypi5.c,
  
  #include <stdio.h>
  #include <math.h>

  int main ( ){
    int i;
    float x,y,h;
    printf("\n How many years old are you? ");
    scanf("%d", &i);
    printf("\n How long is the first side of your right triangle? ");
    scanf("%f", &x);
    printf("\n How long is the second side of your right triangle? ");
    scanf("%f", &y);
    h = sqrtf(x*x + y*y);
    printf("\n You are %d years old and the hypotenuse is %f . \n", i, h);
    return(i);
  }

In terminal on mac (or powershell on PC),
  cd /Users/fredericschoenberg/Documents/2012/202aF12/CinR,
  gcc mypi5.c -o mypi5.out  
  ./mypi5.out

6. Reading in from a file.

You first need a file pointer
  FILE *fp;

FILE is a structure that holds information about the file. 
Use FILE * to pass the information around by reference.

  fp = fopen(filename, mode);

filename is a string that holds the name of the file on disk 
(including a path like /cs/course if necessary).
mode is a string representing how you want to open the file. 
Most often you'll open a file for reading ("r") or writing ("w").

In data.txt,
  1
  14
  5
  37
  7
  2
  1
  100

In readin1.c,
  
  #include <stdio.h>
  #include <R.h>
  #include <Rmath.h>

  void readssq (int *n){
    /* Read in the first n terms in "data.txt"
    and compute the sum of squares.
    Write the output to "output.txt" */
    FILE *myfile1, *myfile2;
    int i;
    float a1, a2;
    
    a2 = 0.0;
    myfile1 = fopen("data.txt", "r");
    myfile2 = fopen("output.txt", "w");
    for(i = 0; i < *n; i++){
	fscanf(myfile1, "%f", &a1);
	a2 += pow(a1, 2);
    }
    fclose(myfile1);
    Rprintf("The sum of squares is %f \n", a2);
    fprintf(myfile2, "%f \n", a2);
    fclose(myfile2);
  }

In R,
    system("R CMD SHLIB readin1.c")
    dyn.load("readin1.so")
    a = .C("readssq",as.integer(5))

The input file that we are opening for reading ("r") 
must already exist, 
but the output file we are opening for writing ("w") 
does not have to exist. 
If it doesn't, it will be created. If this output file 
does already exist, 
its previous contents will be thrown away and lost.
Use "a" for append to write something to the end of a 
file without losing 
its contents.

In fprintf,
Integer values use %d,
Floating point values, i.e. decimals, use %f,
%lf for double precision numbers,
Single characters use %c,
Character strings use %s.

Coming up.
## inputting and outputting matrices.
## Generalized additive models.
## Making R packages.
## Calling R from C.
## C++ and calling C++ from R.