cd Documents/NewtestsourcesRdeHal
R

library(fields)

#setting for simulation of a Gaussian Matern field on a grid n1*n1
n1<-27
#each realization of the field will be observed at the sites
#    (x_1,x_2)(i), i=1,..,n1*n1     equispaced on the unit square:
x<- (1./n1)*matrix( c(rep(1: n1, n1),rep(1: n1,each= n1)), ncol=2)
n<- nrow(x)    # n=n1*n1
# true values for the 4 parameters :
bTrue<- 1000.0    #variance (or power) of the field
# noise level =1
trueRange<-0.2
nu<-1./2    #this “smoothness index” will be assumed known
#
lags<-(1./n1)*(0:(n1-1))
autocorrelationFunct<- Matern( lags, range= trueRange, 
	smoothness=nu)
plot( lags, autocorrelationFunct, type="l")


ut<-system.time(
{
Cov.mat<-  bTrue* Matern( rdist(x,x), 
	smoothness=nu, range= trueRange)
A<- chol( Cov.mat)
})
ut

######### simulation of 9 realizations  #################
set.seed(987)
Z<- array(NA,c(n1*n1,9))
ut<-system.time(
for(indexReplcitate in 1:9){
#one simulates 9 realizations
gtrue<- t(A) %*% rnorm(n)
Z[,indexReplcitate]<-c(gtrue)
}
)
#
######### plot of these 9 realizations  #################
set.panel( 3,3)
#
x1 <- x2 <- seq(1, n1, 1)
for(indexReplcitate in 1:9){
Ztrue<-Z[,indexReplcitate]
image(x1,x2,matrix(Ztrue,n1,n1),asp=1)}

######### analysis of the 1rst of these realizations ####
set.seed(678)
indexReplcitate <-1
y <- Z[,indexReplcitate]+rnorm(n)
# variance empirique corrigée du biais
bEV<-sum(y**2)/n -1.
bEV

########  use of CGEMevalOnGrid  #########################
#generation of w used in “randomized-trace”:
set.seed(345)
w <- rnorm(n)
w<- c( w)
#
candidateRanges.Grid<- trueRange * 10**seq(-1.1,1.1,,15)
#
#ds les 2 lignes suivantes, on "source" la fonction  
# CGEMevalOnGrid,, on la charge et  on l'execute
source("CGEMevalOnGrid.R")
out <- CGEMevalOnGrid(y,w,candidateRanges.Grid,nu)
out

############ impact of w    ##############################
set.panel(1,1)
{plot( candidateRanges.Grid , out$values, type="l", 
	               col=1, lty=2,log="x")
    abline(h= bEV)
   }

set.seed(345)
for(indexReplcitateOfW in 1:6){
w <- rnorm(n)
out <- CGEMevalOnGrid(y,w,candidateRanges.Grid,nu)
candidateCGEMvalues <- out$values
lines( candidateRanges.Grid , candidateCGEMvalues, type="l",
               col=1, lty=2)
}



# use of CGEMevalOnGrid on the 9 realizations ###############
candidateRanges.Grid<- trueRange * 10**seq(-1.1,1.1,,15)
niterForY <- matrix(NA,length(candidateRanges.Grid),9)
niterForW <- matrix(NA,length(candidateRanges.Grid),9)

ut<-system.time({
	set.seed(321)
for(indexReplcitate in 1:9){
y <- Z[,indexReplcitate]+rnorm(n)
bEV<-sum(y**2)/n -1.
#
w <- rnorm(n)
w<- c( w)
# evEqualTo1.w <- c( w)*sqrt((n/sum( w *w)))
#
out <- CGEMevalOnGrid(y,w,candidateRanges.Grid,nu)
candidateCGEMvalues <- out$values
niterForY[, indexReplcitate] <-out$niterForY
niterForW[, indexReplcitate] <-out$niterForW
#
if (indexReplcitate==1) 
   {plot( candidateRanges.Grid , candidateCGEMvalues/bEV, 
   			type="l", col=1, lty=2,log="x")
   abline(h= 1, lwd=2)
   }
else lines( candidateRanges.Grid , candidateCGEMvalues/bEV, 
			type="l", col=1, lty=2)
#
}
})
niterForY
niterForW
##################################################################



# another use of CGEMevalOnGrid: thes estimating equation #######
#     "  candidateCGEMvalues  / bEV  = 1      " #################
# can be solved wrt the microergodic parameter  #################
#
#a grid centered around the true microergodic parameter
# note the factor 2 in place of 10
microergodic.Coef.Grid<- 2**seq(-1,1,,15)
Coef.true<- bTrue*(1/trueRange)**(2*nu)
microergodic.Coef.Grid<-Coef.true  *microergodic.Coef.Grid
#
set.seed(321)
set.panel(1,1)
for(indexReplcitate in 1:9){
y <- Z[,indexReplcitate]+rnorm(n)
bEV<-sum(y**2)/n -1.
#
w <- rnorm(n)
w<- c( w)
# definition of a new grid
candidateRanges.Grid<- (bEV/microergodic.Coef.Grid)**(1/(2* nu))
out <- CGEMevalOnGrid(y,w,candidateRanges.Grid,nu)
niterForY[, indexReplcitate] <-out$niterForY
niterForW[, indexReplcitate] <-out$niterForW
#
candidateCGEMvalues <- out$values
if (indexReplcitate==1) 
   {plot( microergodic.Coef.Grid/Coef.true , 
   			candidateCGEMvalues/bEV, type="l", col=1, lty=2,log="x")
   abline(h= 1, lwd=2)
   }
else lines( microergodic.Coef.Grid/Coef.true , 
		candidateCGEMvalues/bEV, type="l", col=1, lty=2)
#
}
################################################################




# one repeats the above computation but with Normalized w's  ####
# another use of CGEMevalOnGrid: thes estimating equation #######
#     "  candidateCGEMvalues  / bEV  = 1      " #################
# can be solved wrt the microergodic parameter  #################
#
#a grid centered around the true microergodic parameter
# note the factor 2 in place of 10
microergodic.Coef.Grid<- 2**seq(-1,1,,15)
Coef.true<- bTrue*(1/trueRange)**(2*nu)
microergodic.Coef.Grid<-Coef.true  *microergodic.Coef.Grid
#
set.seed(321)
set.panel(1,1)
for(indexReplcitate in 1:9){
y <- Z[,indexReplcitate]+rnorm(n)
bEV<-sum(y**2)/n -1.
#
w <- rnorm(n)
w<- c( w)
#
# normalization of the vector w
w <- w * sqrt((n/sum( w *w)))
#
# definition of a new grid
candidateRanges.Grid<- (bEV/microergodic.Coef.Grid)**(1/(2* nu))
out <- CGEMevalOnGrid(y,w,candidateRanges.Grid,nu)
niterForY[, indexReplcitate] <-out$niterForY
niterForW[, indexReplcitate] <-out$niterForW
#
candidateCGEMvalues <- out$values
if (indexReplcitate==1) 
   {plot( microergodic.Coef.Grid/Coef.true , 
   			candidateCGEMvalues/bEV, type="l", col=1, lty=2,log="x")
   abline(h= 1, lwd=2)
   }
else lines( microergodic.Coef.Grid/Coef.true , 
		candidateCGEMvalues/bEV, type="l", col=1, lty=2)
#
}
################################################################





############ impact of w (normalized) for the 1rst realiz ###########
set.panel(1,1)
set.seed(321)
indexReplcitate<-1
y <- Z[,indexReplcitate]+rnorm(n)
bEV<-sum(y**2)/n -1.
#
candidateRanges.Grid<- (bEV/microergodic.Coef.Grid)**(1/(2* nu))
w <- rnorm(n)
w <- c( w)*sqrt((n/sum( w *w)))
out <- CGEMevalOnGrid(y,w,candidateRanges.Grid,nu)
{plot( microergodic.Coef.Grid/Coef.true  , out$values, type="l", 
	               col=1, lty=2,log="x")
    abline(h= bEV)
   }

set.seed(345)
for(indexReplcitateOfW in 1:9){
w <- rnorm(n)
w <- c( w)*sqrt((n/sum( w *w)))
out <- CGEMevalOnGrid(y,w,candidateRanges.Grid,nu)
candidateCGEMvalues <- out$values
lines( microergodic.Coef.Grid/Coef.true  , candidateCGEMvalues, type="l",
               col=1, lty=2)
}
############      end of "impact of w (normalized)"        ###########



############ impact of w (un-normalized) for the 1rst realiz ###########
set.panel(1,1)
set.seed(321)
indexReplcitate<-1
y <- Z[,indexReplcitate]+rnorm(n)
bEV<-sum(y**2)/n -1.
#
candidateRanges.Grid<- (bEV/microergodic.Coef.Grid)**(1/(2* nu))
w <- rnorm(n)
out <- CGEMevalOnGrid(y,w,candidateRanges.Grid,nu)
{plot( microergodic.Coef.Grid/Coef.true  , out$values, type="l", 
	               col=1, lty=2,log="x")
    abline(h= bEV)
   }

set.seed(345)
for(indexReplcitateOfW in 1:8){
w <- rnorm(n)
out <- CGEMevalOnGrid(y,w,candidateRanges.Grid,nu)
candidateCGEMvalues <- out$values
lines( microergodic.Coef.Grid/Coef.true  , candidateCGEMvalues, type="l",
               col=1, lty=2)
}
############      end of "impact of w (un-normalized)"      ###########






################################################################
#########  use of CGEMbisectionLogScaleSearch            #######
################################################################
set.seed(678)
indexReplcitate <-1
y <- Z[,indexReplcitate]+rnorm(n)

bEV<-sum(y**2)/n -1.
bEV

w <- rnorm(n)
w<- c( w)
w <- w * sqrt(n/sum( w *w))

#
#ds les 2 lignes suivantes, on "source" la fonction CGEMevalOnGrid, 
# on la charge et  on l'execute
source("CGEMbisectionLogScaleSearch.R")
out <- CGEMbisectionLogScaleSearch(y,w,nu,0.01,30.,0.0001)
out


rangeHatCGEMEV <- out$root
cHatCGEMEV<-  bEV*(1/rangeHatCGEMEV)**(2*nu)
cHatCGEMEV


################################################################
#########  idem as above but with       bTrue=10         #######
################################################################
bTrueOld<-bTrue
bTrue<- 10.0 
Z<- Z*sqrt(bTrue/bTrueOld)

####################################
set.seed(678)
indexReplcitate <-1
y <- Z[,indexReplcitate]+rnorm(n)
bEV<-sum(y**2)/n -1.
bEV
 
w <- rnorm(n)
w<- c( w)
w <- c( w)*sqrt((n/sum( w *w)))
#
out <- CGEMbisectionLogScaleSearch(y,w,nu,0.01,30.,0.0001)
out

################################################################
################################################################
################################################################
#########  idem as above but with       trueRange=1.         ###
################################################################
trueRange=1. 
autocorrelationFunct<- Matern( lags, range= trueRange, 
	smoothness=nu)
plot( lags, autocorrelationFunct, type="l")

ut<-system.time(Cov.mat<-  bTrue* Matern( rdist(x,x), 
		smoothness=nu, range= trueRange))
ut<-system.time(A<- chol( Cov.mat))

######### simulation of 9 realizations  #################
set.seed(987)
Z<- array(NA,c(n1*n1,9))
ut<-system.time(
for(indexReplcitate in 1:9){
#one simulates 9 realizations
gtrue<- t(A) %*% rnorm(n)
Z[,indexReplcitate]<-c(gtrue)
}
)
#
######### plot of these 9 realizations  #################
set.panel( 3,3)
#
x1 <- x2 <- seq(1, n1, 1)
for(indexReplcitate in 1:9){
Ztrue<-Z[,indexReplcitate]
image(x1,x2,matrix(Ztrue,n1,n1),asp=1)}


############################################################
######### analysis of the 1rt realization  #################
############################################################
set.seed(678)
indexReplcitate <-1
y <- Z[,indexReplcitate]+rnorm(n)
bEV<-sum(y**2)/n -1.
bEV
# 
w <- rnorm(n)
w<- c( w)
w <- c( w)*sqrt((n/sum( w *w)))


out <- CGEMbisectionLogScaleSearch(y,w,nu,0.01,30.,0.0001)
out


############################ case bTrue = 1000.   #########
bTrueOld<-bTrue
bTrue<- 1000.0 
Z<- Z*sqrt(bTrue/bTrueOld)

####################################
set.seed(678)
indexReplcitate <-1
y <- Z[,indexReplcitate]+rnorm(n)
bEV<-sum(y**2)/n -1.
bEV
# 
w <- rnorm(n)
w<- c( w)
w <- c( w)*sqrt((n/sum( w *w)))
out <- CGEMbisectionLogScaleSearch(y,w,nu,0.01,30.,0.0001)
out

rangeHatCGEMEV <- out$root
cHatCGEMEV<-  bEV*(1/rangeHatCGEMEV)**(2*nu)
cHatCGEMEV





#################### ''extensive'' simulation ###########
#
## use of CGEMbisectionLogScaleSearch ###################
########### for 500 replicates              #############
#########################################################
#### cholesky decomp of the covariance matrix  ##########
ut<-system.time(Cov.mat<-  bTrue* Matern( rdist(x,x), 
	smoothness=nu, range= trueRange))
ut
ut<-system.time(A<- chol( Cov.mat))
ut
#########################################################
####### generation of 500 replicates            #########
#######  and  CGEM-EV estimates for each one   ##########
bHatEV<-matrix(NA,500)
cHatCGEMEVLogScaleSearch<-matrix(NA,500)
nCGiterationsMaxForYLogScaleSearch<-matrix(NA,500)
nCGiterationsMaxForWLogScaleSearch<-matrix(NA,500)
ut<-system.time(
for(indexReplcitate in 1:500){
set.seed(indexReplcitate)  # so that it is reproducible #
y <- c(t(A) %*% rnorm(n))  +rnorm(n)  
bEV<-sum(y**2)/n -1.
#
w <- rnorm(n)
w<- c( w)
w <- c( w)*sqrt((n/sum( w *w)))
#
out <- CGEMbisectionLogScaleSearch(y,w,
				nu,0.01,30.,0.0001)
cHatCGEMEVLogScaleSearch[indexReplcitate]<-  
				bEV*(1/out$root)**(2*nu)
nCGiterationsMaxForYLogScaleSearch[indexReplcitate]<- 
	 max(out$niterCGiterationsHistory[,1])
nCGiterationsMaxForWLogScaleSearch[indexReplcitate]<-  
	max(out$niterCGiterationsHistory[,2])
#
}
)
ut
ctrue<-bTrue*(1/trueRange)**(2*nu)
summary(cHatCGEMEVLogScaleSearch/ctrue)
sd(cHatCGEMEVLogScaleSearch/ctrue)
nCGiterationsMaxForYLogScaleSearch