options(object.size=Inf)
source("react-funcs.S")
scaleab <- function(x,a=0,b=1)
{
 ##This functions takes any vector an transforms it to be between a and b
  if(a>=b) stop("a must be strictly less than b")
  o <- !is.na(x)
  aux <- rep(NA, length(x))
  if(min(x[o])==max(x[o])) aux[o] <- rep(b,length(o))
    else { 
      aux[o] <- (x[o] - min(x[o]))/diff(range(x[o]))
      aux[o] <- aux[o]*(b-a)
      aux[o] <- aux[o]+a
    }
  aux
}
scale01 <- scaleab

redact <- function(z,U,sigma2,delta=0.0000001,maxiter=50){ 
  ##returns the lambda, and fs
  ###SCALE,y and js need to be defined globaly
  yhatold <- 1
  FLAG <- T
  Count <- 1
  start.val <- 1
  while(FLAG){
    assign("sigma2",sigma2,frame=1)
    lambda <- nlminb(start.val,risk,lower=0)$param/SCALE
    f <- c(1,(1+lambda*(2*pi*js)^(2*m) )^{-1},(1+lambda*(2*pi*js)^(2*m) )^{-1},
	   (1+0.5*lambda*(pi*n)^(2*m) )^{-1})
    zhat <- f*z
    yhat <- U%*%zhat
    sigma2 <- 1/(n-sum(f^2))*sum((y-yhat)^2)
    FLAG <- sum((yhat-yhatold)^2/sum(yhatold^2)) > delta & Count < maxiter
    Count <- Count +1
    yhatold <- yhat
    start.val <- lambda*SCALE
  }
  list(lambda=lambda,fitted=yhat,sigma2=sigma2)
}

redact2 <- function(z,U,sigma2,delta=0.0000001,maxiter=50){ 
  ##returns the lambda, and fs
  ###SCALE,y and js need to be defined globaly
  yhatold <- 1
  FLAG <- T
  Count <- 1
  start.vals <- c(1,2)
  while(FLAG){
    assign("sigma2",sigma2,frame=1)
    aux <- nlminb(start.vals,risk2,lower=c(0,1))
    lambda <- aux$param[1]/SCALE
    m <- aux$param[2]
    f <- c(1,(1+lambda*(2*pi*js)^(2*m) )^{-1},(1+lambda*(2*pi*js)^(2*m) )^{-1},
	   (1+0.5*lambda*(pi*n)^(2*m) )^{-1})
    zhat <- f*z
    yhat <- U%*%zhat
    sigma2 <- 1/(n-sum(f^2))*sum((y-yhat)^2)
    FLAG <- sum((yhat-yhatold)^2/sum(yhatold^2)) > delta & Count < maxiter
    Count <- Count +1
    yhatold <- yhat
    start.vals <- c(lambda*SCALE,m)
  }
  list(lambda=lambda,m=m,fitted=yhat,sigma2=sigma2)
}

library(MASS)
xx <- mcycle$times;yy <- mcycle$accel;xx <- scale01(xx) 
aux <- loess(yy~xx,span=.25,data=mcycle,degree=1)$fitted
##objective functions:
risk <- function(lambda){
  lambda <- lambda/SCALE
  sum(  ( (1+lambda*(2*pi*js)^(2*m) )^{-1} - (1-sigma2/as^2)  )^2*as^2)+ sum( ( (1+lambda*(2*pi*js)^(2*m) )^{-1} - (1-sigma2/bs^2)  )^2*bs^2)+ ( (1+0.5*lambda*(pi*n)^(2*m) )^{-1} - (1-sigma2/an^2)  )^2*an^2
}

risk2 <- function(x){
  lambda<-x[1]/SCALE
  m<-x[2]
  sum(  ( (1+lambda*(2*pi*js)^(2*m) )^(-1) - (1-sigma2/as^2)  )^2*as^2)+ sum( ( (1+lambda*(2*pi*js)^(2*m) )^(-1) - (1-sigma2/bs^2)  )^2*bs^2)+ ( (1+0.5*lambda*(pi*n)^(2*m) )^(-1) - (1-sigma2/an^2)  )^2*an^2
}
cveval <- function(lambda){
  lambda <- lambda/SCALE
  
  sum(   (as^2+bs^2)*(1 -  (1+lambda*(2*pi*js)^(2*m) )^(-1))^2 +
      an^2*(1 - (1+0.5*lambda*(pi*n)^(2*m) )^(-1))^2) /
	(1 - sum(c(1,2*(1+lambda*(2*pi*js)^(2*m) )^(-1),
		   (1+0.5*lambda*(pi*n)^(2*m) )^(-1)))/n)^2
}
gmleval <- function(lambda){
  lambda<-lambda/SCALE
  f <- c(1,(1+lambda*(2*pi*js)^(2*m) )^{-1},(1+lambda*(2*pi*js)^(2*m) )^{-1},(1+0.5*lambda*(pi*n)^(2*m) )^{-1})
  A <- U%*%diag(f)%*%t(U)
  aux <- Re(eigen(diag(n)-A)$values)
  aux <- exp(sum(log(aux[aux>.Machine$single.eps]))/(n-m))
  matrix(y,1,n)%*%(diag(n)-A)%*%matrix(y,n,1)/aux
}

##global defintions
m <- 2
delta <- 0.000000001
SCALE <- 10^6 ##this is so that lambda isnt in the order of 10^-6
NN <- 100

###changing parameters
ns <- c(50,100,250)
sigmas <- c(.025,.05,.1,.5)


##storing arrays
mses <- array(0,dim=c(6,4,length(sigmas),length(ns)))
msevars <- array(0,dim=c(6,4,length(sigmas),length(ns)))
percentage <- array(0,dim=c(5,4,length(sigmas),length(ns)))
Lambdas <- array(0,dim=c(6,4,length(sigmas),length(ns)))
Lambdavars <- array(0,dim=c(6,4,length(sigmas),length(ns)))
Ms <- array(0,dim=c(2,4,length(sigmas),length(ns)))
Mvars <-  array(0,dim=c(2,4,length(sigmas),length(ns)))
Sigmas <- array(0,dim=c(4,4,length(sigmas),length(ns)))
Sigmavars <- array(0,dim=c(4,4,length(sigmas),length(ns)))
mse1 <- rep(0,NN);mse2<-rep(0,NN);mse3<-rep(0,NN);mse4<-rep(0,NN);mse5 <- rep(0,NN);mse6 <- rep(0,NN)
lambda1 <- rep(0,NN);lambda2 <- rep(0,NN);lambda3 <- rep(0,NN);lambda4 <- rep(0,NN);lambda5 <- rep(0,NN);lambda6 <- rep(0,NN)
sigma21 <- rep(0,NN); sigma22 <- rep(0,NN); sigma23 <- rep(0,NN);sigma24 <- rep(0,NN)
M1 <- rep(0,NN);M2 <- rep(0,NN)
for(l in 1:2){#length(ns)){
  n <- ns[l]
  ss <- matrix(0,4,n)
  x<- c(1:n)/n
  ss[1,] <- (1-abs(2*x-1)^3)^3
  ss[2,] <- sin(2*pi*x)
  ss[3,] <- approx(xx,aux,x,rule=2)$y
  ss[4,] <- 10*cos(2*pi*x) + 9.5*cos(4*pi*x+.5) + 9*cos(6*pi*x+.3) + 8.5*cos(8*pi*x+.1) + 7*sin(10*pi*x) + 4*cos(12*pi*x-.3) + 1*cos(14*pi*x-.8) + .5*cos(16*pi*x+1.1)
  ss <- t(apply(ss,1,scale01)) ##everything between 0,1
###these are global things: m is degree of sobolev space, js are the indeces 
###for the cosines and delta is for the iterative procedure stop
  js <- 1:(n/2-1)
  
  ##DEFINE THE TRANSFORM
  U <- cbind(rep(1/sqrt(n),n),
	     sapply(1:(n/2-1),function(k) sqrt(2/n)*cos(2*pi*k*x)),
	     sapply(1:(n/2-1),function(k) sqrt(2/n)*sin(2*pi*k*x)),
	     1/sqrt(n)*cos(pi*n*x))
  
  
  #for(k in 1:length(sigmas)){
  k <- 5
  sigma <- sigmas[k]
    for(i in 1:4){
      s <- ss[i,]
      for(j in 1:NN){
	y <- s + rnorm(n,0,sigma)
	##these are used for both splines and react
	z <- t(U)%*%y
	a0 <- z[1]
	as <- z[2:(n/2)]
	bs <- z[(n/2+1):(n-1)]
	an <- z[n] 
	
	sigma2.start <- 1/(n-2)/2*sum((diff(y))^2)
	aux1 <- redact(z,U,sigma2.start,delta=0.00001,maxiter=1) 
	aux2 <- redact(z,U,sigma2.start,delta=0.0000001)
	aux3 <- redact2(z,U,sigma2.start,delta=0.0000001,maxiter=1)
	aux4 <- redact2(z,U,sigma2.start,delta=0.0000001)
	lambdacv <- nlminb(1,cveval,lower=0)$param/SCALE
	lambdagml <- nlminb(1,gmleval,lower=0)$param/SCALE
	fcv <- c(1,(1+lambdacv*(2*pi*js)^(2*m) )^{-1},(1+lambdacv*(2*pi*js)^(2*m) )^{-1},(1+0.5*lambdacv*(pi*n)^(2*m) )^{-1})
	yhatcv <- U%*%(fcv*(t(U)%*%y))
	fgml <- c(1,(1+lambdagml*(2*pi*js)^(2*m) )^{-1},(1+lambdagml*(2*pi*js)^(2*m) )^{-1},(1+0.5*lambdagml*(pi*n)^(2*m) )^{-1})
	yhatgml <- U%*%(fgml*(t(U)%*%y))
	
	mse1[j] <- mean((aux1$fitted-s)^2)
	mse2[j] <- mean((aux2$fitted-s)^2)
	mse3[j] <- mean((aux3$fitted-s)^2)
	mse4[j] <- mean((aux4$fitted-s)^2)
	mse5[j] <- mean((yhatcv-s)^2)
	mse6[j] <- mean((yhatgml-s)^2)
	lambda1[j] <- aux1$l
	lambda2[j] <- aux2$l
	lambda3[j] <- aux3$l
	lambda4[j] <- aux4$l
	lambda5[j] <- lambdacv
	lambda6[j] <- lambdagml
	sigma21[j] <- sigma2.start
	sigma22[j] <- aux2$sigma
	sigma23[j] <- sigma2.start
	sigma24[j] <- aux4$sigma
	M1[j] <- aux3$m;M2[j] <- aux4$m
      }
      mses[1,i,k,l] <- mean(mse1);mses[2,i,k,l] <- mean(mse2);
      mses[3,i,k,l] <- mean(mse3);mses[4,i,k,l] <- mean(mse4);
      mses[5,i,k,l] <- mean(mse5);      mses[6,i,k,l] <- mean(mse6);
      msevars[1,i,k,l] <- var(mse1);msevars[2,i,k,l] <- var(mse2)
      msevars[3,i,k,l] <- var(mse3);msevars[4,i,k,l] <- var(mse4)
       msevars[5,i,k,l] <- var(mse5) ;      msevars[6,i,k,l] <- var(mse6)
      Lambdas[1,i,k,l] <- mean(lambda1);Lambdas[2,i,k,l] <- mean(lambda2)
      Lambdas[3,i,k,l] <- mean(lambda3);Lambdas[4,i,k,l] <- mean(lambda4)
      Lambdas[5,i,k,l] <- mean(lambda5);    Lambdas[6,i,k,l] <- mean(lambda6);
      Lambdavars[1,i,k,l] <- var(lambda1);Lambdavars[2,i,k,l] <- var(lambda2)
      Lambdavars[3,i,k,l] <- var(lambda3);Lambdavars[4,i,k,l] <- var(lambda4)
      Lambdavars[5,i,k,l] <- var(lambda5); Lambdavars[6,i,k,l] <- var(lambda6);
      Ms[1,i,k,l] <- mean(M1); Ms[2,i,k,l] <- mean(M2); 
      Mvars[1,i,k,l] <- var(M1); Mvars[2,i,k,l] <- var(M2)
      Sigmas[1,i,k,l] <- mean(sigma21);Sigmas[2,i,k,l] <- mean(sigma22);
      Sigmas[3,i,k,l] <- mean(sigma23); Sigmas[4,i,k,l] <- mean(sigma24);
      Sigmavars[1,i,k,l] <- var(sigma21);Sigmavars[2,i,k,l] <- var(sigma22);
      Sigmavars[3,i,k,l] <- var(sigma23);Sigmavars[4,i,k,l] <- var(sigma24);
      percentage[1,i,k,l] <- sum(mse1<=mse5)
      percentage[2,i,k,l] <- sum(mse2<=mse5)
      percentage[3,i,k,l] <- sum(mse3<=mse5)
      percentage[4,i,k,l] <- sum(mse4<=mse5)
      percentage[5,i,k,l] <- sum(mse6<=mse5)
      print(signif(c(i,sigma,n, mses[,i,k,l],percentage[,i,k,l]),2))
    }
}