library(R6) library(pspline) Warp <- R6Class("Warp", public = list( y = NULL, t = NULL, b = NULL, lambda=NULL, ker= NULL, tw = 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("warp ",self$ker$greet(),".\n")) }, showker = function () { cat(paste0("ker is", self$ker$greet(), ".\n")) }, warpLoss = function( par,len,p0,eps ) { tor = as.matrix(self$t) n = max(dim(tor)) lam=par[1] self$ker$k_par=exp(len) lambda_t = self$lambda y_r = self$y self$b = par[2:(n+1)] kw=(2*pi/ (p0 + eps*tanh(lam) ) )^2 tw = matrix(c(0),ncol = n,nrow=1) for (i in 1:n) { tw[1,i] = self$ker$kern(tor[i,1],t(tor[,1]))%*%self$b^2 } t=t(tw) h=diff(t) z_t1 = matrix(c(0),ncol = n,nrow=1) z_t2 = matrix(c(0),ncol = n,nrow=2) y_t = matrix(c(0),ncol = n,nrow=2) ### calculate Euler 1st and 2nd order, if 1st order is enornouse do not do spline for(i in 2:(n-1) ) { z_t1[1,i] = (y_r[i+1]-y_r[i-1])/ (h[i]+h[i-1]) } z_t1[1,1]=(y_r[2]-y_r[1])/h[1] z_t1[1,n]=(y_r[n]-y_r[n-1])/h[n-1] for(i in 2:(n-1) ) { z_t2[1,i] = (z_t1[i+1]-z_t1[i-1])/ (h[i]+h[i-1]) } z_t2[1,1]=(z_t1[2]-z_t1[1])/h[1] z_t2[1,n]=(z_t1[n]-z_t1[n-1])/h[n-1] y_t[1,]= y_r res= sum( (z_t2[1,] + kw*y_t[1,])^2 ) + lambda_t*( (tor[1,1]- t[1,1])^2+ (tor[n,1]- t[n,1])^2 ) if( is.nan( sum(z_t1) ) ) { print(len) return(res= 100000) } if( sum(abs(z_t1) )<1000 ) { z_t2[1,] = tryCatch( { predict(sm.spline(t, y_r,df=n-5),t, 2) }, warning = function(war) { print(paste("MY_WARNING: ",war)) }, error = function(err) { # error handler picks up where error was generated print(paste("warp_ERROR: ",err,i)) return(list(-10000,0) ) #return(optim(par1,lossWarpTruelenFtanh,gr=NULL,pick,lambda_t,y_use[1,],len,p0,eps,method="L-BFGS-B") ) },finally = { } ) if(z_t2[[1]][1]==-10000){ return(res= 100000) } y_t[1,] = predict(sm.spline(t, y_r,df=n-5),t, 0) res= sum( (z_t2[1,] + kw*y_t[1,])^2 ) + lambda_t*( (tor[1,1]- t[1,1])^2+ (tor[n,1]- t[n,1])^2 ) } res }, warpLossLen=function(par,lam,p0,eps ){ tor = as.matrix(self$t) n = max(dim(tor)) lam = lam#par[1] self$ker$k_par=exp(par) lambda_t = self$lambda y_r = self$y #bw1= self$b kw=(2*pi/ (p0 + eps*tanh(lam) ) )^2 tw = matrix(c(0),ncol = n,nrow=1) for (i in 1:n) { tw[1,i] = self$ker$kern(tor[i,1],t(tor[,1]))%*%self$b^2 } t=t(tw) h=diff(t) z_t1 = matrix(c(0),ncol = n,nrow=1) z_t2 = matrix(c(0),ncol = n,nrow=2) y_t = matrix(c(0),ncol = n,nrow=2) ### calculate Euler 1st and 2nd order, if 1st order is enornouse do not do spline for(i in 2:(n-1) ) { z_t1[1,i] = (y_r[i+1]-y_r[i-1])/ (h[i]+h[i-1]) } z_t1[1,1]=(y_r[2]-y_r[1])/h[1] z_t1[1,n]=(y_r[n]-y_r[n-1])/h[n-1] for(i in 2:(n-1) ) { z_t2[1,i] = (z_t1[i+1]-z_t1[i-1])/ (h[i]+h[i-1]) } z_t2[1,1]=(z_t1[2]-z_t1[1])/h[1] z_t2[1,n]=(z_t1[n]-z_t1[n-1])/h[n-1] y_t[1,]= y_r res= sum( (z_t2[1,] + kw*y_t[1,])^2 ) + lambda_t*( (tor[1,1]- t[1,1])^2+ (tor[n,1]- t[n,1])^2 ) if( is.nan( sum(z_t1) ) ) { print(len) return(res= 100000) } if( sum(abs(z_t1) )<1000 ){ z_t2[1,] = predict(sm.spline(t, y_r),t, 2) y_t[1,] = predict(sm.spline(t, y_r),t, 0) res= sum( (z_t2[1,] + kw*y_t[1,])^2 ) + lambda_t*( (tor[1,1]- t[1,1])^2+ (tor[n,1]- t[n,1])^2 ) } res }, warpSin = function( len ,lop ) { kw1= 1 #lop=3 lout = length(self$t) bw = self$b for(i in 1:lop) { par<-c(kw1,bw) test <- tryCatch({ optim(par,self$warpLoss,gr=NULL,len,p0,eps,method="BFGS") }, warning = function(war) { print(paste("MY_WARNING: ",war)) }, error = function(err) { # error handler picks up where error was generated print(paste("warp_ERROR: ",err,i)) return(list(-10000,0) ) #return(optim(par1,lossWarpTruelenFtanh,gr=NULL,pick,lambda_t,y_c,len,p0,eps,method="L-BFGS-B") ) },finally = { } ) if(test[[1]][1]==-10000) { wscore=100000 break }else { self$b= test$par[2:(lout+1)] kw1 = test$par[1] par2=len bw= self$b #bbbl$warpLossLen(par2,kw1,p0,eps) test2<-tryCatch({ optim(par2,self$warpLossLen,gr=NULL,kw1,p0,eps,method="L-BFGS-B") }, warning = function(war) { print(paste("MY_WARNING: ",war)) }, error = function(err) { # error handler picks up where error was generated print(paste("warp_ERROR: ",err,i)) return(list(-10000,0) ) #return(optim(par1,lossWarpTruelenFtanh,gr=NULL,pick,lambda_t,y_c,len,p0,eps,method="L-BFGS-B") ) },finally = { } ) if(test2[[1]][1]==-10000){ wscore= 100000 break } wscore=test2$value len=test2$par } } self$ker$k_par = exp(len) n = length(self$t) tw = matrix(c(0),ncol = n,nrow=1) for (i in 1:n) { tw[1,i] = self$ker$kern(self$t[i],t(self$t))%*%self$b^2 } self$tw = tw return( list( "rescore"=wscore,"len"=len) ) }, slowWarp = function(str,end,p0,eps) { bw = rep(c(1),length(self$t)) kw = 1 loscore=c(0) lens=seq(str,end,0.1) #seq(4.7,6.5,0.1) for(iii in 1:length(lens)) { par<-c(kw,bw) len=lens[iii] ptm <- proc.time() test<-optim(par,bbb$warpLoss,gr=NULL,len,p0,eps,method="BFGS")#optim(par,lossWarpTruelenFtanh,gr=NULL,pick,lambda_t,y_c,len,p0,eps,method="BFGS") proc.time() - ptm loscore[iii]= test$value } len=lens[which(loscore==min(loscore) )][1] return(list("len"=len,"loscore"=loscore)) } ) ) lossWarpTruelenFtanh<-function( par,times,lambda_t,y_r,l,p0,eps) { tor = as.matrix(times) n = max(dim(tor)) lam=par[1] len=exp(l) bw1 = par[2:(n+1)] sigmd = Sigmoid$new(len) kw=(2*pi/ (p0 + eps*tanh(lam) ) )^2 tw = matrix(c(0),ncol = n,nrow=1) for (i in 1:n) { tw[1,i] = sigmd$kern(tor[i,1],t(tor[,1]))%*%bw1^2 } t=t(tw) h=diff(t) z_t1 = matrix(c(0),ncol = n,nrow=1) z_t2 = matrix(c(0),ncol = n,nrow=2) y_t = matrix(c(0),ncol = n,nrow=2) ### calculate Euler 1st and 2nd order, if 1st order is enornouse do not do spline for(i in 2:(n-1) ) { z_t1[1,i] = (y_r[i+1]-y_r[i-1])/ (h[i]+h[i-1]) } z_t1[1,1]=(y_r[2]-y_r[1])/h[1] z_t1[1,n]=(y_r[n]-y_r[n-1])/h[n-1] for(i in 2:(n-1) ) { z_t2[1,i] = (z_t1[i+1]-z_t1[i-1])/ (h[i]+h[i-1]) } z_t2[1,1]=(z_t1[2]-z_t1[1])/h[1] z_t2[1,n]=(z_t1[n]-z_t1[n-1])/h[n-1] y_t[1,]= y_r res= sum( (z_t2[1,] + kw*y_t[1,])^2 ) + lambda_t*( (tor[1,1]- t[1,1])^2+ (tor[n,1]- t[n,1])^2 ) if( is.nan( sum(z_t1) ) ) { print(len) return(res= 100000) } if( sum(abs(z_t1) )<1000 ){ z_t2[1,] = tryCatch({ #predict(sm.spline(t, y_r),t, 2) predict(sm.spline(t, y_r,df=n-5),t, 2) }, warning = function(war) { print(paste("MY_WARNING: ",war)) }, error = function(err) { # error handler picks up where error was generated print(paste("warp_ERROR: ",err,i)) return(list(-10000,0) ) #return(optim(par1,lossWarpTruelenFtanh,gr=NULL,pick,lambda_t,y_use[1,],len,p0,eps,method="L-BFGS-B") ) },finally = { } ) if(z_t2[[1]][1]==-10000){ return(res= 100000) } #y_t[1,] = predict(sm.spline(t, y_r),t, 0) y_t[1,] = predict(sm.spline(t, y_r,df=n-5),t, 0) res= sum( (z_t2[1,] + kw*y_t[1,])^2 ) + lambda_t*( (tor[1,1]- t[1,1])^2+ (tor[n,1]- t[n,1])^2 ) } res }