library(R6)
library(deSolve)
library(pracma)
library(mvtnorm)


ode<-R6Class("ode",
  public = list(
    sample = NULL,
    t = NULL,
    y_ode = NULL,
    ode_par=NULL,
    ode_fun=NULL,
    gr_lNODE=NULL,

    initialize = function(sample = NULL, fun = NULL , grfun=NULL ,ode_par = NULL) {
      self$sample <- sample
      self$ode_par <- ode_par
      self$ode_fun = fun
      self$gr_lNODE = grfun
      self$greet()
    },

    greet = function() {
      cat(paste0("ode is sample ", self$sample, ".\n"))
    },
## ordinary differential equation in matrix form
## ode solver
    solve_ode = function(par_ode,xinit,tinterv){
        x0 = xinit#as.matrix(c(0.5,1))
		t0 <- tinterv[1]#0; 
		tf <- tinterv[2]#6
		solt <- ode23s(self$ode_fun, t0, tf, x0, par_ode=par_ode,hmax=(tf-t0) )
		pick = seq(1,length(solt$t),self$sample)
		self$t = solt$t[pick]
		rownames(solt$y) <- NULL
		self$y_ode= t(solt$y[pick,])
    },
    
    rmsfun = function(par_ode,state,M1,true_par){
        out <- ode(y =state,times =self$t, func = M1,parms=par_ode,method="ode23")
        funlos = sqrt( sum( ( t(out[,-1]) - self$y_ode)^2 ) / length(c(self$y_ode)) )
        parlos= sqrt( sum( (par_ode - true_par)^2 ) / length(true_par)  )
        relalos= sqrt( sum( ( (par_ode - true_par)/true_par )^2 ) / length(true_par)  )
        c(funlos,parlos,relalos)
    },
## gradient produced by ode
	gradient = function(y_p,par_ode){
	    ## y_p need to be 50*2  apply to each row
        ydif = apply(y_p,2,self$ode_fun,t=self$t,par_ode)
		ydif
	},
## mismatch between ode gradient and interpolant gradient
	##################################### griadient mathcing ##################
	lossNODE= function( par,y_p,z_p) {
		par_ode= exp(par)
        grad_ode = self$gradient(y_p,par_ode)
		lm = max(dim(y_p)) 
		res=0
		start = 2
		res = sum( (z_p[,start:lm]-grad_ode[,start:lm])^2 )
		#res = sum( (z_p-grad_ode)^2 )
		res
	},

	grlNODE= function( par,y_p,z_p) { 
		par_ode= exp(par)   ## ode par
		lm = max(dim(y_p))
        grad_ode = self$gradient(y_p,par_ode)
        start = 2
		dres=self$gr_lNODE( par_ode,grad_ode[,start:lm],y_p[,start:lm],z_p[,start:lm] )
		#dres=self$gr_lNODE( par_ode,grad_ode,y_p,z_p )
		dres
	},
    
    loss32NODE= function( par,y_p,z_p) {
		par_ode= exp(par)
        grad_ode = self$gradient(y_p,par_ode)
		lm = max(dim(y_p)) 
		res=0
		res = sum( (z_p[,]-grad_ode[,])^2 )
		res
	},

	grl32NODE= function( par,y_p,z_p) { 
		par_ode= exp(par)   ## ode par
		lm = max(dim(y_p))
        grad_ode = self$gradient(y_p,par_ode)
		#dres= c(0)
		dres=self$gr_lNODE( par_ode,grad_ode,y_p,z_p )
		dres
	},


    optim_par = function(par,y_p,z_p){
		op = optim(log(par),self$lossNODE,self$grlNODE,y_p,z_p,method="L-BFGS-B")
		self$ode_par = exp(op$par)
		exp(op$par)
	}

  )

)