library(R6) rkhs <- R6Class("rkhs", public = list( y = NULL, t = NULL, b = NULL, lambda=NULL, ker= NULL, initialize = function(y = NULL, t=NULL,b=NULL,lambda=NULL,ker=NULL) { self$y = y self$t = t self$b = b self$lambda = lambda self$ker = ker self$greet() }, greet = function() { cat(paste("RKHS ",self$ker$greet(),".\n")) }, showker = function () { cat(paste0("ker is", self$ker$greet(), ".\n")) }, predict = function() { ########### calculate optimised lengthscale ##################################################### y = scale(as.matrix(self$y), center=TRUE, scale=FALSE) mean_y= apply(as.matrix(self$y),2,mean) t = as.numeric(self$t) n = length(self$t) # ############################## make prediction for derivative ###################################################### y_t = array(c(0),n)#matrix(c(0),ncol = n,nrow=dimstate) z_t = array(c(0),n)#matrix(c(0),ncol = n,nrow=dimstate) for (i in 1:n) { y_t[i] = self$ker$kern(t(t),t[i])%*%self$b +mean_y z_t[i] = self$ker$dkdt( t[i],t(t))%*%self$b #*usgrad[il,i] } return(list("pred"=y_t,"grad"=z_t))#,"y_p"=y_p,"z_p"=z_p,"zt"=z_t,"zt2"=z_t2 ) ) }, predictT = function(testT) { ########### calculate optimised lengthscale ##################################################### y = scale(as.matrix(self$y), center=TRUE, scale=FALSE) mean_y= apply(as.matrix(self$y),2,mean) t = as.numeric(self$t) testn = length(testT) # ############################## make prediction for derivative ###################################################### y_t = array(c(0),testn) z_t = array(c(0),testn) for (i in 1:testn) { y_t[i] = self$ker$kern(t(t),testT[i])%*%self$b +mean_y z_t[i] = self$ker$dkdt( testT[i],t(t))%*%self$b } return(list("pred"=y_t,"grad"=z_t)) }, lossRK = function(par,tl1,y_d,jitter) { self$ker$k_par = exp(par) #self$ker$k_par = (par) res =0 tl1=as.matrix(tl1) n=max(dim(tl1)) K=matrix(c(0),nco=n, nrow=n) for (i in 1:n) { for (j in 1:n) { K[i,j]=self$ker$kern(tl1[i],tl1[j]) } } # calculate the coeffiecient 'b' of basis function ik = solve(K+jitter[1]*diag(n)) b=ik%*%y_d for(i in 1:n) { res=res+ (y_d[i]- self$ker$kern(t(tl1),tl1[i])%*%b)^2 } if(length(par)==1){ res = res + jitter* t(b)%*%K%*%b } else { res = res + jitter[1]* t(b)%*%K%*%b + jitter[2]*self$ker$k_par[2] } res }, grlossRK = function(par,tl1,y_d,jitter) { self$ker$k_par = exp(par) #self$ker$k_par = (par) tl1=as.matrix(tl1) n=max(dim(tl1)) np = length(par) dres = array(c(0),c(1,np)) K=matrix(c(0),nco=n, nrow=n) dK= array(c(0),c(n,n,np)) db=array(c(0),c(n,np)) for (i in 1:n) { for (j in 1:n) { #ts = tl1[i]-tl1[j] dK[i,j,1:np]=self$ker$dkd_kpar(tl1[i],tl1[j]) K[i,j]=self$ker$kern(tl1[i],tl1[j]) } } ik = solve(K+jitter[1]*diag(n)) b=ik%*%y_d for(k in 1:np){ db[,k]= -ik%*%dK[,,k]%*%ik%*%y_d } for(k in 1:np) { for(i in 1:n) { dres[,k]=dres[,k]+ 2* (y_d[i]- self$ker$kern(t(tl1),tl1[i])%*%b) *(-1)*( dK[,i,k]%*%b + K[,i]%*%db[,k] ) } if(k==1){ dres[,k] = dres[,k] + ( t(db[,k])%*%K%*%b + t(b)%*%dK[,,k]%*%b + t(b)%*%K%*%db[,k] ) * jitter[k] } else { dres[,k] = dres[,k] + ( t(db[,k])%*%K%*%b + t(b)%*%dK[,,k]%*%b + t(b)%*%K%*%db[,k] ) * jitter[k-1] + jitter[k]*self$ker$k_par[k] } } dres }, numgrad = function(par,tl1,y_d,jitter) { require(numDeriv) return( grad(self$loss11,par,,,,tl1,y_d,jitter) ) }, skcross = function( init,bounded ) { tbd = array(c(1), c(length(self$ker$k_par)) ) if( missing(bounded) ) { bounded= c(-Inf,Inf) lbound = bounded[1]*tbd ubound = bounded[2]*tbd } else { lbound = log(bounded[1]*tbd) ubound = log(bounded[2]*tbd) } y = scale(as.matrix(self$y), center=TRUE, scale=FALSE) mean_y= apply(as.matrix(self$y),2,mean) t_y = matrix(c(self$t),ncol=dim(y)[1]) n_y=max(dim(t_y)) np = length(self$ker$k_par) jitter=matrix(c(0),ncol=np) fold = 3 oneba = seq(1,n_y,3) twoba = seq(2,n_y,3) thrba = c(1:n_y)[c(-oneba,-twoba)] fd<-list(oneba,twoba,thrba) for(index in 1:np) { crlambda=c(5,1,0.1,0.01,1e-03,1e-04) #crlambda=c(100,10,5,1,0.1,0.01,1e-03,1e-04,1e-05,1e-06,1e-07) for(iii in 1:4) { n_l = length(crlambda) cres = matrix(c(0),ncol=n_l,nrow=1) ibreak=0 for (k in 1:n_l) { jitter[index] = crlambda[k] jitter = as.numeric(jitter) for (c in 1: fold) { t = t_y[1,-fd[[c]]] y_d= y[-fd[[c]],1] n = length(t) s1 = tryCatch({ optim( log(init),self$lossRK,self$grlossRK,t,y_d,jitter,method="L-BFGS-B",lower=lbound,upper=ubound ) #optim( log(self$ker$k_par),self$lossRK,self$grlossRK,t,y_d,jitter,method="L-BFGS-B",lower=lbound,upper=ubound ) }, warning = function(war) { print(paste("MY_WARNING: ",war)) }, error = function(err) { # error handler picks up where error was generated print(paste("MY_ERROR: ",err)) return( optim( log(init),self$lossRK,self$grlossRK,t,y_d,jitter,method="BFGS") ) # optim(log(self$ker$k_par),self$lossRK,self$grlossRK,t,y_d,jitter,method="BFGS") ) },finally = { } ) if(s1[[1]][1]==-100){ ibreak=100 break } self$ker$k_par = exp(s1$par) cK=matrix(c(0),ncol=n, nrow=n) # make new b for (i in 1:n) { for (j in 1:n) { cK[i,j]=self$ker$kern(t[i],t[j]) } } # calculate the coeffiecient 'b' of basis function cik = solve(cK+ jitter[1]*diag(n)) cb=cik%*%y_d for(i in 1:length(fd[[c]])) { cres[k] = cres[k]+ (y[fd[[c]],1][i]- self$ker$kern(t_y[1,fd[[c]]][i],t(t))%*%cb)^2 } } } ank = crlambda[which(cres==min(cres))] pank=ank/10 crlambda = c( ank-2*pank,ank-pank,ank+pank,ank+2*pank) } jitter[index] = ank } self$lambda = jitter s1 = optim( log(init),self$lossRK,self$grlossRK,as.numeric(t_y),as.numeric(y),jitter,method="L-BFGS-B",lower=lbound,upper=ubound ) self$ker$k_par = exp(s1$par) K1=matrix(c(0),ncol=n_y, nrow=n_y) for (i in 1:n_y) { for (j in 1:n_y) { K1[i,j]=self$ker$kern(t_y[i],t_y[j]) } } # calculate the coeffiecient 'b' of basis function ik1 = solve(K1 + self$lambda[1]*diag(n_y) ) self$b= as.numeric( ik1%*%as.numeric(y) ) return(jitter) }, mkcross = function(init) { y = scale(as.matrix(self$y), center=TRUE, scale=FALSE) #mean_y= apply(as.matrix(self$y),2,mean) t_y = matrix(c(self$t),ncol=dim(y)[1]) n_y=max(dim(t_y)) np = length(self$ker$k_par) jitter=matrix(c(0),ncol=np) fold = 3 oneba = seq(1,n_y,3) twoba = seq(2,n_y,3) thrba = c(1:n_y)[c(-oneba,-twoba)] fd<-list(oneba,twoba,thrba) crl1=c(100,10,5,1,0.3,0.1,0.01,1e-03,1e-04,1e-05,1e-06) crl2=c(100,10,5,1,0.1,0.01,1e-03,1e-04,1e-05,1e-06,1e-07) crll<-list(crl1,crl2) for(index in 1:np) { crlambda=crll[[index]] #crlambda=c(100,10,5,1,0.1,0.01,1e-03,1e-04,1e-05,1e-06,1e-07) for(iii in 1:4) { n_l = length(crlambda) cres = matrix(c(0),ncol=n_l,nrow=1) ibreak=0 for (k in 1:n_l) { jitter[index] = crlambda[k] jitter = as.numeric(jitter) for (c in 1: fold) { t = t_y[1,-fd[[c]]] y_d= y[-fd[[c]],1] n = length(t) s1 = optim( log(init),self$lossRK,self$grlossRK,t,y_d,jitter,method="L-BFGS-B",lower=log(c(0.001,0.001)),upper=log(c(1000,1000)) ) self$ker$k_par = exp(s1$par) cK=matrix(c(0),ncol=n, nrow=n) for (i in 1:n) { for (j in 1:n) { cK[i,j]=self$ker$kern(t[i],t[j]) } } # calculate the coeffiecient 'b' of basis function cik = solve(cK+ jitter[1]*diag(n)) cb=cik%*%y_d for(i in 1:length(fd[[c]])) { cres[k] = cres[k]+ (y[fd[[c]],1][i]- self$ker$kern(t_y[1,fd[[c]]][i],t(t))%*%cb)^2 } } } ank = crlambda[which(cres==min(cres))] pank=ank/10 crlambda = c( ank-2*pank,ank-pank,ank+pank,ank+2*pank) } jitter[index] = ank } self$lambda = jitter s1 = optim( log(init),self$lossRK,self$grlossRK,self$t,self$y,jitter,method="L-BFGS-B",lower=log(c(0.001,0.001)),upper=log(c(1000,1000)) ) self$ker$k_par = exp(s1$par) K1=matrix(c(0),ncol=n_y, nrow=n_y) for (i in 1:n_y) { for (j in 1:n_y) { K1[i,j]=self$ker$kern(self$t[i],self$t[j]) } } # calculate the coeffiecient 'b' of basis function ik1 = solve(K1 + self$lambda[1]*diag(n_y) ) self$b= as.numeric( ik1%*%as.numeric(y) ) return(jitter) }, loss11 = function(par,tl1,y_d,jitter) { self$ker$k_par = exp(par) #self$ker$k_par = (par) res =0 tl1=as.matrix(tl1) n=max(dim(tl1)) K=matrix(c(0),nco=n, nrow=n) for (i in 1:n) { for (j in 1:n) { K[i,j]=self$ker$kern(tl1[i],tl1[j]) } } # calculate the coeffiecient 'b' of basis function ik = solve(K+jitter[1]*diag(n)) b=ik%*%y_d for(i in 1:n) { res=res+ (y_d[i]- self$ker$kern(t(tl1),tl1[i])%*%b)^2 } if(length(par)==1){ res = res + jitter* t(b)%*%K%*%b } else { res = res + jitter[1]* t(b)%*%K%*%b + jitter[2]*self$ker$k_par[2] } res }, grloss11 = function(par,tl1,y_d,jitter) { self$ker$k_par = exp(par) #self$ker$k_par = (par) tl1=as.matrix(tl1) n=max(dim(tl1)) np = length(par) dres = array(c(0),c(1,np)) K=matrix(c(0),nco=n, nrow=n) dK= array(c(0),c(n,n,np)) db=array(c(0),c(n,np)) for (i in 1:n) { for (j in 1:n) { #ts = tl1[i]-tl1[j] dK[i,j,1:np]=self$ker$dkd_kpar(tl1[i],tl1[j]) K[i,j]=self$ker$kern(tl1[i],tl1[j]) } } ik = solve(K+jitter[1]*diag(n)) b=ik%*%y_d for(k in 1:np){ db[,k]= -ik%*%dK[,,k]%*%ik%*%y_d } for(k in 1:np) { for(i in 1:n) { dres[,k]=dres[,k]+ 2* (y_d[i]- self$ker$kern(t(tl1),tl1[i])%*%b) *(-1)*( dK[,i,k]%*%b + K[,i]%*%db[,k] ) } if(k==1){ dres[,k] = dres[,k] + ( t(db[,k])%*%K%*%b + t(b)%*%dK[,,k]%*%b + t(b)%*%K%*%db[,k] ) * jitter[k] } else { dres[,k] = dres[,k] + ( t(db[,k])%*%K%*%b + t(b)%*%dK[,,k]%*%b + t(b)%*%K%*%db[,k] ) * jitter[k-1] + jitter[k]*self$ker$k_par[k] } } dres } ) )