require(stats)
norm2 <- function(x) Conj(t(x))%*%x;

#f<- function(t) {sin(2*pi*t);}
f <- function(t) {t*(1-t)}
#f <- function(t) {t*(1-t)*sin(2*pi*t)}

s<-0.3;
t<-seq(0, 1, 0.01); n<-length(t);
theta<-f(t);
#theta<-c(theta, theta[seq(n, 1, -1)]); t<-c(t, 1+t); n<-length(t)

y<-theta+s*rnorm(n);
plot(t, y, type="l"); Sys.sleep(0.5);


ftheta<-fft(theta)/sqrt(n);
fy<-fft(y)/sqrt(n);

kmax <- (n-1)/2; #kmax<- 10;
Cp <- array(0, kmax);
Oracle <- array(0, kmax);
Bic <- array(0, kmax);
for (k in 1:kmax){
  dm=2*k-1;
  Oracle[k] <- norm2(ftheta[(k+1):(n-k+1)]) + s^2*dm;
  
  fth <- fy;
  fth[(k+1):(n-k+1)]<- 0;
  Cp[k] <- -norm2(fth) + 2*s^2*dm;
  Bic[k] <- -norm2(fth)/(2*s^2) + dm/2*log(n);
}
#plot(1:kmax, Cp, type="l")
ko <- which.min(Oracle); print(c("k oracle = ", ko));
kCp <- which.min(Cp); print(c("k de Mallows = ", kCp));
kBic <- which.min(Bic); print(c("k Bic = ", kBic));

k<-ko;
fth <- fy;
fth[(k+1):(n-k+1)]<- 0;
th <- fft(fth, inverse=T)/sqrt(n);
lines(t, th, type="l", lwd=3, col="blue"); Sys.sleep(0.5);
lines(t, theta, type="l", lwd=3, col="red");