forcing-perturbation/eigentechniques/unit_testing/test_EOF.R

source('../CFU_Rfunc.txt')
library(s2dverification)
#
# A little testing script by Virginie
# 

nt=6
ny=5
nx=4
neofs=min(nt,nx*ny)
a=array(rnorm(nt*ny*nx),dim=c(nt,ny,nx))
a=a-InsertDim(Mean1Dim(a,1),1,nt)
lon=c(20,25,30,35)
lat=c(20,25,30,35,40)
a[c(18,22)]=NA

# Weight computation
#wgt=array(cos(lat*pi/180),dim=c(ny,nx))
#wgt <- wgt/sum(wgt,na.rm=T)
#wgt <- sqrt(wgt)
#wgt <- array(1,dim=c(ny,nx))

#Within CFU_EOF
b=CFU_EOF(lon,lat,a,neofs=neofs)

#Without CFU_EOF
#abis=array(a*InsertDim(wgt,1,nt),dim=c(nt,ny*nx))
#b=La.svd(abis)
#EOF=array(b$v,dim=c(neofs,ny,nx))
#PC=b$u

test=array(0,dim=dim(a))
for (i in 1:nt) {
  for (j in 1:neofs) {
    # Check within CFU_EOF
    #test[i,,]=test[i,,]+b$PCs[i,j]*b$EOFs[j,,]*b$D[j]/wgt
    #test[i,,]=test[i,,]+b$PCs[i,j]*b$EOFs[j,,]/wgt # new PC normalization
    test[i,,]=test[i,,]+b$PCs[i,j]*b$EOFs[j,,] # new EOF and PC normalization

    # Check without CFU_EOF
    #test[i,,]=test[i,,]+PC[i,j]*EOF[j,,]*b$d[j]/wgt 
  }
}

print("Restitution of the original field after decomposing by CFU_EOF")
print("The difference below should be of the order of the machine precision:")
print(max(abs(test-a),na.rm=T))

test2=array(0,dim=c(nt,neofs))
for (j in 1:neofs) {
  proj=CFU_MODE(InsertDim(InsertDim(InsertDim(a,1,1),2,1),4,1),b,mode=j)
  for (i in 1:nt) {
    # Check within CFU_MODE
    test2[i,j]=test2[i,j]+proj$PC.ver[1,1,i,1]

    # Check without CFU_MODE
    #test2[i,j]=test2[i,j]+sum(a[i,,]*EOF[j,,]/b$d[j]/wgt)
  }
}
print("Restitution of the PC from CFU_EOF when projecting with CFU_MODE")
print("The difference below should be of the order of the machine precision:")
#diag=test2-PC
diag=test2-b$PCs
#diag[which(1e-14>abs(InsertDim(b$d,1,nt)))]=0 # If the eigenvalue is of the order 
# of the machine precision ~ 1e-15 then dividing by it gives a huge PC by projection
# hence I do not consider the reconstruction in such case. 
diag[which(1e-5>abs(InsertDim(b$Var,1,nt)))]=0 # If the eigenvalue is of the order 
# of the machine precision ~ 1e-15 then dividing by it gives a huge PC by projection
# hence I do not consider the reconstruction in such case. 
print(max(abs(diag)))