options(object.size=10^100)
###You need the function my.nls which can be obtained from same web-page
###First we define functions.. Window is gaussian for ease of computation
###scale01 scales thing to be between 0 and 1
##We made a mistake in Window..we multiply by sigma instead of dividing
##but the rest of program was written with this in mind
source("my.nls.S")
Window <- function(o,middle=pi/2,sigma=3){
  dnorm((o-middle)/max(abs(o-middle))*sigma)
}
scale01 <- 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
}
TRUENH <- 2
A <- 2*c(40,7) ## with A and B we deinf rho and psi.
B <- 2*c(8,25) ## use trigonomatric function cos(A+B)=COS(A)COS(B)+etc...
dur <- .04 
n1 <- 660
n2 <- 660*2^(1/12)*flex
change <- .0250
RATE <- 44100/4 #We divide by four so that program runs 16 times faster!!!
tt <- seq(0,dur,len=RATE*dur)
N <- length(tt)
l <- c(rep(n1,round(change*RATE)),seq(n1,n2,len=N-round(change*RATE)))
COS <- sapply(1:TRUENH,function(k) A[k]*cos(k*2*pi*l*tt))
SIN <- sapply(1:TRUENH,function(k) B[k]*sin(k*2*pi*l*tt))
truemean <- apply(COS,1,sum)+apply(SIN,1,sum)
hsn <- 20 ### to be fair we consider 10 values of h in each of the 3 sections
hs <- c(c(1,3,5,7,8,9,10),seq(11,20,len=hsn),seq(21,50,len=hsn))
NN <- 5000 ###number of simulations
res <- array(0,dim=c(length(hs),3*7,NN)) ##here will save the prelim results
for(i in 1:length(hs)){
  cat("\n",round(i/length(hs)*100))
  h <- hs[i]
  w <- Window(tt,middle=dur/2,h)
  w <- w/w[(N+1)/2]
  W2 <- sum(w^2)/N  #THESE ARE APPROX TO INTEGRALS
  W0 <- sum(w)/N
  U0 <- sum(seq(0,1,len=N)*w)/N
  w <- Window(tt,middle=dur/2,h)
  w[w<dnorm(3)] <- 0
  w <- w/sum(Window(tt,middle=dur/2,min(hs)))*N
  SW <- sum(w)
  res[i,,] <-  sapply(1:NN,function(i){ 
    cat(".")
    y <- truemean + rnorm(N,0,1)
    aux <- data.frame(yy=sqrt(w)*y,ww=sqrt(w),tt=tt)
    parameters(aux) <- list(l0=n1)
    NH <- 1
    fit1 <-my.nls(yy~cbind(ww*cos(2*pi*l0*tt),ww*sin(2*pi*l0*tt)),
	       data=aux,algorithm="plinear",start=list(l0=parameters(aux)$l0))
    resid <- fit1$resid[aux$ww>0]/aux$ww[aux$ww>0]
    NP <- U0/W0*2*(NH+1) 
    AMPS <-fit1$param[seq(2,2*NH,2)]^2 + fit1$param[seq(3,2*NH+1,2)]^2
    BICp <- 2*NH*log(SW/2)+log(N*W2*sum(AMPS*c(1:NH)^2)/2) 
    sigma2hat <- log(sum(resid^2)/(N-2))
    wsigma2hat <- log(sum(w[aux$ww>0]*resid^2)/(SW - NP) )
    MSE1 <- (truemean[(N+1)/2] - fit1$fitted[(N+1)/2]/aux$ww[(N+1)/2])^2 
    AIC1 <- sigma2hat + (2*NH+1)/N
    WAIC1 <- wsigma2hat + NP/SW
    BIC1 <- sigma2hat + log(N)/N
    WBIC1 <- wsigma2hat + BICp/SW
    CAICF1 <- sigma2hat + log(N)/N + (2*NH+1)/N
    WCAICF1 <- wsigma2hat + BICp/SW + NP/SW
    NH<-2########K=2
    fit1<-my.nls(yy~cbind(ww*cos(2*pi*l0*tt),ww*sin(2*pi*l0*tt),
		       ww*cos(4*pi*l0*tt),ww*sin(4*pi*l0*tt)),
	      data=aux,algorithm="plinear",start=list(l0=parameters(aux)$l0))
    resid <- fit1$resid[aux$ww>0]/aux$ww[aux$ww>0]
    NP <- U0/W0*2*(NH+1)
    AMPS <-fit1$param[seq(2,2*NH,2)]^2 + fit1$param[seq(3,2*NH+1,2)]^2
    BICp <- 2*NH*log(SW/2)+log(N*W2*sum(AMPS*c(1:NH)^2)/2) 
    sigma2hat <- log(sum(resid^2)/(N-2))
    wsigma2hat <- log(sum(w[aux$ww>0]*resid^2)/(SW - NP) )
    MSE2 <- (truemean[(N+1)/2] - fit1$fitted[(N+1)/2]/aux$ww[(N+1)/2])^2
    AIC2 <- sigma2hat + (2*NH+1)/N
    WAIC2 <- wsigma2hat + NP/SW
    BIC2 <- sigma2hat + log(N)/N
    WBIC2 <- wsigma2hat + BICp/SW
    CAICF2 <- sigma2hat + log(N)/SW + (2*NH+1)/SW
    WCAICF2 <- wsigma2hat + BICp/SW + NP/SW
    NH <-3 #### K=3
    fit1<-my.nls(yy~cbind(ww*cos(2*pi*l0*tt),ww*sin(2*pi*l0*tt),
		       ww*cos(4*pi*l0*tt),ww*sin(4*pi*l0*tt),
		       ww*cos(6*pi*l0*tt),ww*sin(6*pi*l0*tt)),
	      data=aux,algorithm="plinear",start=list(l0=parameters(aux)$l0))
    resid <- fit1$resid[aux$ww>0]/aux$ww[aux$ww>0]
    NP <- U0/W0*2*(NH+1) 
    AMPS <-fit1$param[seq(2,2*NH,2)]^2 + fit1$param[seq(3,2*NH+1,2)]^2
    BICp <- 2*NH*log(SW/2)+log(N*W2*sum(AMPS*c(1:NH)^2)/2) 
    sigma2hat <- log(sum(resid^2)/(N-2))
    wsigma2hat <- log(sum(w[aux$ww>0]*resid^2)/(SW - NP) )
    MSE3 <- (truemean[(N+1)/2] - fit1$fitted[(N+1)/2]/aux$ww[(N+1)/2])^2
    AIC3 <- sigma2hat + (2*NH+1)/N
    WAIC3 <- wsigma2hat + NP/SW
    BIC3 <- sigma2hat + log(N)/N
    WBIC3 <- wsigma2hat + BICp/SW
    CAICF3 <- sigma2hat + log(N)/SW + (2*NH+1)/N
    WCAICF3 <- wsigma2hat + BICp/SW + NP/SW
    c(MSE1,AIC1,WAIC1,BIC1,WBIC1,CAICF1,WCAICF1,
      MSE2,AIC2,WAIC2,BIC2,WBIC2,CAICF2,WCAICF2,
      MSE3,AIC3,WAIC3,BIC3,WBIC3,CAICF3,WCAICF3)
  })
}
###NOW WE PREPARE THAT TABLE... its be in tab 
aux <- array(0,dim=c(6,2,500))
tab <- array(0,dim=c(3,3,6))
h1 <- 10 ##THESE ARE USED TO DEFINE THE 3 segments of h values
h2 <- 20

for(k in 1:6){
    aux <- apply(res[,(0:2)*7+k+1,],3,function(x){#as.vector(x)}) 
###what is (0:2)*3+k+2? k represents what criterion.2 is AIC, 3 is WAIC, etc..
###we get a value for each simulations and for each model. 0:2 are the the 3
###models, we have to skip 7.they stored like this MSE1,AIC1,...,MSE2,AIC2,....
      x<-order(as.vector(x))[1]
      i <- ceiling(x/length(hs))
      c(i,hs[x - (i-1)*length(hs)])})
    tab[,,k] <- t(sapply(1:3,function(i){
      c(sum(aux[1,]==i & aux[2,]<h1)/NN*100,
	sum(aux[1,]==i & aux[2,] < h2 & aux[2,]>=h1)/NN*100,
	�     ##sum(aux[1,]==i & aux[2,] < h3 & aux[,2]>=h2)/NN*100)})
	sum(aux[1,]==i & aux[2,] >= h2)/NN*100)}))
}

print(tab)
### THE FINAL RESULTS ARE in tab...