1. Hw and final projects.
2. XCode.
3. Compiling C and calling C from R.
4. Hello world.
5. pi.
6. dnorm.
7. Sum of squared differences between observations.

1. Hw2, 3, and final projects.

-- Hw3 is on the course website http://www.stat.ucla.edu/~frederic/202a/F12 . 
Hw2 is due Tue, Oct 30, 1030am. 
Hw3 is due Tue, Nov 13, 1030am.
Late homeworks will not be accepted!

-- No class or office hours Nov 1 or Nov 6.

Final projects.

  Thur, Nov 29
1. MALEKMOHAMMADI, MAHSA .      CHUGH, KARAN .      XIA, FANGTING .      
2. BELOGOLOVA, HELEN .      WANG, YUTING .      ELLIOTT, PETER WILLIAM .      
3. LIANG, MUZHOU .      CHEN, YU-CHING .      CARLITZ, RUTH DENALI .      
4. DANIELSON, AARON .      CUNG, BIANCA .      CARLEN, JANE .      
5. HUANG, MING-CHUN .      KARYEKAR, AMIT SUHAS .      SHI, HAOLUN .      
6. CHANG, HUA-I .      TRICKEY, KIRSTEN ALEXANDRA .      SKIPPER, PETER TYREL .      

  Tues, Dec 4 
7. ZHU, HUA .      KIM, BRIAN .      LACOUR, MICHAEL JULES .     
9. HUANG, YAFEI .      CHOI, ICK KYU .      LIU, JIA LIN .      
10. GLEASON, TANIA MARIE .      LUO, YIYANG .      JASWANI, DRUMIL PARAS .      
11. LI, JINCHAO .      CHEN, RUI .      WANG, YAN .      
12. WANG, CHENHUI .      GAURAV, ADITYA .      TCHEMODANOV, NATALIA .      
13. XIE, JIANWEN .      XIAO, QIAN .      BUDHRANI, VINEESHA .      

  Thur, Dec 6 
14. CAO, RENZHE .      WANG, WENJIA .      ZHOU, YUXI .      
15. VASILYEV, ALEXANDER .      LIAN, XIAOCHEN .      TSUI, JASON .      
16. LI, LU .      NAYAK, SHILPI .  UPPALA, MEDHA .    
17. RICHMOND, ROBERT JAMISON .      LI, ZHONGBO .      ZHOU, YUAN .      
18. ARFA, JONATHAN MICHAEL .      WOLF, JEFFREY ARIEN .   LU, YANG .     
19. KAPADIA, HIRAL .      NIE, XIAOHAN .      CHEN, XIANJIE .
  
2. XCode, for Mac OSX, includes a C and C++ compiler, among other tools.
You can get XCode from 
http://wildfire.stat.ucla.edu/rick/xcode .
Click on the one choice, xcode_3.2.6_and_ios, to get it. I will only leave 
it up there for a week or so. It takes a couple hours to download.
If you have a PC, you will need a compiler. There are many free ones. 
See for instance http://www.compilers.net/dir/free/compilers/ccpp.htm .

3. Compiling C and calling C from R.

The first step to writing C code is opening a text editor.
Write your C code, call it something.c, then compile it, to create an 
object file called something.so. Then you can load that into R.
It seems to primarily work with R from terminal, but you can run C in R in 
terminal, do save.image(), and then start R, use the same working 
directory, and then do load(".RData") .

4. Hello world.

Create a C function to print "Hello world!" and call this function n 
times in R.

4a. In a text editor, create a C file.  Say it's called hello.c  
The file looks like:

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

	/* Start a comment.
	   Continue your comment.
        */

       void hello (int *n)
          {
	  int i,j;
	  double a2;
	  a2 = 4.5;
	  j = 3;
	  for(i = 0; i < *n; i++)
	     Rprintf("The %d rd integral is %f .\n", j, a2);
	  }


4b. In UNIX, in same directory where hello.c is, type

	R CMD SHLIB hello.c

or, in R, do system("R CMD SHLIB hello.c")

4c. In R, in the same directory, do

	dyn.load("hello.so")
	
	hello2 <- function(n){
		.C("hello",as.integer(n))
	}

	y = hello2(10)

5. Approximate pi in C.

In mypi.c,
    #include <R.h>
    #include <Rmath.h>
    void pi2 (int *n, double *y){
      int i;
      double x[*n];
      x[0] = 1.0;
      y[0] = sqrt(6.0);
      for(i = 1; i < *n; i++) {
          x[i] = x[i-1] + 1.0 / ((i+1.0)*(i+1.0));
          /* or x[i] = x[i-1] + 1.0 / pow(i+1.0,2.0); */
          y[i] = sqrt(6.0 * x[i]);
      }
    }
      

In R,
    system("R CMD SHLIB mypi.c")
    dyn.load("mypi.so")
    pi3 = function(n){
		.C("pi2",as.integer(n), y = double(n))
	}
    b = pi3(1000000)
    b$y[1000000]

Note that you have to be incredibly careful in C when doing arithmetic between 
integers and non-integers.
If instead of 1.0 / ((i+1.0)*(i+1.0));
you do         1  / ((i+1.0)*(i+1.0));
or            1.0 / (i+1.0)^2;
crazy stuff happens.
^ is a bitwise operator meaning "XOR", i.e. X^Y = 1 if X=1 or Y=1 but not both.
    
6. dnorm.
You can access C versions of many basic R functions, including for 
instance dnorm(), rnorm(), etc.
A nice reference is 
http://www.biostat.jhsph.edu/~bcaffo/statcomp/files/cprog3_ho.pdf .
The syntax in C of dnorm is 
double dnorm(double x, double mu, double sigma, int give_log)
 
in mydn.c,
    #include <R.h>
    #include <Rmath.h>
    void norm2 (int *n, double *upper, double *bw, double *y){
      int i;
      double x, inc;
      x = -1.0 * *upper;
      inc = 2.0 * *upper / *n;
      for(i = 0; i < *n; i++) {
      y[i] = dnorm(x / *bw, 0,1,0);
      x += inc;
      }
    }
    
In R, 
    system("R CMD SHLIB mydn.c")
    dyn.load("mydn.so")
    norm3 = function(n, u, b){
		d = .C("norm2", as.integer(n), as.double(u), 
		as.double(b), y = double(n))
	d$y
	}
    b = 100
    n = 10000
    u = 5*b
    a = norm3(n,u,b)
    title2 = paste("normal density with sd ", as.character(b))
    plot(seq(-u,u,length=n), a, type="l", main=title2,xlab="x", ylab="f(x)")

7. Sum of squared differences between observations.

In sumsq.c,

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

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

in R,
	system("R CMD SHLIB sumsq.c")
	dyn.load("sumsq.so")
	sum3 = function(data2){
		n = length(data2)
		a = .C("ss2", as.double(data2), as.integer(n),
			y=double(n))
		a$y  ## or equivalently a[[3]]
	}

	b = c(1,3,4)
        sum3(b)

	n = c(100, 1000, 2000, 3000, 5000, 7000, 8000, 10000)
	t2 = rep(0,8)
	for(i in 1:8){
	  b = runif(n[i])
	  timea = Sys.time()
	  d = sum3(b)
	  timeb = Sys.time()
	  t2[i] = timeb-timea
	  cat(n[i]," ")
	}
	plot(n,t2,ylab="time (sec)")
	
	## Now try the same thing in R, without C.
	
	sum4 = function(data2){
	  n = length(data2)
	  x = rep(0,n)
	  for(i in 1:n){
	    for(j in 1:n){
	      x[i] = x[i] + (data2[i] - data2[j])^2
	      }
	    }
	  x
	  }

	b = c(1,3,4)
        sum4(b)

	n = c(1:8)*100
	t3 = rep(0,8)
	for(i in 1:8){
	  b = runif(n[i])
	  timea = Sys.time()
	  d = sum4(b)
	  timeb = Sys.time()
	  t3[i] = timeb-timea
	cat(n[i]," ")
	}
	plot(n,t3,ylab="time (sec)")