R/GQD.mcmc.R

In DiffusionRgqd: Inference and Analysis for Generalized Quadratic Diffusions

globalVariables('priors')
GQD.mcmc <-
function(X,time,mesh=10,theta,sds,updates,burns=min(round(updates/2),25000),Dtype='Saddle',Trunc=c(4,4),RK.order=4,P=200,alpha=0,lower=min(na.omit(X))/2,upper=max(na.omit(X))*2,exclude=NULL,plot.chain=TRUE,Tag=NA,wrt=FALSE,print.output=TRUE,palette='mono')
{
  solver   =function(Xs, Xt, theta, N , delt , N2, tt  , P , alpha, lower , upper, tro  ){}
  rm(list =c('solver'))
   theta = theta+runif(length(theta),0.01,0.02)*sign(theta)
   adapt=0
   check_for_model=function()
  {
    txt=''
    namess=c('G0','G1','G2','Q0','Q1','Q2')
    func.list=rep(0,length(namess))
    obs=objects(pos=1)
    for(i in 1:length(namess))
    {
    if(sum(obs==namess[i])){func.list[i]=1}
    }

    check=F
    if(sum(func.list)==0)
   {
   txt='
  --------------------------------------------------------------------------------
  No model has been defined yet! Try for example:
  --------------------------------------------------------------------------------
  GQD.remove()
  G0=function(t){theta[1]*theta[2]}
  G1=function(t){-theta[1]}
  Q1=function(t){theta[3]*theta[3]}

  model=GQD.mcmc(X,time,10,theta =rep(1,3),sds=rep(0.1,3),updates=10000)
  --------------------------------------------------------------------------------
   '
   check=T
   }
   if((sum(func.list)>0)&&(sum(func.list[-c(1:3)])==0))
   {
   txt='
  --------------------------------------------------------------------------------
  At least one diffusion coefficient has to be defined! Try for example:
  --------------------------------------------------------------------------------
  GQD.remove()
  G0=function(t){theta[1]*theta[1]}
  model=GQD.mcmc(X,time,10,theta =c(0.1),sds=0.1,updates=10000)
  --------------------------------------------------------------------------------
   '
     check=T
    }
    return(list(check=check,txt=txt))
  }
  check_for=check_for_model()
  if(check_for[[1]]){stop(check_for[[2]])}
  theta = theta+runif(length(theta),0.001,0.002)*sign(theta)
  pow=function(x,p)
  {
    x^p
  }
  prod=function(a,b){a*b}
  T.seq=time
  # Warning Module

   TR.order=Trunc[1]
   DTR.order=Trunc[2]
   Dtypes =c('Saddle','Normal','Gamma','InvGamma','Beta')
   Dindex = which(Dtypes==Dtype)
   IntRange = c(lower,upper)
   IntDelta =1/100
   Xtr = min(IntRange)
    b1 = '\n==============================================================================\n'
   b2 = '==============================================================================\n'
  warn=c(
    '1.  Missing input: Argument {X} is missing.\n'
    ,'2.  Missing input: Argument {time} is missing.\n'
    ,'3.  Missing input: Argument {theta} is missing.\n'
    ,'4.  Missing input: Argument {sds} is missing.\n'
    ,'5.  Input type: Argument {X} must be of type vector!.\n'
    ,'6.  Input type: Argument {time} must be of type vector!.\n'
    ,'7.  Input: Less starting parameters than model parameters.\n'
    ,'8.  Input: More starting parameters than model parameters.\n'
    ,'9.  Input: length(X) must be > 10.\n'
    ,'10. Input: length(time) must be > 10.\n'
    ,'11. Input: length(lower)!=1.\n'
    ,'12. Input: length(upper)!=1.\n'
    ,'13. Input: length(P)!=1.\n'
    ,'14. Input: length(mesh)!=1.\n'
    ,'15. Input: length(alpha)!=1.\n'
    ,'16. Input: length(Trunc)!=1.\n'
    ,'17. Input: length(RK.order)!=1.\n'
    ,'18. Density: Dtype has to be one of Saddle, Normal, Gamma, InvGamma or Beta.\n'
    ,'19. Density: Range [lower,upper] must be strictly positive for Dtype Gamma or InvGamma.\n'
    ,'20. Density: Dtype cannot be Beta for observations not in (0,1).\n'
    ,'21. Density: Argument {upper} must be > {lower}.\n'
    ,'22. Density: P must be >= 10.\n'
    ,'23. Density: Trunc[2] must be <= Trunc[1].\n'
    ,'24. ODEs : Large max(diff(time))/mesh may result in poor approximations. Try larger mesh.\n'
    ,'25. ODEs : max(diff(time))/mesh must be <1.\n'
    ,'26. ODEs : Runge Kutta scheme must be of order 4 or 10.\n'
    ,'27. ODEs : Argument {mesh} must be >= 5.\n'
    ,'28. Input: length(X)!=length(time).\n'
    ,'29. MCMC : Argument {burns} must be < {updates}.\n'
    ,'30. MCMC : Argument {updates} must be > 2.\n'
    ,'31. MCMC : length(theta)!=length(sds).\n'
    ,'32. Model: There has to be at least one model coefficient.\n'
    ,'33. Input: length(updates)!=1.\n'
    ,'34. Input: length(burns)!=1.\n'
    ,'35. Prior: priors(theta) must return a single value.\n'
    ,'36. Input: NAs not allowed.\n'
    ,'37. Input: length(Dtype)!=1.\n'
    ,'38. Input: NAs not allowed.\n'
    ,'39. Input: Time series contains values of small magnitude.\n    This may result in numerical instabilities.\n    It may be advisable to scale the data by a constant factor.\n'
  )

   warntrue = rep(F,40)
      check.thetas = function(theta,tt)
   {
     t=tt
     theta = theta+runif(length(theta),0.001,0.002)*sign(theta)
     namess=c('G0','G1','G2','Q0','Q1','Q2')
     func.list=rep(0,length(namess))
     obs=objects(pos=1)
     for(i in 1:length(namess))
     {
       if(sum(obs==namess[i])){func.list[i]=1}
     }
     pers.represented = rep(0,length(theta))
     for(i in which(func.list==1))
     {
       for(j in 1:length(theta))
       {
         dresult1=eval(body(namess[i]))
         theta[j] = theta[j]+runif(1,0.1,0.2)
         dresult2=eval(body(namess[i]))
         dff = abs(dresult1-dresult2)
         if(any(round(dff,6)!=0)){pers.represented[j]=pers.represented[j]+1}
       }
     }
     return(pers.represented)
   }

   check.thetas2 = function(theta)
   {
     namess=c('G0','G1','G2','Q0','Q1','Q2')
     func.list=rep(0,length(namess))
     obs=objects(pos=1)
     for(i in 1:length(namess))
     {
       if(sum(obs==namess[i])){func.list[i]=1}
     }

     l=0
     for(k in which(func.list==1))
     {
       str=body(namess[k])[2]
       for(i in 1:length(theta))
       {
         for(j in 1:20)
         {
           str=sub(paste0('theta\\[',i,'\\]'),'clear',str)
         }
       }
       l=l+length(grep('theta',str))
       l
     }
     return(l)
   }

   # Check missing values first:
   if(missing(X))                                                {warntrue[1]=T}
   if(missing(time))                                             {warntrue[2]=T}
   if(missing(theta))                                            {warntrue[3]=T}
   if(missing(sds))                                              {warntrue[4]=T}
   if(!is.vector(X))                                             {warntrue[5]=T}
   if(!is.vector(time))                                          {warntrue[6]=T}
   # Check model parameters:
   if(check.thetas2(theta)!=0)                                   {warntrue[7]=T}
   if(!warntrue[7]){if(any(check.thetas(theta,T.seq)==0))        {warntrue[8]=T}}

   # Check input length:
   if(length(X)<10)                                              {warntrue[9]=T}
   if(length(time)<10)                                          {warntrue[10]=T}
   if(length(lower)>1)                                          {warntrue[11]=T}
   if(length(upper)>1)                                          {warntrue[12]=T}
   if(length(P)!=1)                                             {warntrue[13]=T}
   if(length(mesh)!=1)                                        {warntrue[14]=T}
   if(length(alpha)!=1)                                         {warntrue[15]=T}
   if(length(Trunc)!=2)                                         {warntrue[16]=T}
   if(length(RK.order)!=1)                                      {warntrue[17]=T}
   if(length(updates)!=1)                                     {warntrue[33]=T}
   if(length(burns)!=1)                                         {warntrue[34]=T}
   if(length(Dtype)!=1)                                         {warntrue[37]=T}


   # Check density approx parameters:
   if(sum(Dindex)==0)                                          {warntrue[18] =T}
   if(!warntrue[18])
   {
    if((Dindex==3)|(Dindex==4)){if(lower[1]<=0)                {warntrue[19] =T}}
    if(Dindex==5){if(any(X<=0)|any(X>=1))                      {warntrue[20] =T}}
   }
   if(!any(warntrue[c(11,12)])){if(upper<=lower)                {warntrue[21] =T}}
   if(!warntrue[13]){if(P<10)                                   {warntrue[22] =T}}
   if(!warntrue[16]){if(Trunc[2]>Trunc[1])                      {warntrue[23] =T}}

   #  Miscelaneous checks:
   excl=0
   if(is.null(exclude)){excl=length(T.seq)-1+200}
   if(!is.null(exclude)){excl=exclude}
   test.this =max(diff(T.seq)[-excl])/mesh
   if(test.this>0.1)                                            {warntrue[24]=T}
   if(test.this>=1)                                             {warntrue[25]=T}
   if(!warntrue[17]){if(!((RK.order==4)|(RK.order==10)))        {warntrue[26]=T}}
   if(!warntrue[14]){if(mesh<5)                               {warntrue[27]=T}}
   if(length(X)!=length(time))                                  {warntrue[28]=T}
   if(!any(warntrue[c(33,34)])){if(burns>updates)             {warntrue[29]=T}}
   if(!warntrue[33]){if(updates<2)                            {warntrue[30]=T}}
   if(length(theta)!=length(sds))                               {warntrue[31]=T}
   if(any(is.na(X))||any(is.na(time)))                          {warntrue[36]=T}

   # Print output:
   if(any(warntrue))
   {
      prnt = b1
      for(i in which(warntrue))
      {
         prnt = paste0(prnt,warn[i])
      }
      prnt = paste0(prnt,b2)
      stop(prnt)
   }
     if(any(X<10^-2)){warntrue[39]=T}
   # Print warnings:
   if(any(warntrue))
   {
      prnt = b1
      for(i in which(warntrue))
      {
         prnt = paste0(prnt,warn[i])
      }
      prnt = paste0(prnt,b2)
      warning(prnt)
   }

   nnn=length(X)

   #==============================================================================
   #                           Data resolution Module
   #==============================================================================

   homo=T
   homo.res=T
   delt=(diff(T.seq)/mesh)[1]
   t=T.seq
   if(is.null(exclude))
   {
     if(sum(round(diff(T.seq)-diff(T.seq)[1],10)==0)!=length(T.seq)-1){homo.res=F}
   }
   if(!is.null(exclude))
   {
     if(sum(round(diff(T.seq)[-excl]-c(diff(T.seq)[-excl])[1],10)==0)!=length(T.seq)-1-length(excl)){homo.res=F}
   }
   #==============================================================================
   #                           Prior distribution Module
   #==============================================================================
   if(sum(objects(pos=1)=='priors')==1)
   {
      pp=function(theta){}
      body(pp)=parse(text =body(priors)[2])
      prior.list=paste0('d(theta)',':',paste0(body(priors)[2]))
      if(length(priors(theta))!=1){stop(" ==============================================================================
   Incorrect input: Prior distribution must return a single value only!
   ==============================================================================");}
   }
   if(sum(objects(pos=1)=='priors')==0)
   {
     prior.list=paste0('d(theta)',':',' None.')
     pp=function(theta){1}
   }
   
   
   #==============================================================================
   #                           Solution TYPE Module
   #==============================================================================
   namess=c('G0','G1','G2','Q0','Q1','Q2')
   func.list=rep(0,6)
   obs=objects(pos=1)
   for(i in 1:6)
   {
     if(sum(obs==namess[i])){func.list[i]=1}
   }

   state1=(sum(func.list[c(3,5,6)]==1)==0)
   if(state1){DTR.order=2;TR.order=2;sol.state='Normally distributed diffusion.';}
   if((state1&(Dtype!='Saddle'))){TR.order=2;DTR.order=2;sol.state='2nd Ord. Truncation + Std Normal Dist.';}
   state2=!state1
   if(state2)
   {
     state2.types=c('Saddlepoint Appr.   ',
                  ' Ext. Normal Appr.    ',
                  ' Ext. Gamma Appr.     ',
                  ' Ext. Inv. Gamma Appr.')
     sol.state=paste0(TR.order,' Ord. Truncation +',DTR.order,'th Ord. ',state2.types[Dindex])
   }



if(TR.order==2)
{
    fpart=
    '
    #include <RcppArmadillo.h>
    #include <math.h>
    #include <Rcpp.h>
    #define pi           3.14159265358979323846  /* pi */
    using namespace arma;
    using namespace Rcpp;
    using namespace R;
    // [[Rcpp::depends("RcppArmadillo")]]
    // [[Rcpp::export]]
    vec prod(vec a,vec b)
    {
     return(a%b);
    }

    mat f(mat a,vec theta,vec t,int N2)
    {

       mat atemp(N2,2);'
}
if(TR.order==4)
{
    fpart=
    '#include <RcppArmadillo.h>
    #include <RcppArmadillo.h>
    #include <math.h>
    #include <Rcpp.h>
    #define pi           3.14159265358979323846  /* pi */
    using namespace arma;
    using namespace Rcpp;
    using namespace R;

    // [[Rcpp::depends("RcppArmadillo")]]
    // [[Rcpp::export]]
    vec prod(vec a,vec b)
    {
     return(a%b);
    }

    mat f(mat a,vec theta,vec t,int N2)
    {

       mat atemp(N2,4);'
}
if(TR.order==6)
{
    fpart=
    '
    #include <RcppArmadillo.h>
    #include <math.h>
    #include <Rcpp.h>
    #define pi           3.14159265358979323846  /* pi */
    using namespace arma;
    using namespace Rcpp;
    using namespace R;

    // [[Rcpp::depends("RcppArmadillo")]]
    // [[Rcpp::export]]
    vec prod(vec a,vec b)
    {
     return(a%b);
    }

    mat f(mat a,vec theta,vec t,int N2)
    {

       mat atemp(N2,6);'
}
if(TR.order==8)
{
    fpart=
    '
    #include <RcppArmadillo.h>
    #include <math.h>
    #include <Rcpp.h>
    #define pi           3.14159265358979323846  /* pi */
    using namespace arma;
    using namespace Rcpp;
    using namespace R;
    // [[Rcpp::depends("RcppArmadillo")]]
    // [[Rcpp::export]]
    vec prod(vec a,vec b)
    {
     return(a%b);
    }

    mat f(mat a,vec theta,vec t,int N2)
    {

       mat atemp(N2,8);'
}
if(RK.order==10)
{
if(homo.res)
{
ODEpart= '
   return atemp;
}

// [[Rcpp::export]]
List  solver(vec Xs,vec Xt,vec theta,int N,double delt,int N2,vec tt,int P,double alpha,double lower,double upper,int tro)
{
    mat x0(N2,tro);
    mat xa(N2,tro);
    mat xe(N2,tro);
  mat fx0(N2,tro);
  mat fx1(N2,tro);
  mat fx2(N2,tro);
  mat fx3(N2,tro);
  mat fx4(N2,tro);
  mat fx5(N2,tro);
  mat fx6(N2,tro);
  mat fx7(N2,tro);
  mat fx8(N2,tro);
  mat fx9(N2,tro);
  mat fx10(N2,tro);
  mat fx11(N2,tro);
  mat fx12(N2,tro);
  mat fx13(N2,tro);
  mat fx14(N2,tro);
  mat fx15(N2,tro);
  mat fx16(N2,tro);

    double whch =0;
    x0.fill(0);
    x0.col(0)=Xs;
    vec d=tt;
    for (int i = 1; i < N+1; i++)
{
  fx0=f(x0,theta,d,N2)*delt;
  fx1=f(x0+0.1*fx0,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt,N2)*delt;
  fx2=f(x0+-0.915176561375291*fx0+1.45453440217827*fx1,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt,N2)*delt;
  fx3=f(x0+0.202259190301118*fx0+0.606777570903354*fx2,theta,d+0.809036761204472681298727796821953655285910174551513523258250*delt,N2)*delt;
  fx4=f(x0+0.184024714708644*fx0+0.197966831227192*fx2-0.0729547847313633*fx3,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt,N2)*delt;
  fx5=f(x0+0.0879007340206681*fx0+0.410459702520261*fx3+0.482713753678866*fx4,theta,d+0.981074190219795268254879548310562080489056746118724882027805*delt,N2)*delt;
  fx6=f(x0+0.085970050490246*fx0+0.330885963040722*fx3+0.48966295730945*fx4-0.0731856375070851*fx5,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt,N2)*delt;
  fx7=f(x0+0.120930449125334*fx0+0.260124675758296*fx4+0.0325402621549091*fx5-0.0595780211817361*fx6,theta,d+0.354017365856802376329264185948796742115824053807373968324184*delt,N2)*delt;
  fx8=f(x0+0.110854379580391*fx0-0.0605761488255006*fx5+0.321763705601778*fx6+0.510485725608063*fx7,theta,d+0.882527661964732346425501486979669075182867844268052119663791*delt,N2)*delt;
  fx9=f(x0+0.112054414752879*fx0-0.144942775902866*fx5-0.333269719096257*fx6+0.49926922955688*fx7+0.509504608929686*fx8,theta,d+0.642615758240322548157075497020439535959501736363212695909875*delt,N2)*delt;
  fx10=f(x0+0.113976783964186*fx0-0.0768813364203357*fx5+0.239527360324391*fx6+0.397774662368095*fx7+0.0107558956873607*fx8-0.327769124164019*fx9,theta,d+0.357384241759677451842924502979560464040498263636787304090125*delt,N2)*delt;
  fx11=f(x0+0.0798314528280196*fx0-0.0520329686800603*fx5-0.0576954146168549*fx6+0.194781915712104*fx7+0.145384923188325*fx8-0.0782942710351671*fx9-0.114503299361099*fx10,theta,d+0.117472338035267653574498513020330924817132155731947880336209*delt,N2)*delt;
  fx12=f(x0+0.985115610164857*fx0+0.330885963040722*fx3+0.48966295730945*fx4-1.37896486574844*fx5-0.861164195027636*fx6+5.78428813637537*fx7+3.28807761985104*fx8-2.38633905093136*fx9-3.25479342483644*fx10-2.16343541686423*fx11,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt,N2)*delt;
  fx13=f(x0+0.895080295771633*fx0+0.197966831227192*fx2-0.0729547847313633*fx3-0.851236239662008*fx5+0.398320112318533*fx6+3.63937263181036*fx7+1.5482287703983*fx8-2.12221714704054*fx9-1.58350398545326*fx10-1.71561608285936*fx11-0.0244036405750127*fx12,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt,N2)*delt;
  fx14=f(x0+-0.915176561375291*fx0+1.45453440217827*fx1+0*fx2+0*fx3-0.777333643644968*fx4+0*fx5-0.0910895662155176*fx6+0.0910895662155176*fx12+0.777333643644968*fx13,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt,N2)*delt;
  fx15=f(x0+0.1*fx0-0.157178665799771*fx2+0.157178665799771*fx14,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt,N2)*delt;
  fx16=f(x0+0.181781300700095*fx0+0.675*fx1+0.34275815984719*fx2+0*fx3+0.259111214548323*fx4-0.358278966717952*fx5-1.04594895940883*fx6+0.930327845415627*fx7+1.77950959431708*fx8+0.1*fx9-0.282547569539044*fx10-0.159327350119973*fx11-0.145515894647002*fx12-0.259111214548323*fx13-0.34275815984719*fx14-0.675*fx15,theta,d+delt,N2)*delt;

  x0=x0+(0.0333333333333333333333333333333333333333333333333333333333333*fx0
        +0.0250000000000000000000000000000000000000000000000000000000000*fx1
        +0.0333333333333333333333333333333333333333333333333333333333333*fx2
        +0.000000000000000000000000000000000000000000000000000000000000*fx3
        +0.0500000000000000000000000000000000000000000000000000000000000*fx4
        +0.000000000000000000000000000000000000000000000000000000000000*fx5
        +0.0400000000000000000000000000000000000000000000000000000000000*fx6
        +0.000000000000000000000000000000000000000000000000000000000000*fx7
        +0.189237478148923490158306404106012326238162346948625830327194*fx8
        +0.277429188517743176508360262560654340428504319718040836339472*fx9
        +0.277429188517743176508360262560654340428504319718040836339472*fx10
        +0.189237478148923490158306404106012326238162346948625830327194*fx11
        -0.0400000000000000000000000000000000000000000000000000000000000*fx12
        -0.0500000000000000000000000000000000000000000000000000000000000*fx13
        -0.0333333333333333333333333333333333333333333333333333333333333*fx14
        -0.0250000000000000000000000000000000000000000000000000000000000*fx15
        +0.0333333333333333333333333333333333333333333333333333333333333*fx16);
    xe  =  abs(fx1.col(1)-fx15.col(1))/360.0;
    if(xe.max()>whch)
    {
      whch = xe.max();
    }
    d=d+delt;
}
  '
}
if(!homo.res)
{
  ODEpart= '
   return atemp;
}

// [[Rcpp::export]]
List  solver(vec Xs,vec Xt,vec theta,int N,double delt,int N2,vec tt,int P,double alpha,double lower,double upper,int tro)
{
    mat x0(N2,tro);
    mat xa(N2,tro);
    mat xe(N2,tro);
  mat fx0(N2,tro);
  mat fx1(N2,tro);
  mat fx2(N2,tro);
  mat fx3(N2,tro);
  mat fx4(N2,tro);
  mat fx5(N2,tro);
  mat fx6(N2,tro);
  mat fx7(N2,tro);
  mat fx8(N2,tro);
  mat fx9(N2,tro);
  mat fx10(N2,tro);
  mat fx11(N2,tro);
  mat fx12(N2,tro);
  mat fx13(N2,tro);
  mat fx14(N2,tro);
  mat fx15(N2,tro);
  mat fx16(N2,tro);

    double whch =0;
    x0.fill(0);
    x0.col(0)=Xs;
    vec d=tt;
    for (int i = 1; i < N+1; i++)
{
  fx0=f(x0,theta,d,N2)%delt;
  fx1=f(x0+0.1*fx0,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt.col(0),N2)%delt;
  fx2=f(x0+-0.915176561375291*fx0+1.45453440217827*fx1,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt.col(0),N2)%delt;
  fx3=f(x0+0.202259190301118*fx0+0.606777570903354*fx2,theta,d+0.809036761204472681298727796821953655285910174551513523258250*delt.col(0),N2)%delt;
  fx4=f(x0+0.184024714708644*fx0+0.197966831227192*fx2-0.0729547847313633*fx3,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt.col(0),N2)%delt;
  fx5=f(x0+0.0879007340206681*fx0+0.410459702520261*fx3+0.482713753678866*fx4,theta,d+0.981074190219795268254879548310562080489056746118724882027805*delt.col(0),N2)%delt;
  fx6=f(x0+0.085970050490246*fx0+0.330885963040722*fx3+0.48966295730945*fx4-0.0731856375070851*fx5,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt.col(0),N2)%delt;
  fx7=f(x0+0.120930449125334*fx0+0.260124675758296*fx4+0.0325402621549091*fx5-0.0595780211817361*fx6,theta,d+0.354017365856802376329264185948796742115824053807373968324184*delt.col(0),N2)%delt;
  fx8=f(x0+0.110854379580391*fx0-0.0605761488255006*fx5+0.321763705601778*fx6+0.510485725608063*fx7,theta,d+0.882527661964732346425501486979669075182867844268052119663791*delt.col(0),N2)%delt;
  fx9=f(x0+0.112054414752879*fx0-0.144942775902866*fx5-0.333269719096257*fx6+0.49926922955688*fx7+0.509504608929686*fx8,theta,d+0.642615758240322548157075497020439535959501736363212695909875*delt.col(0),N2)%delt;
  fx10=f(x0+0.113976783964186*fx0-0.0768813364203357*fx5+0.239527360324391*fx6+0.397774662368095*fx7+0.0107558956873607*fx8-0.327769124164019*fx9,theta,d+0.357384241759677451842924502979560464040498263636787304090125*delt.col(0),N2)%delt;
  fx11=f(x0+0.0798314528280196*fx0-0.0520329686800603*fx5-0.0576954146168549*fx6+0.194781915712104*fx7+0.145384923188325*fx8-0.0782942710351671*fx9-0.114503299361099*fx10,theta,d+0.117472338035267653574498513020330924817132155731947880336209*delt.col(0),N2)%delt;
  fx12=f(x0+0.985115610164857*fx0+0.330885963040722*fx3+0.48966295730945*fx4-1.37896486574844*fx5-0.861164195027636*fx6+5.78428813637537*fx7+3.28807761985104*fx8-2.38633905093136*fx9-3.25479342483644*fx10-2.16343541686423*fx11,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt.col(0),N2)%delt;
  fx13=f(x0+0.895080295771633*fx0+0.197966831227192*fx2-0.0729547847313633*fx3-0.851236239662008*fx5+0.398320112318533*fx6+3.63937263181036*fx7+1.5482287703983*fx8-2.12221714704054*fx9-1.58350398545326*fx10-1.71561608285936*fx11-0.0244036405750127*fx12,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt.col(0),N2)%delt;
  fx14=f(x0+-0.915176561375291*fx0+1.45453440217827*fx1+0*fx2+0*fx3-0.777333643644968*fx4+0*fx5-0.0910895662155176*fx6+0.0910895662155176*fx12+0.777333643644968*fx13,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt.col(0),N2)%delt;
  fx15=f(x0+0.1*fx0-0.157178665799771*fx2+0.157178665799771*fx14,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt.col(0),N2)%delt;
  fx16=f(x0+0.181781300700095*fx0+0.675*fx1+0.34275815984719*fx2+0*fx3+0.259111214548323*fx4-0.358278966717952*fx5-1.04594895940883*fx6+0.930327845415627*fx7+1.77950959431708*fx8+0.1*fx9-0.282547569539044*fx10-0.159327350119973*fx11-0.145515894647002*fx12-0.259111214548323*fx13-0.34275815984719*fx14-0.675*fx15,theta,d+delt.col(0),N2)%delt;

  x0=x0+(0.0333333333333333333333333333333333333333333333333333333333333*fx0
        +0.0250000000000000000000000000000000000000000000000000000000000*fx1
        +0.0333333333333333333333333333333333333333333333333333333333333*fx2
        +0.000000000000000000000000000000000000000000000000000000000000*fx3
        +0.0500000000000000000000000000000000000000000000000000000000000*fx4
        +0.000000000000000000000000000000000000000000000000000000000000*fx5
        +0.0400000000000000000000000000000000000000000000000000000000000*fx6
        +0.000000000000000000000000000000000000000000000000000000000000*fx7
        +0.189237478148923490158306404106012326238162346948625830327194*fx8
        +0.277429188517743176508360262560654340428504319718040836339472*fx9
        +0.277429188517743176508360262560654340428504319718040836339472*fx10
        +0.189237478148923490158306404106012326238162346948625830327194*fx11
        -0.0400000000000000000000000000000000000000000000000000000000000*fx12
        -0.0500000000000000000000000000000000000000000000000000000000000*fx13
        -0.0333333333333333333333333333333333333333333333333333333333333*fx14
        -0.0250000000000000000000000000000000000000000000000000000000000*fx15
        +0.0333333333333333333333333333333333333333333333333333333333333*fx16);
    xe  =  abs(fx1.col(1)-fx15.col(1))/360.0;
    if(xe.max()>whch)
    {
      whch = xe.max();
    }
    d=d+delt.col(0);

  }
  '
}
}
if(RK.order==4)
{
 if(homo.res)
 {
 ODEpart= '
    return atemp;
}

// [[Rcpp::export]]
List  solver(vec Xs,vec Xt,vec theta,int N,double delt,int N2,vec tt,int P,double alpha,double lower,double upper,int tro)
{
    mat x0(N2,tro);
    mat xa(N2,tro);
    mat xe(N2,tro);
    mat fx1(N2,tro);
    mat fx2(N2,tro);
    mat fx3(N2,tro);
    mat fx4(N2,tro);
    mat fx5(N2,tro);
    mat fx6(N2,tro);
    double whch =0;
    x0.fill(0);
    x0.col(0)=Xs;
    vec d=tt;
    for (int i = 1; i < N+1; i++)
{

    fx1 = f(x0,theta,d,N2)*delt;
    fx2 = f(x0+0.25*fx1,theta,d+0.25*delt,N2)*delt;
    fx3 = f(x0+0.09375*fx1+0.28125*fx2,theta,d+0.375*delt,N2)*delt;
    fx4 = f(x0+0.879381*fx1-3.277196*fx2+ 3.320892*fx3,theta,d+0.9230769*delt,N2)*delt;
    fx5 = f(x0+2.032407*fx1-8*fx2+7.173489*fx3-0.2058967*fx4,theta,d+delt,N2)*delt;
    fx6 = f(x0-0.2962963*fx1+2*fx2-1.381676*fx3+0.4529727*fx4-0.275*fx5,theta,d+0.5*delt,N2)*delt;

    xa  =  x0+0.1185185*fx1+0.5189864*fx3+0.5061315*fx4-0.18*fx5+0.03636364*fx6;
    x0  =  x0+0.1157407*fx1+0.5489279*fx3+0.5353314*fx4-0.2*fx5;
    xe  =  abs(x0.col(1)-xa.col(1));
    if(xe.max()>whch)
    {
      whch = xe.max();
    }
    d=d+delt;
}
    '
 }
 if(!homo.res)
 {
 ODEpart= '
    return atemp;
}

// [[Rcpp::export]]
List  solver(vec Xs,vec Xt,vec theta,int N, mat delt,int N2,vec tt,int P,double alpha,double lower,double upper,int tro)
{
    mat x0(N2,tro);
    mat xa(N2,tro);
    mat xe(N2,tro);
    mat fx1(N2,tro);
    mat fx2(N2,tro);
    mat fx3(N2,tro);
    mat fx4(N2,tro);
    mat fx5(N2,tro);
    mat fx6(N2,tro);
    double whch =0;
    vec d=tt;
    x0.fill(0);
    x0.col(0)=Xs;
    for (int i = 1; i < N+1; i++)
{
    fx1 = f(x0,theta,d,N2)%delt;
    fx2 = f(x0+0.25*fx1,theta,d+0.25*delt.col(0),N2)%delt;
    fx3 = f(x0+0.09375*fx1+0.28125*fx2,theta,d+0.375*delt.col(0),N2)%delt;
    fx4 = f(x0+0.879381*fx1-3.277196*fx2+ 3.320892*fx3,theta,d+0.9230769*delt.col(0),N2)%delt;
    fx5 = f(x0+2.032407*fx1-8*fx2+7.173489*fx3-0.2058967*fx4,theta,d+delt.col(0),N2)%delt;
    fx6 = f(x0-0.2962963*fx1+2*fx2-1.381676*fx3+0.4529727*fx4-0.275*fx5,theta,d+0.5*delt.col(0),N2)%delt;
    xa  =  x0+0.1185185*fx1+0.5189864*fx3+0.5061315*fx4-0.18*fx5+0.03636364*fx6;
    x0  =  x0+0.1157407*fx1+0.5489279*fx3+0.5353314*fx4-0.2*fx5;
    xe  =  abs(x0.col(1)-xa.col(1));
    if(xe.max()>whch)
    {
      whch = xe.max();
    }
    d=d+delt.col(0);
}
    '
 }
}

if(DTR.order==2)
{
    Inv=
    '
  mat u(N2,2);
  u.col(0)=x0.col(0);
  u.col(1)=x0.col(1)+x0.col(0)%u.col(0);

  vec det=u.col(1)-u.col(0)%u.col(0);

  vec b11=(u.col(1))/det;
  vec b12=-u.col(0)/det;

  vec b21=-u.col(0)/det;
  vec b22=(1.0)/det;
  '
  switch(Dindex,
  {
    Dpart='
    vec val = -0.5*log(2*3.141592653589793)-0.5*log(x0.col(1))-0.5*((Xt-x0.col(0))%(Xt-x0.col(0))/x0.col(1));

    List ret;
    ret["like"] = val;
    ret["max"] = whch;
    return(ret);
}
    '
  })
}

if(DTR.order==4)
{
  Inv=
    '
  mat u(N2,4);
  u.col(0)=x0.col(0);
  u.col(1)=x0.col(1)+x0.col(0)%u.col(0);
  u.col(2)=x0.col(2)+x0.col(0)%u.col(1)+2*x0.col(1)%u.col(0);
  u.col(3)=x0.col(3)+x0.col(0)%u.col(2)+3*x0.col(1)%u.col(1)+3*x0.col(2)%u.col(0);
  
  vec det=u.col(1)%u.col(3)+u.col(0)%u.col(2)%u.col(1)+u.col(1)%u.col(0)%u.col(2)-u.col(2)%u.col(2)-u.col(1)%u.col(1)%u.col(1)-u.col(0)%u.col(0)%u.col(3);
  
  vec b11=(u.col(1)%u.col(3)-u.col(2)%u.col(2))/det;
  vec b12=(u.col(1)%u.col(2)-u.col(0)%u.col(3))/det;
  vec b13=(u.col(0)%u.col(2)-u.col(1)%u.col(1))/det;
  
  vec b21=(u.col(2)%u.col(1)-u.col(0)%u.col(3))/det;
  vec b22=(u.col(3)-u.col(1)%u.col(1))/det;
  vec b23=(u.col(1)%u.col(0)-u.col(2))/det;
  
  vec b31=(u.col(0)%u.col(2)-u.col(1)%u.col(1))/det;
  vec b32=(u.col(0)%u.col(1)-u.col(2))/det;
  vec b33=(u.col(1)-u.col(0)%u.col(0))/det;
  '
  switch(Dindex,
{
  Dpart='
  vec p=(1.0/3.0) *(3*(x0.col(3)/6.0)%x0.col(1) - pow(x0.col(2)/2.0,2))/pow(x0.col(3)/6.0,2);
  vec q=(1.0/27.0)*(27*pow(x0.col(3)/6.0,2)%(x0.col(0)-Xt) - 9*(x0.col(3)/6.0)%(x0.col(2)/2.0)%x0.col(1) + 2*pow(x0.col(2)/2.0,3))/pow(x0.col(3)/6.0,3);
  vec chk=pow(q,2)/4.0 + pow(p,3)/27.0;
  vec th=-(x0.col(2)/2.0)/(3*(x0.col(3)/6.0))+pow(-q/2.0+sqrt(chk),(1.0/3.0))-pow(q/2.0+sqrt(chk),(1.0/3.0));
  
  vec K =x0.col(0)%th+(x0.col(1)%th%th)/2.0+(x0.col(2)%th%th%th)/6.0 +(x0.col(3)%th%th%th%th)/24.0;
  vec K1=x0.col(0)   +(x0.col(1)%th)       +(x0.col(2)%th%th)/2.0    +(x0.col(3)%th%th%th)/6.0;
  vec K2=x0.col(1)   +(x0.col(2)%th)       +(x0.col(3)%th%th)/2.0;
  vec val=-0.5*log(2*3.141592653589793*K2)+(K-th%K1);

  List ret;
  ret["like"] = val;
  ret["max"] = whch;
  return(ret);
}
  '
},
{
  Dpart=
    '
  vec betas1 =+b12+2*b13%u.col(0);
  vec betas2 =+b22+2*b23%u.col(0);
  vec betas3 =+b32+2*b33%u.col(0);
  vec K=exp(-betas1*pow(Xtr,1)-0.5*betas2*pow(Xtr,2)-0.333333333333333333*betas3*pow(Xtr,3));
  
  vec lo =-(0.5/(u.col(0)-lower))%(sqrt(exp(2*alpha)+4*pow((u.col(0)-lower),2))-exp(alpha));
  vec up =-(0.5/(u.col(0)-upper))%(sqrt(exp(2*alpha)+4*pow((u.col(0)-upper),2))-exp(alpha));
  vec DT = (up-lo)/(1.00*P);
  vec tau   = exp(alpha)*lo/(1-pow(lo,2))+u.col(0);
  vec rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);
  vec K=exp(-betas1%pow(tau,1)-0.5*betas2%pow(tau,2)-0.333333333333333333*betas3%pow(tau,3))%rho%DT;
  for (int i = 1; i <= P; i++)
  {
    lo=lo+DT;
    tau   = exp(alpha)*lo/(1-pow(lo,2))+u.col(0);
    rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);
    K=K+exp(-betas1%pow(tau,1)-0.5*betas2%pow(tau,2)-0.333333333333333333*betas3%pow(tau,3))%rho%DT;
  }
  vec val = ((-betas1%pow(Xt,1)-0.5*betas2%pow(Xt,2)-0.333333333333333333*betas3%pow(Xt,3))-log(K));
  List ret;
  ret["like"] = val;
  ret["max"] = whch;
  return(ret);
}
  '
},
{
  Dpart='
  vec betas1 =+b11+2*b12%u.col(0)+3*b13%u.col(1);
  vec betas2 =+b21+2*b22%u.col(0)+3*b23%u.col(1);
  vec betas3 =+b31+2*b32%u.col(0)+3*b33%u.col(1);
  
  vec lo =-(0.5/(u.col(0)-lower))%(sqrt(exp(2*alpha)+4*pow((u.col(0)-lower),2))-exp(alpha));
  vec up =-(0.5/(u.col(0)-upper))%(sqrt(exp(2*alpha)+4*pow((u.col(0)-upper),2))-exp(alpha));
  vec DT = (up-lo)/(1.00*P);
  vec tau   = exp(alpha)*lo/(1-pow(lo,2))+u.col(0);
  vec rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);
  vec K=exp(-betas1%log(tau)+(-betas2%tau-0.5*betas3%pow(tau,2)))%rho%DT;
  for (int i = 1; i <= P; i++)
  {
    lo=lo+DT;
    tau   = exp(alpha)*lo/(1-pow(lo,2))+u.col(0);
    rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);
    K=K+exp(-betas1%log(tau)+(-betas2%tau-0.5*betas3%pow(tau,2)))%rho%DT;
  }
  vec val =((-betas1%log(Xt)+(-betas2%Xt-0.5*betas3%pow(Xt,2)))-log(K));
  List ret;
  ret["like"] = val;
  ret["max"] = whch;
  return(ret);
}
  '
  },
{
  Dpart='
  vec betas1 =+2*b11%u.col(0)+3*b12%u.col(1)+4*b13%u.col(2);
  vec betas2 =+2*b21%u.col(0)+3*b22%u.col(1)+4*b23%u.col(2);
  vec betas3 =+2*b31%u.col(0)+3*b32%u.col(1)+4*b33%u.col(2);
  
  vec lo =-(0.5/(u.col(0)-lower))%(sqrt(exp(2*alpha)+4*pow((u.col(0)-lower),2))-exp(alpha));
  vec up =-(0.5/(u.col(0)-upper))%(sqrt(exp(2*alpha)+4*pow((u.col(0)-upper),2))-exp(alpha));
  vec DT = (up-lo)/(1.00*P);
  vec tau   = exp(alpha)*lo/(1-pow(lo,2))+u.col(0);
  vec rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);
  vec K=exp(-betas2%log(tau)+(betas1/tau-betas3%tau))%rho%DT;
  for (int i = 1; i <= P; i++)
  {
    lo=lo+DT;
    tau   = exp(alpha)*lo/(1-pow(lo,2))+u.col(0);
    rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);
    K=K+exp(-betas2%log(tau)+(betas1/tau-betas3%tau))%rho%DT;
  }
  vec val((-betas2%log(Xt)+(betas1/Xt-betas3%Xt))-log(K));
  List ret;
  ret["like"] = val;
  ret["max"] = whch;
  return(ret);
}
  '
  },
{
  Dpart='
  
  vec betas1 =+b11%(1-2*u.col(0))+b12%(2*u.col(0)-3*u.col(1))+b13%(3*u.col(1)-4*u.col(2));
  vec betas2 =+b21%(1-2*u.col(0))+b22%(2*u.col(0)-3*u.col(1))+b23%(3*u.col(1)-4*u.col(2));
  vec betas3 =+b31%(1-2*u.col(0))+b32%(2*u.col(0)-3*u.col(1))+b33%(3*u.col(1)-4*u.col(2));
  
  vec lo =-(0.5/(u.col(0)-lower))%(sqrt(exp(2*alpha)+4*pow((u.col(0)-lower),2))-exp(alpha));
  vec up =-(0.5/(u.col(0)-upper))%(sqrt(exp(2*alpha)+4*pow((u.col(0)-upper),2))-exp(alpha));
  vec DT = (up-lo)/(1.00*P);
  vec tau   = exp(alpha)*lo/(1-pow(lo,2))+u.col(0);
  vec rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);
  vec K=exp(-betas2%log(tau)+(betas1/tau-betas3%tau))%rho%DT;
  for (int i = 1; i <= P; i++)
  {
    lo=lo+DT;
    tau   = exp(alpha)*lo/(1-pow(lo,2))+u.col(0);
    rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);
    K=K+exp(-betas2%log(tau)+(betas1/tau-betas3%tau))%rho%DT;
  }
  vec val = ((-betas2%log(Xt)+(betas1/Xt-betas3%Xt))-log(K));
  List ret;
  ret["like"] = val;
  ret["max"] = whch;
  return(ret);
}'
})
if(Dindex!=1)
{
  Dpart=paste(Inv,Dpart)
}
}
if(DTR.order==6)
{
  Inv=
    '
  vec u1 = x0.col(0);
  vec u2 = x0.col(1)+x0.col(0)%u1;
  vec u3 = x0.col(2)+x0.col(0)%u2+2*x0.col(1)%u1;
  vec u4 = x0.col(3)+x0.col(0)%u3+3*x0.col(1)%u2+3*x0.col(2)%u1 ;
  vec u5 = x0.col(4)+x0.col(0)%u4+4*x0.col(1)%u3+6*x0.col(2)%u2+4*x0.col(3)%u1 ;
  vec u6 = x0.col(5)+x0.col(0)%u5+5*x0.col(1)%u4+10*x0.col(2)%u3+10*x0.col(3)%u2+5*x0.col(4)%u1    ;
  
  vec det=-u1%(u1%(u4%u6-u5%u5)-u3%(u2%u6-u3%u5)+u4%(u2%u5-u3%u4))+u2%(u1%(u3%u6-u4%u5)-u2%(u2%u6-u3%u5)+u4%(u2%u4-u3%u3))+u2%(u4%u6-u5%u5)-u3%(u3%u6-u4%u5)-u3%(u1%(u3%u5-u4%u4)-u2%(u2%u5-u3%u4)+u3%(u2%u4-u3%u3))+u4%(u3%u5-u4%u4);
  
  vec b11 = (u2%(u4%u6-u5%u5)-u3%(u3%u6-u4%u5)+u4%(u3%u5-u4%u4))/det ;
  vec b12 = (-u1%(u4%u6-u5%u5)+u2%(u3%u6-u4%u5)-u3%(u3%u5-u4%u4))/det ;
  vec b13 = (u1%(u3%u6-u4%u5)-u2%(u2%u6-u4%u4)+u3%(u2%u5-u3%u4))/det  ;
  vec b14 = (-u1%(u3%u5-u4%u4)+u2%(u2%u5-u3%u4)-u3%(u2%u4-u3%u3))/det ;
  
  vec b21 = (-u1%(u4%u6-u5%u5)+u3%(u2%u6-u3%u5)-u4%(u2%u5-u3%u4))/det;
  vec b22 = (-u2%(u2%u6-u3%u5)+u4%u6-u5%u5+u3%(u2%u5-u3%u4))/det ;
  vec b23 = (u2%(u1%u6-u3%u4)-u3%u6-u3%(u1%u5-u3%u3)+u4%u5)/det  ;
  vec b24 = (-u2%(u1%u5-u2%u4)+u3%u5-u4%u4+u3%(u1%u4-u2%u3))/det ;
  
  vec b31 = (u1%(u3%u6-u4%u5)-u2%(u2%u6-u3%u5)+u4%(u2%u4-u3%u3))/det;
  vec b32 = (u1%(u2%u6-u3%u5)-u3%u6+u4%u5-u3%(u2%u4-u3%u3))/det ;
  vec b33 = (-u1%(u1%u6-u3%u4)+u2%u6-u4%u4+u3%(u1%u4-u2%u3))/det;
  vec b34 = (u1%(u1%u5-u2%u4)-u2%u5+u3%u4-u3%(u1%u3-u2%u2))/det ;
  
  vec b41 = (-u1%(u3%u5-u4%u4)+u2%(u2%u5-u3%u4)-u3%(u2%u4-u3%u3))/det;
  vec b42 = (-u1%(u2%u5-u3%u4)+u3%u5-u4%u4+u2%(u2%u4-u3%u3))/det;
  vec b43 = (u1%(u1%u5-u3%u3)-u2%u5-u2%(u1%u4-u2%u3)+u3%u4)/det ;
  vec b44 = (-u1%(u1%u4-u2%u3)+u2%u4-u3%u3+u2%(u1%u3-u2%u2))/det ;
  '
  switch(Dindex,
{
  Dpart='
  vec p=(1.0/3.0) *(3*(x0.col(3)/6.0)%x0.col(1) - pow(x0.col(2)/2.0,2))/pow(x0.col(3)/6.0,2);
  vec q=(1.0/27.0)*(27*pow(x0.col(3)/6.0,2)%(x0.col(0)-Xt) - 9*(x0.col(3)/6.0)%(x0.col(2)/2.0)%x0.col(1) + 2*pow(x0.col(2)/2.0,3))/pow(x0.col(3)/6.0,3);
  vec chk=pow(q,2)/4.0 + pow(p,3)/27.0;
  vec th=-(x0.col(2)/2.0)/(3*(x0.col(3)/6.0))+pow(-q/2.0+sqrt(chk),(1.0/3.0))-pow(q/2.0+sqrt(chk),(1.0/3.0));
  
  vec K =x0.col(0)%th+(x0.col(1)%th%th)/2.0+(x0.col(2)%th%th%th)/6.0 +(x0.col(3)%th%th%th%th)/24.0;
  vec K1=x0.col(0)   +(x0.col(1)%th)       +(x0.col(2)%th%th)/2.0    +(x0.col(3)%th%th%th)/6.0;
  vec K2=x0.col(1)   +(x0.col(2)%th)       +(x0.col(3)%th%th)/2.0;
  vec val=-0.5*log(2*3.141592653589793*K2)+(K-th%K1);
  List ret;
  ret["like"] = val;
  ret["max"] = whch;
  return(ret);
}
  '
},
{
  Dpart=
    '
  vec betas1 =+b12+2*b13%u1+3*b14%u2;
  vec betas2 =+b22+2*b23%u1+3*b24%u2;
  vec betas3 =+b32+2*b33%u1+3*b34%u2;
  vec betas4 =+b42+2*b43%u1+3*b44%u2;
  
  vec lo =-(0.5/(u1-lower))%(sqrt(exp(2*alpha)+4*pow((u1-lower),2))-exp(alpha));
  vec up =-(0.5/(u1-upper))%(sqrt(exp(2*alpha)+4*pow((u1-upper),2))-exp(alpha));
  vec DT = (up-lo)/(1.00*P);
  vec tau   = exp(alpha)*lo/(1-pow(lo,2))+u1;
  vec rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);
  vec K=exp(-betas1%pow(tau,1)-0.5*betas2%pow(tau,2)-0.333333333333333333*betas3%pow(tau,3)-0.25*betas4%pow(tau,4))%rho%DT;
  for (int i = 1; i <= P; i++)
  {
    lo=lo+DT;
    tau   = exp(alpha)*lo/(1-pow(lo,2))+u1;
    rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);
    K=K+exp(-betas1%pow(tau,1)-0.5*betas2%pow(tau,2)-0.333333333333333333*betas3%pow(tau,3)-0.25*betas4%pow(tau,4))%rho%DT;
  }
  vec val=((-betas1%pow(Xt,1)-0.5*betas2%pow(Xt,2)-0.333333333333333333*betas3%pow(Xt,3)-0.25*betas4%pow(Xt,4))-log(K));
  List ret;
  ret["like"] = val;
  ret["max"] = whch;
  return(ret);
}
  '
},
{
Dpart='
  
  vec betas1 =+b11+2*b12%u1+3*b13%u2+4*b14%u3;
  vec betas2 =+b21+2*b22%u1+3*b23%u2+4*b24%u3;
  vec betas3 =+b31+2*b32%u1+3*b33%u2+4*b34%u3;
  vec betas4 =+b41+2*b42%u1+3*b43%u2+4*b44%u3;

  vec lo =-(0.5/(u1-lower))%(sqrt(exp(2*alpha)+4*pow((u1-lower),2))-exp(alpha));
  vec up =-(0.5/(u1-upper))%(sqrt(exp(2*alpha)+4*pow((u1-upper),2))-exp(alpha));
  vec DT = (up-lo)/(1.00*P);
  vec tau   = exp(alpha)*lo/(1-pow(lo,2))+u1;
  vec rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);
  vec K=exp(-betas1%log(tau)+(-betas2%tau-0.5*betas3%pow(tau,2)-0.33333333333333333333333*betas4%pow(tau,3)))%rho%DT;
  for (int i = 1; i <= P; i++)
  {
    lo=lo+DT;
    tau   = exp(alpha)*lo/(1-pow(lo,2))+u1;
    rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);
    K=K+exp(-betas1%log(tau)+(-betas2%tau-0.5*betas3%pow(tau,2)-0.33333333333333333333333*betas4%pow(tau,3)))%rho%DT;
  }
  vec val = ((-betas1%log(Xt)+(-betas2%Xt-0.5*betas3%pow(Xt,2)-0.33333333333333333333333*betas4%pow(Xt,3)))-log(K));
  List ret;
  ret["like"] = val;
  ret["max"] = whch;
  return(ret);
}
'
},
{
Dpart='

vec betas1 =+2*b11%u1+3*b12%u2+4*b13%u3+5*b14%u4;
vec betas2 =+2*b21%u1+3*b22%u2+4*b23%u3+5*b24%u4;
vec betas3 =+2*b31%u1+3*b32%u2+4*b33%u3+5*b34%u4;
vec betas4 =+2*b41%u1+3*b42%u2+4*b43%u3+5*b44%u4;

  vec lo =-(0.5/(u1-lower))%(sqrt(exp(2*alpha)+4*pow((u1-lower),2))-exp(alpha));
  vec up =-(0.5/(u1-upper))%(sqrt(exp(2*alpha)+4*pow((u1-upper),2))-exp(alpha));
  vec DT = (up-lo)/(1.00*P);
  vec tau   = exp(alpha)*lo/(1-pow(lo,2))+u1;
  vec rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);

  vec K=exp(-betas2%log(tau)+(betas1/tau-betas3%tau-0.5*betas4%pow(tau,2)))%rho%DT;
  for (int i = 1; i <= P; i++)
  {
    lo=lo+DT;
    tau   = exp(alpha)*lo/(1-pow(lo,2))+u1;
    rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);
    K=K+exp(-betas2%log(tau)+(betas1/tau-betas3%tau-0.5*betas4%pow(tau,2)))%rho%DT;
  }
  vec val = ((-betas2%log(Xt)+(betas1/Xt-betas3%Xt-0.5*betas4%pow(Xt,2)))-log(K));
  List ret;
  ret["like"] = val;
  ret["max"] = whch;
  return(ret);
}
'
},
{
 Dpart='
 vec betas1 =+b11%(1-2*u1)+b12%(2*u1-3*u2)+b13%(3*u2-4*u3)+b14%(4*u3-5*u4);
 vec betas2 =+b21%(1-2*u1)+b22%(2*u1-3*u2)+b23%(3*u2-4*u3)+b24%(4*u3-5*u4);
 vec betas3 =+b31%(1-2*u1)+b32%(2*u1-3*u2)+b33%(3*u2-4*u3)+b34%(4*u3-5*u4);
 vec betas4 =+b41%(1-2*u1)+b42%(2*u1-3*u2)+b43%(3*u2-4*u3)+b44%(4*u3-5*u4);
 
 vec K=exp(-betas1*log(Xtr)+(betas1+betas2+betas3+betas4)*log(1-Xtr)+(betas3+betas4)*Xtr+0.5*betas4*pow(Xtr,2));
 for (int i = 1; i <= lim; i++)
{
 Xtr=Xtr+delt2;
 K=K+exp(-betas1*log(Xtr)+(betas1+betas2+betas3+betas4)*log(1-Xtr)+(betas3+betas4)*Xtr+0.5*betas4*pow(Xtr,2))*delt2;
}
 vec val =((-betas1%log(Xt)+(betas1+betas2+betas3+betas4)%log(1-Xt)+(betas3+betas4)%Xt+0.5*betas4%pow(Xt,2))-log(K));
  List ret;
  ret["like"] = val;
  ret["max"] = whch;
  return(ret);
 }
 '
})
 if(Dindex!=1)
 {
    Dpart=paste(Inv,Dpart)
 }
}
if(DTR.order==8)
{
  Inv='
  vec u1=                                                                                  x0.col(0);
  vec u2=                                                                       x0.col(1)+1*x0.col(0)%u1;
  vec u3=                                                            x0.col(2)+1*x0.col(0)%u2+2*x0.col(1)%u1;
  vec u4=                                                 x0.col(3)+1*x0.col(0)%u3+3*x0.col(1)%u2+3*x0.col(2)%u1;
  vec u5=                                      x0.col(4)+1*x0.col(0)%u4+4*x0.col(1)%u3+6*x0.col(2)%u2+4*x0.col(3)%u1;
  vec u6=                         x0.col(5)+1*x0.col(0)%u5+5*x0.col(1)%u4+10*x0.col(2)%u3+10*x0.col(3)%u2+5*x0.col(4)%u1;
  vec u7=             x0.col(6)+1*x0.col(0)%u6+6*x0.col(1)%u5+15*x0.col(2)%u4+20*x0.col(3)%u3+15*x0.col(4)%u2+6*x0.col(5)%u1;
  vec u8= x0.col(7)+1*x0.col(0)%u7+7*x0.col(1)%u6+21*x0.col(2)%u5+35*x0.col(3)%u4+35*x0.col(4)%u3+21*x0.col(5)%u2+7*x0.col(6)%u1;
  
  vec det=-u1%(u1%(u4%(u6%u8-u7%u7)-u5%(u5%u8-u6%u7)+u6%(u5%u7-u6%u6))-u3%(u2%(u6%u8-u7%u7)-u5%(u3%u8-u4%u7)+u6%(u3%u7-u4%u6))+u4%
  (u2%(u5%u8-u6%u7)-u4%(u3%u8-u4%u7)+u6%(u3%u6-u4%u5))-u5%(u2%(u5%u7-u6%u6)-u4%(u3%u7-u4%u6)+u5%(u3%u6-u4%u5)))+u2%(u1%
  (u3%(u6%u8-u7%u7)-u5%(u4%u8-u5%u7)+u6%(u4%u7-u5%u6))-u2%(u2%(u6%u8-u7%u7)-u5%(u3%u8-u4%u7)+u6%(u3%u7-u4%u6))+u4%
  (u2%(u4%u8-u5%u7)-u3%(u3%u8-u4%u7)+(u3%u5-u4%u4)%u6)-u5%(u2%(u4%u7-u5%u6)-u3%(u3%u7-u4%u6)+u5%(u3%u5-u4%u4)))-u3%(u1%
  (u3%(u5%u8-u6%u7)-u4%(u4%u8-u5%u7)+u6%(u4%u6-u5%u5))-u2%(u2%(u5%u8-u6%u7)-u4%(u3%u8-u4%u7)+u6%(u3%u6-u4%u5))+u3%
  (u2%(u4%u8-u5%u7)-u3%(u3%u8-u4%u7)+(u3%u5-u4%u4)%u6)-u5%(u2%(u4%u6-u5%u5)-u3%(u3%u6-u4%u5)+u4%(u3%u5-u4%u4)))+u2%
  (u4%(u6%u8-u7%u7)-u5%(u5%u8-u6%u7)+u6%(u5%u7-u6%u6))-u3%(u3%(u6%u8-u7%u7)-u5%(u4%u8-u5%u7)+u6%(u4%u7-u5%u6))+u4%
  (u3%(u5%u8-u6%u7)-u4%(u4%u8-u5%u7)+u6%(u4%u6-u5%u5))+u4%(u1%(u3%(u5%u7-u6%u6)-u4%(u4%u7-u5%u6)+u5%(u4%u6-u5%u5))-u2%
  (u2%(u5%u7-u6%u6)-u4%(u3%u7-u4%u6)+u5%(u3%u6-u4%u5))+u3%(u2%(u4%u7-u5%u6)-u3%(u3%u7-u4%u6)+u5%(u3%u5-u4%u4))-u4%
  (u2%(u4%u6-u5%u5)-u3%(u3%u6-u4%u5)+u4%(u3%u5-u4%u4)))-u5%(u3%(u5%u7-u6%u6)-u4%(u4%u7-u5%u6)+u5%(u4%u6-u5%u5));
  
  vec b11= (u2%(u4%(u6%u8-u7%u7)-u5%(u5%u8-u6%u7)+u6%(u5%u7-u6%u6))-u3%(u3%(u6%u8-u7%u7)-u5%(u4%u8-u5%u7)+u6%(u4%u7-u5%u6))+u4%(u3%(u5%u8-u6%u7)-u4%(u4%u8-u5%u7)+u6%(u4%u6-u5%u5))-u5%(u3%(u5%u7-u6%u6)-u4%(u4%u7-u5%u6)+u5%(u4%u6-u5%u5)))/det;
  vec b12= (-u1%(u4%(u6%u8-u7%u7)-u5%(u5%u8-u6%u7)+u6%(u5%u7-u6%u6))+u2%(u3%(u6%u8-u7%u7)-u5%(u4%u8-u5%u7)+u6%(u4%u7-u5%u6))-u3%(u3%(u5%u8-u6%u7)-u4%(u4%u8-u5%u7)+u6%(u4%u6-u5%u5))+u4%(u3%(u5%u7-u6%u6)-u4%(u4%u7-u5%u6)+u5%(u4%u6-u5%u5)))/det;
  vec b13= (u1%(u3%(u6%u8-u7%u7)-u4%(u5%u8-u6%u7)+u5%(u5%u7-u6%u6))-u2%(u2%(u6%u8-u7%u7)-u4%(u4%u8-u5%u7)+u5%(u4%u7-u5%u6))+u3%(u2%(u5%u8-u6%u7)-u3%(u4%u8-u5%u7)+u5%(u4%u6-u5%u5))-u4%(u2%(u5%u7-u6%u6)-u3%(u4%u7-u5%u6)+u4%(u4%u6-u5%u5)))/det;
  vec b14= (-u1%(u3%(u5%u8-u6%u7)-u4%(u4%u8-u6%u6)+u5%(u4%u7-u5%u6))+u2%(u2%(u5%u8-u6%u7)-u4%(u3%u8-u5%u6)+u5%(u3%u7-u5%u5))-u3%(u2%(u4%u8-u6%u6)-u3%(u3%u8-u5%u6)+u5%(u3%u6-u4%u5))+u4%(u2%(u4%u7-u5%u6)-u3%(u3%u7-u5%u5)+u4%(u3%u6-u4%u5)))/det;
  vec b15= (u1%(u3%(u5%u7-u6%u6)-u4%(u4%u7-u5%u6)+u5%(u4%u6-u5%u5))-u2%(u2%(u5%u7-u6%u6)-u4%(u3%u7-u4%u6)+u5%(u3%u6-u4%u5))+u3%(u2%(u4%u7-u5%u6)-u3%(u3%u7-u4%u6)+u5%(u3%u5-u4%u4))-u4%(u2%(u4%u6-u5%u5)-u3%(u3%u6-u4%u5)+u4%(u3%u5-u4%u4)))/det;
  
  vec b21= (-u1%(u4%(u6%u8-u7%u7)-u5%(u5%u8-u6%u7)+u6%(u5%u7-u6%u6))+u3%(u2%(u6%u8-u7%u7)-u5%(u3%u8-u4%u7)+u6%(u3%u7-u4%u6))-u4%(u2%(u5%u8-u6%u7)-u4%(u3%u8-u4%u7)+u6%(u3%u6-u4%u5))+u5%(u2%(u5%u7-u6%u6)-u4%(u3%u7-u4%u6)+u5%(u3%u6-u4%u5)))/det;
  vec b22= (-u2%(u2%(u6%u8-u7%u7)-u5%(u3%u8-u4%u7)+u6%(u3%u7-u4%u6))+u3%(u2%(u5%u8-u6%u7)-u4%(u3%u8-u4%u7)+u6%(u3%u6-u4%u5))+u4%(u6%u8-u7%u7)-u5%(u5%u8-u6%u7)-u4%(u2%(u5%u7-u6%u6)-u4%(u3%u7-u4%u6)+u5%(u3%u6-u4%u5))+u6%(u5%u7-u6%u6))/det;
  vec b23= (u2%(u1%(u6%u8-u7%u7)-u4%(u3%u8-u4%u7)+u5%(u3%u7-u4%u6))-u3%(u1%(u5%u8-u6%u7)-u3%(u3%u8-u4%u7)+u5%(u3%u6-u4%u5))-u3%(u6%u8-u7%u7)+u4%(u5%u8-u6%u7)+u4%(u1%(u5%u7-u6%u6)-u3%(u3%u7-u4%u6)+u4%(u3%u6-u4%u5))-u5%(u5%u7-u6%u6))/det;
  vec b24= (-u2%(u1%(u5%u8-u6%u7)-u4%(u2%u8-u4%u6)+u5%(u2%u7-u4%u5))+u3%(u1%(u4%u8-u6%u6)-u3%(u2%u8-u4%u6)+u5%(u2%u6-u4%u4))+u3%(u5%u8-u6%u7)-u4%(u4%u8-u6%u6)-u4%(u1%(u4%u7-u5%u6)-u3%(u2%u7-u4%u5)+u4%(u2%u6-u4%u4))+u5%(u4%u7-u5%u6))/det;
  vec b25= (u2%(u1%(u5%u7-u6%u6)-u4%(u2%u7-u3%u6)+u5%(u2%u6-u3%u5))-u3%(u1%(u4%u7-u5%u6)-u3%(u2%u7-u3%u6)+u5%(u2%u5-u3%u4))-u3%(u5%u7-u6%u6)+u4%(u4%u7-u5%u6)+u4%(u1%(u4%u6-u5%u5)-u3%(u2%u6-u3%u5)+u4%(u2%u5-u3%u4))-u5%(u4%u6-u5%u5))/det;
  
  vec b31= (u1%(u3%(u6%u8-u7%u7)-u5%(u4%u8-u5%u7)+u6%(u4%u7-u5%u6))-u2%(u2%(u6%u8-u7%u7)-u5%(u3%u8-u4%u7)+u6%(u3%u7-u4%u6))+u4%(u2%(u4%u8-u5%u7)-u3%(u3%u8-u4%u7)+(u3%u5-u4%u4)%u6)-u5%(u2%(u4%u7-u5%u6)-u3%(u3%u7-u4%u6)+u5%(u3%u5-u4%u4)))/det;
  vec b32= (u1%(u2%(u6%u8-u7%u7)-u5%(u3%u8-u4%u7)+u6%(u3%u7-u4%u6))-u3%(u2%(u4%u8-u5%u7)-u3%(u3%u8-u4%u7)+(u3%u5-u4%u4)%u6)-u3%(u6%u8-u7%u7)+u5%(u4%u8-u5%u7)+u4%(u2%(u4%u7-u5%u6)-u3%(u3%u7-u4%u6)+u5%(u3%u5-u4%u4))-u6%(u4%u7-u5%u6))/det;
  vec b33= (-u1%(u1%(u6%u8-u7%u7)-u4%(u3%u8-u4%u7)+u5%(u3%u7-u4%u6))+u3%(u1%(u4%u8-u5%u7)-u2%(u3%u8-u4%u7)+u5%(u3%u5-u4%u4))+u2%(u6%u8-u7%u7)-u4%(u4%u8-u5%u7)-u4%(u1%(u4%u7-u5%u6)-u2%(u3%u7-u4%u6)+u4%(u3%u5-u4%u4))+u5%(u4%u7-u5%u6))/det;
  vec b34= (u1%(u1%(u5%u8-u6%u7)-u4%(u2%u8-u4%u6)+u5%(u2%u7-u4%u5))-u3%(u1%(u3%u8-u5%u6)-u2%(u2%u8-u4%u6)+u5%(u2%u5-u3%u4))-u2%(u5%u8-u6%u7)+u4%(u3%u8-u5%u6)+u4%(u1%(u3%u7-u5%u5)-u2%(u2%u7-u4%u5)+u4%(u2%u5-u3%u4))-u5%(u3%u7-u5%u5))/det;
  vec b35= (-u1%(u1%(u5%u7-u6%u6)-u4%(u2%u7-u3%u6)+u5%(u2%u6-u3%u5))+u3%(u1%(u3%u7-u4%u6)-u2%(u2%u7-u3%u6)+(u2%u4-u3%u3)%u5)+u2%(u5%u7-u6%u6)-u4%(u3%u7-u4%u6)-u4%(u1%(u3%u6-u4%u5)-u2%(u2%u6-u3%u5)+u4%(u2%u4-u3%u3))+u5%(u3%u6-u4%u5))/det;
  
  vec b41= (-u1%(u3%(u5%u8-u6%u7)-u4%(u4%u8-u5%u7)+u6%(u4%u6-u5%u5))+u2%(u2%(u5%u8-u6%u7)-u4%(u3%u8-u4%u7)+u6%(u3%u6-u4%u5))-u3%(u2%(u4%u8-u5%u7)-u3%(u3%u8-u4%u7)+(u3%u5-u4%u4)%u6)+u5%(u2%(u4%u6-u5%u5)-u3%(u3%u6-u4%u5)+u4%(u3%u5-u4%u4)))/det;
  vec b42= (-u1%(u2%(u5%u8-u6%u7)-u4%(u3%u8-u4%u7)+u6%(u3%u6-u4%u5))+u2%(u2%(u4%u8-u5%u7)-u3%(u3%u8-u4%u7)+(u3%u5-u4%u4)%u6)+u3%(u5%u8-u6%u7)-u4%(u4%u8-u5%u7)-u4%(u2%(u4%u6-u5%u5)-u3%(u3%u6-u4%u5)+u4%(u3%u5-u4%u4))+u6%(u4%u6-u5%u5))/det;
  vec b43= (u1%(u1%(u5%u8-u6%u7)-u3%(u3%u8-u4%u7)+u5%(u3%u6-u4%u5))-u2%(u1%(u4%u8-u5%u7)-u2%(u3%u8-u4%u7)+u5%(u3%u5-u4%u4))-u2%(u5%u8-u6%u7)+u3%(u4%u8-u5%u7)+u4%(u1%(u4%u6-u5%u5)-u2%(u3%u6-u4%u5)+u3%(u3%u5-u4%u4))-u5%(u4%u6-u5%u5))/det;
  vec b44= (-u1%(u1%(u4%u8-u6%u6)-u3%(u2%u8-u4%u6)+u5%(u2%u6-u4%u4))+u2%(u1%(u3%u8-u5%u6)-u2%(u2%u8-u4%u6)+u5%(u2%u5-u3%u4))+u2%(u4%u8-u6%u6)-u3%(u3%u8-u5%u6)-u4%(u1%(u3%u6-u4%u5)-u2%(u2%u6-u4%u4)+u3%(u2%u5-u3%u4))+u5%(u3%u6-u4%u5))/det;
  vec b45= (u1%(u1%(u4%u7-u5%u6)-u3%(u2%u7-u3%u6)+u5%(u2%u5-u3%u4))-u2%(u1%(u3%u7-u4%u6)-u2%(u2%u7-u3%u6)+(u2%u4-u3%u3)%u5)-u2%(u4%u7-u5%u6)+u3%(u3%u7-u4%u6)+u4%(u1%(u3%u5-u4%u4)-u2%(u2%u5-u3%u4)+u3%(u2%u4-u3%u3))-u5%(u3%u5-u4%u4))/det;
  
  vec b51= (u1%(u3%(u5%u7-u6%u6)-u4%(u4%u7-u5%u6)+u5%(u4%u6-u5%u5))-u2%(u2%(u5%u7-u6%u6)-u4%(u3%u7-u4%u6)+u5%(u3%u6-u4%u5))+u3%(u2%(u4%u7-u5%u6)-u3%(u3%u7-u4%u6)+u5%(u3%u5-u4%u4))-u4%(u2%(u4%u6-u5%u5)-u3%(u3%u6-u4%u5)+u4%(u3%u5-u4%u4)))/det;
  vec b52= (u1%(u2%(u5%u7-u6%u6)-u4%(u3%u7-u4%u6)+u5%(u3%u6-u4%u5))-u2%(u2%(u4%u7-u5%u6)-u3%(u3%u7-u4%u6)+u5%(u3%u5-u4%u4))-u3%(u5%u7-u6%u6)+u4%(u4%u7-u5%u6)+u3%(u2%(u4%u6-u5%u5)-u3%(u3%u6-u4%u5)+u4%(u3%u5-u4%u4))-u5%(u4%u6-u5%u5))/det;
  vec b53= (-u1%(u1%(u5%u7-u6%u6)-u3%(u3%u7-u4%u6)+u4%(u3%u6-u4%u5))+u2%(u1%(u4%u7-u5%u6)-u2%(u3%u7-u4%u6)+u4%(u3%u5-u4%u4))+u2%(u5%u7-u6%u6)-u3%(u4%u7-u5%u6)-u3%(u1%(u4%u6-u5%u5)-u2%(u3%u6-u4%u5)+u3%(u3%u5-u4%u4))+u4%(u4%u6-u5%u5))/det;
  vec b54= (u1%(u1%(u4%u7-u5%u6)-u3%(u2%u7-u4%u5)+u4%(u2%u6-u4%u4))-u2%(u1%(u3%u7-u5%u5)-u2%(u2%u7-u4%u5)+u4%(u2%u5-u3%u4))-u2%(u4%u7-u5%u6)+u3%(u3%u7-u5%u5)+u3%(u1%(u3%u6-u4%u5)-u2%(u2%u6-u4%u4)+u3%(u2%u5-u3%u4))-u4%(u3%u6-u4%u5))/det;
  vec b55= (-u1%(u1%(u4%u6-u5%u5)-u3%(u2%u6-u3%u5)+u4%(u2%u5-u3%u4))+u2%(u1%(u3%u6-u4%u5)-u2%(u2%u6-u3%u5)+u4%(u2%u4-u3%u3))+u2%(u4%u6-u5%u5)-u3%(u3%u6-u4%u5)-u3%(u1%(u3%u5-u4%u4)-u2%(u2%u5-u3%u4)+u3%(u2%u4-u3%u3))+u4%(u3%u5-u4%u4))/det;
  '
  
  
  switch(Dindex,
{
  Dpart='
  vec p=(1.0/3.0) *(3*(x0.col(3)/6.0)%x0.col(1) - pow(x0.col(2)/2.0,2))/pow(x0.col(3)/6.0,2);
  vec q=(1.0/27.0)*(27*pow(x0.col(3)/6.0,2)%(x0.col(0)-Xt) - 9*(x0.col(3)/6.0)%(x0.col(2)/2.0)%x0.col(1) + 2*pow(x0.col(2)/2.0,3))/pow(x0.col(3)/6.0,3);
  vec chk=pow(q,2)/4.0 + pow(p,3)/27.0;
  vec th=-(x0.col(2)/2.0)/(3*(x0.col(3)/6.0))+pow(-q/2.0+sqrt(chk),(1.0/3.0))-pow(q/2.0+sqrt(chk),(1.0/3.0));
  
  vec K =x0.col(0)%th+(x0.col(1)%th%th)/2.0+(x0.col(2)%th%th%th)/6.0 +(x0.col(3)%th%th%th%th)/24.0;
  vec K1=x0.col(0)   +(x0.col(1)%th)       +(x0.col(2)%th%th)/2.0    +(x0.col(3)%th%th%th)/6.0;
  vec K2=x0.col(1)   +(x0.col(2)%th)       +(x0.col(3)%th%th)/2.0;
  vec val=-0.5*log(2*3.141592653589793*K2)+(K-th%K1);
  List ret;
  ret["like"] = val;
  ret["max"] = whch;
  return(ret);
}
  '
},
{
  Dpart=
    '
  vec betas1 =+b12+2*b13%u1+3*b14%u2+4*b15%u3;
  vec betas2 =+b22+2*b23%u1+3*b24%u2+4*b25%u3;
  vec betas3 =+b32+2*b33%u1+3*b34%u2+4*b35%u3;
  vec betas4 =+b42+2*b43%u1+3*b44%u2+4*b45%u3;
  vec betas5 =+b52+2*b53%u1+3*b54%u2+4*b55%u3;
  
  vec lo =-(0.5/(u1-lower))%(sqrt(exp(2*alpha)+4*pow((u1-lower),2))-exp(alpha));
  vec up =-(0.5/(u1-upper))%(sqrt(exp(2*alpha)+4*pow((u1-upper),2))-exp(alpha));
  vec DT = (up-lo)/(1.00*P);
  vec tau   = exp(alpha)*lo/(1-pow(lo,2))+u1;
  vec rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);
  vec K=exp(-betas1%pow(tau,1)-0.5*betas2%pow(tau,2)-0.333333333333333333*betas3%pow(tau,3)-0.25*betas4%pow(tau,4)-0.2*betas5%pow(tau,5))%rho%DT;
  for (int i = 1; i <= P; i++)
  {
    lo=lo+DT;
    tau   = exp(alpha)*lo/(1-pow(lo,2))+u1;
    rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);
    K=K+exp(-betas1%pow(tau,1)-0.5*betas2%pow(tau,2)-0.333333333333333333*betas3%pow(tau,3)-0.25*betas4%pow(tau,4)-0.2*betas5%pow(tau,5))%rho%DT;
  }

  vec val =((-betas1%pow(Xt,1)-0.5*betas2%pow(Xt,2)-0.333333333333333333*betas3%pow(Xt,3)-0.25*betas4%pow(Xt,4)-0.2*betas5%pow(Xt,5))-log(K));
  List ret;
  ret["like"] = val;
  ret["max"] = whch;
  return(ret);
}
  '
},
{
  Dpart='
  vec betas1 =+b11+2*b12%u1+3*b13%u2+4*b14%u3+4*b15%u4;
  vec betas2 =+b21+2*b22%u1+3*b23%u2+4*b24%u3+4*b25%u4;
  vec betas3 =+b31+2*b32%u1+3*b33%u2+4*b34%u3+4*b35%u4;
  vec betas4 =+b41+2*b42%u1+3*b43%u2+4*b44%u3+4*b45%u4;
  vec betas5 =+b51+2*b52%u1+3*b53%u2+4*b54%u3+4*b55%u4;
  
  vec lo =-(0.5/(u1-lower))%(sqrt(exp(2*alpha)+4*pow((u1-lower),2))-exp(alpha));
  vec up =-(0.5/(u1-upper))%(sqrt(exp(2*alpha)+4*pow((u1-upper),2))-exp(alpha));
  vec DT = (up-lo)/(1.00*P);
  vec tau   = exp(alpha)*lo/(1-pow(lo,2))+u1;
  vec rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);
  
  vec K=exp(-betas1%log(tau)+(-betas2%tau-0.5*betas3%pow(tau,2)-0.33333333333333333333333*betas4%pow(tau,3)-0.25*betas5%pow(tau,4)))%rho%DT;
  for (int i = 1; i <= P; i++)
  {
    lo=lo+DT;
    tau   = exp(alpha)*lo/(1-pow(lo,2))+u1;
    rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);
    K=K+exp(-betas1%log(tau)+(-betas2%tau-0.5*betas3%pow(tau,2)-0.33333333333333333333333*betas4%pow(tau,3)-0.25*betas5%pow(tau,4)))%rho%DT;
  }
  vec val =((-betas1%log(Xt)+(-betas2%Xt-0.5*betas3%pow(Xt,2)-0.33333333333333333333333*betas4%pow(Xt,3)-0.25*betas5%pow(Xt,4)))-log(K));
  List ret;
  ret["like"] = val;
  ret["max"] = whch;
  return(ret);
}
  '
},
{
  Dpart='
  vec betas1 =+2*b11%u1+3*b12%u2+4*b13%u3+5*b14%u4+6*b15%u5;
  vec betas2 =+2*b21%u1+3*b22%u2+4*b23%u3+5*b24%u4+6*b25%u5;
  vec betas3 =+2*b31%u1+3*b32%u2+4*b33%u3+5*b34%u4+6*b35%u5;
  vec betas4 =+2*b41%u1+3*b42%u2+4*b43%u3+5*b44%u4+6*b45%u5;
  vec betas5 =+2*b51%u1+3*b52%u2+4*b53%u3+5*b54%u4+6*b55%u5;
  
  vec lo =-(0.5/(u1-lower))%(sqrt(exp(2*alpha)+4*pow((u1-lower),2))-exp(alpha));
  vec up =-(0.5/(u1-upper))%(sqrt(exp(2*alpha)+4*pow((u1-upper),2))-exp(alpha));
  vec DT = (up-lo)/(1.00*P);
  vec tau   = exp(alpha)*lo/(1-pow(lo,2))+u1;
  vec rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);

  vec K=exp(-betas2%log(tau)+(betas1/tau-betas3%tau-0.5*betas4%pow(tau,2)-0.3333333333333333*betas5%pow(tau,3)))%rho%DT;
  for (int i = 1; i <= P; i++)
  {
    lo=lo+DT;
    tau   = exp(alpha)*lo/(1-pow(lo,2))+u1;
    rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);
    K=K+exp(-betas2%log(tau)+(betas1/tau-betas3%tau-0.5*betas4%pow(tau,2)-0.3333333333333333*betas5%pow(tau,3)))%rho%DT;
  }
  vec val =((-betas2%log(Xt)+(betas1/Xt-betas3%Xt-0.5*betas4%pow(Xt,2)-0.3333333333333333*betas5%pow(Xt,3)))-log(K));
  List ret;
  ret["like"] = val;
  ret["max"] = whch;
  return(ret);
}
  '
  },
{
  Dpart='
  
  vec betas1 =+b11%(1-2*u1)+b12%(2*u1-3*u2)+b13%(3*u2-4*u3)+b14%(4*u3-5*u4)+b15%(5*u4-6*u5);
  vec betas2 =+b21%(1-2*u1)+b22%(2*u1-3*u2)+b23%(3*u2-4*u3)+b24%(4*u3-5*u4)+b25%(5*u4-6*u5);
  vec betas3 =+b31%(1-2*u1)+b32%(2*u1-3*u2)+b33%(3*u2-4*u3)+b34%(4*u3-5*u4)+b35%(5*u4-6*u5);
  vec betas4 =+b41%(1-2*u1)+b42%(2*u1-3*u2)+b43%(3*u2-4*u3)+b44%(4*u3-5*u4)+b45%(5*u4-6*u5);
  vec betas5 =+b51%(1-2*u1)+b52%(2*u1-3*u2)+b53%(3*u2-4*u3)+b54%(4*u3-5*u4)+b55%(5*u4-6*u5);
  
  vec lo =-(0.5/(u1-lower))%(sqrt(exp(2*alpha)+4*pow((u1-lower),2))-exp(alpha));
  vec up =-(0.5/(u1-upper))%(sqrt(exp(2*alpha)+4*pow((u1-upper),2))-exp(alpha));
  vec DT = (up-lo)/(1.00*P);
  vec tau   = exp(alpha)*lo/(1-pow(lo,2))+u1;
  vec rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);

  vec K=exp(-betas2%log(tau)+(betas1/tau-betas3%tau-0.5*betas4%pow(tau,2)-0.3333333333333333*betas5%pow(tau,3)))%rho%DT;
  for (int i = 1; i <= P; i++)
  {
    lo=lo+DT;
    tau   = exp(alpha)*lo/(1-pow(lo,2))+u1;
    rho   = exp(alpha)*(1+pow(lo,2))/pow((1-pow(lo,2)),2);
    K=K+exp(-betas2%log(tau)+(betas1/tau-betas3%tau-0.5*betas4%pow(tau,2)-0.3333333333333333*betas5%pow(tau,3)))%rho%DT;
  }
  vec val = ((-betas2%log(Xt)+(betas1/Xt-betas3%Xt-0.5*betas4%pow(Xt,2)-0.3333333333333333*betas5%pow(Xt,3)))-log(K));
  List ret;
  ret["like"] = val;
  ret["max"] = whch;
  return(ret);
}
  '
  })
if(Dindex!=1)
{
  Dpart=paste(Inv,Dpart)
}
}






   #==============================================================================
   #                                Syntax Check
   #==============================================================================

    strcheck=function(str)
    {
    k=0
    if(length(grep('/',str))!=0)
    {
     warning('C++ uses integer division when denominators are of type int. Use decimal values if possible (e.g. 0.5 i.s.o. 1/2).',call. = F)
     k=k+1
    }
    if(length(grep('\\^',str))!=0)
    {
     stop('^-Operator not defined in C++: Use pow(x,p) instead (e.g. pow(x,2) i.s.o. x^2).',call. = F)
     k=k+1
    }
    return(k)
   }
   for(i in which(func.list==1))
   {
     strcheck(body(namess[i])[2])
   }

   #==============================================================================
   #                   Generate TYPE of Solution
   #==============================================================================
   if(state1)
   {
   # DATA RESOLUTION -------------------------------------------------------------
   if(!homo.res)
   {
    delt=cbind(diff(T.seq)/mesh,diff(T.seq)/mesh)
    if(!is.null(exclude))
    {
      diffs=diff(T.seq)/mesh
      diffs[excl]=1/20/mesh
      delt=cbind(diffs,diffs)
    }
   }
   # REFERENCE MATRIX ------------------------------------------------------------
   MAT=rbind(
   c('(1+0*a.col(0))','a.col(0)' ,''),
   c('','2*a.col(1)','(1+0*a.col(0))'),
   c('','3*a.col(2)',''),
   c('','4*a.col(3)',''))

   # HOMOGENEITY -----------------------------------------------------------------
   namess2=c('G0','G1','Q0')
   func.list2=rep(0,3)
   obs=objects(pos=1)
   for(i in 1:3)
   {
    if(sum(obs==namess2[i])){func.list2[i]=1}
   }

   func.list.timehomo=func.list2*0

   for(i in which(func.list2==1))
   {
     # which expressions vary over time
     result=eval(body(namess2[i]))
     func.list.timehomo[i]=2-(sum(diff(result)==0)==(length(result)-1))
   }
   if(any(func.list.timehomo==2)){homo=F}
   #print(func.list.timehomo)
   #BUILD ODE --------------------------------------------------------------------
   dims=rep('(',2)
   for(i in 1:2)
   {
    for(j in which(func.list2==1))
    {
      if(MAT[i,j]!='')
      {
        dims[i]=paste0(dims[i],'+(',body(namess2[j])[2],')',c('*','%')[func.list.timehomo[j]],MAT[i,j])
      }
    }
    dims[i]=paste0(dims[i],')')
   }

   if(any(dims=='()')){dims[which(dims=='()')]='0*a.col(0)'}


   for(i in 1:2)
   {
      dims[i]=paste0(paste0(paste0('   atemp.col(',i-1,')='),dims[i]),';')
   }
   # WRIGHT AND SOURCE -----------------------------------------------------------

   txt.full=paste(fpart,'\n',dims[1],'\n',dims[2],ODEpart,Dpart)
   type.sol ="                 Generalized Ornstein-Uhlenbeck "
   #library(Rcpp)
   #library(RcppArmadillo)
   if(wrt){write(txt.full,file='GQD.mcmc.cpp')}
   stre="Compiling C++ code. Please wait."
   cat(stre, " \r")
   flush.console()
   sourceCpp(code=txt.full)
   cat('                                     ','\r')
   }

   if(state2)
   {
   # DATA RESOLUTION -------------------------------------------------------------
   if((!homo.res)&(TR.order==4))
   {
    delt=cbind(diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh)
    if(!is.null(exclude))
    {
      diffs=diff(T.seq)/mesh
      diffs[excl]=1/20/mesh
      delt=cbind(diffs,diffs,diffs,diffs)
    }
   }
   if((!homo.res)&(TR.order==6))
   {
    delt=cbind(diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh)
    if(!is.null(exclude))
    {
      diffs=diff(T.seq)/mesh
      diffs[excl]=1/20/mesh
      delt=cbind(diffs,diffs,diffs,diffs,diffs,diffs)
    }
   }
   if((!homo.res)&(TR.order==8))
   {
     delt=cbind(diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh)
     if(!is.null(exclude))
     {
       diffs=diff(T.seq)/mesh
       diffs[excl]=1/20/mesh
       delt=cbind(diffs,diffs,diffs,diffs,diffs,diffs,diffs,diffs)
     }
   }
   # REFERENCE MATRIX ------------------------------------------------------------
   if(TR.order==4)
   {
     MAT=rbind(
     c('(1+0*a.col(0))','a.col(0)','(a.col(1)+a.col(0)%a.col(0))','','',''),
     c('','(2*a.col(1))','(2*a.col(2)+4*a.col(0)%a.col(1))','(1+0*a.col(0))','a.col(0)','(a.col(1)+a.col(0)%a.col(0))'),
     c('','(3*a.col(2))','(3*a.col(3)+6*a.col(0)%a.col(2)+6*a.col(1)%a.col(1))','','(3*a.col(1))','(3*a.col(2)+6*a.col(0)%a.col(1))'),
     c('','(4*a.col(3))','(8*a.col(0)%a.col(3)+24*a.col(1)%a.col(2))','','(6*a.col(2))','(6*a.col(3)+12*a.col(0)%a.col(2)+12*a.col(1)%a.col(1))'))
   }
   if(TR.order==6)
   {
      MAT = rbind(
     c('(1+0*a.col(0))','(a.col(0))'  ,'(1*a.col(1)+1*a.col(0)%a.col(0))'                                            ,'','',''),
     c('','(2*a.col(1))','(2*a.col(2)+4*a.col(0)%a.col(1))'                                            ,'(1+0*a.col(0))','(   a.col(0))','(a.col(1)+a.col(0)%a.col(0))'),
     c('','(3*a.col(2))','(3*a.col(3)+6*a.col(0)%a.col(2)+6*a.col(1)%a.col(1))'                        ,'','( 3*a.col(1))','(3*a.col(2)+6*a.col(0)%a.col(1))'),
     c('','(4*a.col(3))','(4*a.col(4)+8*a.col(0)%a.col(3)+24*a.col(1)%a.col(2))'                       ,'','( 6*a.col(2))','(6*a.col(3)+12*a.col(0)%a.col(2)+12*a.col(1)%a.col(1))'),
     c('','(5*a.col(4))','(5*a.col(5)+10*a.col(0)%a.col(4)+40*a.col(1)%a.col(3)+30*a.col(2)%a.col(2))' ,'','(10*a.col(3))','(10*a.col(4)+20*a.col(0)%a.col(3)+60*a.col(1)%a.col(2))'),
     c('','(6*a.col(5))','(          +12*a.col(0)%a.col(5)+60*a.col(1)%a.col(4)+120*a.col(2)%a.col(3))','','(15*a.col(4))','(15*a.col(5)+30*a.col(0)%a.col(4)+120*a.col(1)%a.col(3)+90*a.col(2)%a.col(2))'))
   }
   if(TR.order==8)
   {
     MAT = rbind(
       c('(1+0*a.col(0))','  (a.col(0))','(1*a.col(1) +1*a.col(0)%a.col(0))'                                                                          ,'','',''),
       c('','(2*a.col(1))','(2*a.col(2) +4*a.col(0)%a.col(1))'                                                                          ,'(1+0*a.col(0))','(   a.col(0))','(a.col(1)+a.col(0)%a.col(0))'),
       c('','(3*a.col(2))','(3*a.col(3) +6*a.col(0)%a.col(2)  +6*a.col(1)%a.col(1))'                                                  ,'','( 3*a.col(1))','( 3*a.col(2) +6*a.col(0)%a.col(1))'),
       c('','(4*a.col(3))','(4*a.col(4) +8*a.col(0)%a.col(3) +24*a.col(1)%a.col(2))'                                                  ,'','( 6*a.col(2))','( 6*a.col(3)+12*a.col(0)%a.col(2) +12*a.col(1)%a.col(1))'),
       c('','(5*a.col(4))','(5*a.col(5)+10*a.col(0)%a.col(4) +40*a.col(1)%a.col(3)   +30*a.col(2)%a.col(2))'                        ,'','(10*a.col(3))','(10*a.col(4)+20*a.col(0)%a.col(3) +60*a.col(1)%a.col(2))'),
       c('','(6*a.col(5))','(6*a.col(6)+12*a.col(0)%a.col(5) +60*a.col(1)%a.col(4)  +120*a.col(2)%a.col(3))'                        ,'','(15*a.col(4))','(15*a.col(5)+30*a.col(0)%a.col(4)+120*a.col(1)%a.col(3) +90*a.col(2)%a.col(2))'),
       c('','(7*a.col(6))','(7*a.col(7)+14*a.col(0)%a.col(6) +84*a.col(1)%a.col(5)  +210*a.col(2)%a.col(4)+140*a.col(3)%a.col(3))','','(21*a.col(5))','(21*a.col(6)+42*a.col(0)%a.col(5)+210*a.col(1)%a.col(4)+420*a.col(2)%a.col(3))'),
       c('','(8*a.col(7))','(           +16*a.col(0)%a.col(7)+112*a.col(1)%a.col(6)  +336*a.col(2)%a.col(5)+560*a.col(4)%a.col(3))','','(28*a.col(6))','(28*a.col(7)+56*a.col(0)%a.col(6)+336*a.col(1)%a.col(5)+840*a.col(2)%a.col(4)+560*a.col(3)%a.col(3))'))
   }
   
   # HOMOGENEITY -----------------------------------------------------------------
   namess2=c('G0','G1','G2','Q0','Q1','Q2')
   func.list2=rep(0,6)
   obs=objects(pos=1)
   for(i in 1:6)
   {
    if(sum(obs==namess[i])){func.list2[i]=1}
   }

   func.list.timehomo=func.list2*0

   for(i in which(func.list2==1))
   {
     # which expressions vary over time
     result=eval(body(namess2[i]))
     func.list.timehomo[i]=2-(sum(diff(result)==0)==(length(result)-1))
   }
    #print(t)
    
    #print(func.list.timehomo)
   if(any(func.list.timehomo==2)){homo=F}

   #BUILD ODE --------------------------------------------------------------------
   dims=rep('(',TR.order)
   for(i in 1:TR.order)
   {
    for(j in which(func.list2==1))
    {
      if(MAT[i,j]!='')
      {
        dims[i]=paste0(dims[i],'+(',body(namess2[j])[2],')',c('*','%')[func.list.timehomo[j]],MAT[i,j])
      }
    }
    dims[i]=paste0(dims[i],')')
   }

   if(any(dims=='()')){dims[which(dims=='()')]='0*a.col(0)'}


   for(i in 1:TR.order)
   {
      dims[i]=paste0(paste0(paste0('   atemp.col(',i-1,')='),dims[i]),';')
   }
   # WRIGHT AND SOURCE -----------------------------------------------------------


   if(TR.order==4)
   {
     txt.full=paste(fpart,'\n',dims[1],'\n',dims[2],'\n',dims[3],'\n',dims[4],ODEpart,Dpart)
   }
   if(TR.order==6)
   {
     txt.full=paste(fpart,'\n',dims[1],'\n',dims[2],'\n',dims[3],'\n',dims[4],'\n',dims[5],'\n',dims[6],ODEpart,Dpart)
   }
   if(TR.order==8)
   {
     txt.full=paste(fpart,'\n',dims[1],'\n',dims[2],'\n',dims[3],'\n',dims[4],'\n',dims[5],'\n',dims[6],'\n',dims[7],'\n',dims[8],ODEpart,Dpart)
   }
   type.sol ="                 Generalized Quadratic Diffusion (GQD) "

   #library(Rcpp)
   #library(RcppArmadillo)
   if(wrt){write(txt.full,file='GQD.mcmc.cpp')}
   
   stre="Compiling C++ code. Please wait."
   cat(stre, " \r")
   flush.console()
   sourceCpp(code=txt.full)
   cat('                                     ','\r')
   
   
   }
   #==============================================================================
   #                           Prior distribution Module
   #==============================================================================
   if(sum(objects(pos=1)=='priors')==1)
   {
     pp=function(theta){}
     body(pp)=parse(text =body(priors)[2])
     prior.list=paste0('d(theta)',':',paste0(body(priors)[2]))
   }
   if(sum(objects(pos=1)=='priors')==0) 
   {
     prior.list=paste0('d(theta)',':',' None.')
     pp=function(theta){1}
   }
   if(length(pp(theta))!=1){stop("Prior density must return only a single value!")}


   #==============================================================================
   #                           Interface Module
   #==============================================================================

    namess4=namess
    trim <- function (x) gsub("([[:space:]])", "", x)

    for(i in 1:6)
    {
        if(sum(obs==namess4[i]))
        {
          namess4[i]=paste0(namess4[i],' : ',trim(body(namess4[i])[2]))
        }
    }
    namess4=matrix(namess4,length(namess4),1)
    prior.list = trim(prior.list)
    buffer0=c('================================================================')
    buffer1=c('----------------------------------------------------------------')
    buffer2=c('................................................................')
    buffer3=c('...   ...   ...   ...   ...   ...   ...   ...   ...   ...   ... ')
    buffer4=c('_____________________ Drift Coefficients _______________________')
    buffer5=c('___________________ Diffusion Coefficients _____________________')
    buffer6=c('_____________________ Prior Distributions ______________________')
    buffer7=c('_______________________ Model/Chain Info _______________________')

    Info=c(buffer0,type.sol,buffer0,buffer4,namess4[1:3],buffer5,namess4[4:6],buffer6,'',prior.list)
    Info=data.frame(matrix(Info,length(Info),1))
    colnames(Info)=''
    if(print.output)
    {
     print(Info,row.names = FALSE,right=F)
    }



    ############################################################################
    ############################################################################
    ############################################################################

    tme=Sys.time()
    
    par.matrix=matrix(0,length(theta),updates)
    muvec  = theta
    covvec = diag(sds^2)
    errs   = rep(0,updates)
    ll     = rep(0,updates)
    acc=ll
    kk=0
    par.matrix[,1]=theta
    prop.matrix = par.matrix
    retries = 0
    max.retries = 0
    retry.count   = 0
    retry.indexes = c()
    success = TRUE


    rs=solver(X[-nnn],X[-1],c(0,theta),mesh,delt,nnn-1,T.seq[-nnn],P,alpha,lower,upper,TR.order)
    if(is.na(sum(rs$like)))
    {
          retry.count =1
          while(is.na(sum(rs$like))&&(retry.count<=10))
          {
            theta.start=theta+rnorm(length(theta),sd=sds)
            rs = solver(X[-nnn],X[-1],c(0,theta.start),mesh,delt,nnn-1,T.seq[-nnn],P,alpha,lower,upper,TR.order)
            if(is.na(sum(rs$like[-excl]))){retry.count=retry.count+1}
          }
    }

    lold=sum(rs$like)
    if(rs$max>0.1){warning('Probable that starting values were chosen poorly.')}
    errs[1] =rs$max
    ll[1]=lold
    if(adapt==0)
    {

    pb <- txtProgressBar(1,updates,1,style = 1,width = 65)

    failed.chain=F
    i=2
    while(i<=updates)
    {
        theta.temp=theta
        theta=theta+rnorm(length(theta),sd=sds)
        prop.matrix[,i] = theta
        rs = solver(X[-nnn],X[-1],c(0,theta),mesh,delt,nnn-1,T.seq[-nnn],P,alpha,lower,upper,TR.order)
        lnew=sum(rs$like[-excl])
        rat=min(exp(lnew-lold)*pp(theta)/pp(theta.temp),1)
        if(is.na(rat))
        {
          retry.count =1
          retries=retries+1
          retry.indexes[retries] = i
          max.retries=max.retries+1
          while(is.na(rat)&&(retry.count<=10))
          {
            i = max(i-10,2)
            theta = par.matrix[,i]
            theta=theta+rnorm(length(theta),sd=sds)
            prop.matrix[,i] = theta
            rs = solver(X[-nnn],X[-1],c(0,theta),mesh,delt,nnn-1,T.seq[-nnn],P,alpha,lower,upper,TR.order)
            lnew=sum(rs$like[-excl])
            rat=min(exp(lnew-lold)*pp(theta)/pp(theta.temp),1)
            if(is.na(rat)){retry.count=retry.count+1}
          }
        }


        u=runif(1)
        is.true =(rat>u)
        is.false=!is.true
        theta=theta*is.true+theta.temp*is.false
        lold=lnew*is.true +lold*is.false
        errs[i]=rs$max*is.true+errs[i-1]*is.false
        par.matrix[,i]=theta
        ll[i]=lold
        kk=kk+is.true
        acc[i]=is.true
        if(max.retries>2000){print('Fail: Failed evaluation limit exceeded!');failed.chain=T;break;}
        if(any(is.na(theta))){print('Fail: Samples were NA! ');failed.chain=T;break;}
        setTxtProgressBar(pb, i)
        i = i+1
    }
    close(pb)
    acc = cumsum(acc)/(1:updates)
    }
    if(adapt!=0)
    {
    pb <- txtProgressBar(1,updates,1,style = 1,width = 65)
    failed.chain=F
    for(i in 2:updates)
    {
        theta.temp=theta
        if(i>min(5000,round(burns/2)))
        {
          theta=theta+(1-adapt)*rnorm(length(theta),sd=sqrt(2.38^2/length(theta)*diag(covvec)))+adapt*rnorm(length(theta),sd=sqrt(0.1^2/length(theta)*diag(covvec)))
        }else
        {
          theta=theta+rnorm(length(theta),sd=sds)
        }
        prop.matrix[,i] = theta
        rs =solver(X[-nnn],X[-1],c(0,theta),mesh,delt,nnn-1,T.seq[-nnn],P,alpha,lower,upper,TR.order)
        lnew=sum(rs$like[-excl])
        rat=min(exp(lnew-lold)*pp(theta)/pp(theta.temp),1)
        u=runif(1)
        is.true =(rat>u)
        is.false=!is.true
        theta=theta*is.true+theta.temp*is.false
        lold=lnew*is.true +lold*is.false
        errs[i]=rs$max*is.true+errs[i-1]*is.false
        par.matrix[,i]=theta
        ll[i]=lold
        kk=kk+is.true
        acc[i]=kk/i
        muvec=muvec +1/(i)*(theta-muvec)
        covvec = covvec+1/(i)*((theta-muvec)%*%t(theta-muvec)-covvec)
        if(any(is.na(theta))){print('Fail');failed.chain=T;break;}
        setTxtProgressBar(pb, i)

    }
    close(pb)
    #print(sds)
    }
    tme.eval = function(start_time)
    {
      start_time = as.POSIXct(start_time)
      dt = difftime(Sys.time(), start_time, units="secs")
      format(.POSIXct(dt,tz="GMT"), "%H:%M:%S")
     }
    tme=tme.eval(tme)

    if(dim(par.matrix)[1]>1)
    {
      theta =apply(par.matrix[,-c(1:burns)],1,mean)
    }
    if(dim(par.matrix)[1]==1)
    {
      theta=mean(par.matrix[,-c(1:burns)])
    }
    meanD=mean(-2*ll[-c(1:burns)])
    rs=solver(X[-nnn],X[-1],c(0,theta),mesh,delt,nnn-1,T.seq[-nnn],P,alpha,lower,upper,TR.order)
    pd=meanD-(-2*sum(rs$like[-excl]))
    DIC=meanD+pd
    actual.p=length(theta)

    model.inf=list(elapsed.time=tme,time.homogeneous=c('Yes','No')[2-homo],p=actual.p,DIC=DIC,pd=pd,N=length(X)-length(excl)+1,Tag=Tag,burns=burns)
    Info2=c( buffer7,
             paste0("Chain Updates       : ",updates),
             paste0("Burned Updates      : ",burns),
             paste0("Time Homogeneous    : ",c('Yes','No')[2-homo]),
             paste0("Data Resolution     : ",c(paste0('Homogeneous: dt=',round(max(diff(T.seq)[-excl]),4)),paste0('Variable: min(dt)=',round(min(diff(T.seq)[-excl]),4),', max(dt)=',round(max(diff(T.seq)[-excl]),4)))[2-homo.res]),
             paste0("# Removed Transits. : ",c("None",length(excl))[2-is.null(exclude)]),
             paste0("Density approx.     : ",sol.state),
             paste0('Elapsed time        : ',tme),
             buffer3,
             paste0("dim(theta)          : ",round(actual.p,3)),
             paste0("DIC                 : ",round(DIC,3)),
             paste0("pd (eff. dim(theta)): ",round(pd,3)),
             buffer1)
    Info2=data.frame(matrix(Info2,length(Info2),1))
    colnames(Info2)=''
    if(print.output)
    {
    print(Info2,row.names = FALSE,right=F)
    }
    if(plot.chain)
    {
      nper=length(theta)
      if(nper==1){par(mfrow=c(1,2))}
      if(nper==2){par(mfrow=c(2,2))}
      if(nper==3){par(mfrow=c(2,2))}
      if(nper>3)
      {
        d1=1:((nper)+1)
        d2=d1
        O=outer(d1,d2)
        test=O-((nper)+1)
        test[test<0]=100
        test=test[1:4,1:4]
        test
        wh=which(test==min(test))
        
        d1=d1[col(test)[wh[1]]]
        d2=d2[row(test)[wh[1]]]
        par(mfrow=c(d1,d2))
      }
      if(palette=='mono')
      {
        cols =rep('#222299',nper)
      }else{
        cols=rainbow_hcl(nper, start = 10, end = 275,c=100,l=70)
      }
      ylabs=paste0('theta[',1:nper,']')
      for(i in 1:nper)
      {
          plot(prop.matrix[i,],col='gray90',type='s',main=ylabs[i],xlab='Iteration',ylab='')
          lines(par.matrix[i,],col=cols[i],type='s')
          abline(v=burns,lty='dotdash')
          if(adapt!=0){abline(v=min(5000,round(burns/2)),lty='dotted',col='red')}
      }
      plot(acc,type='l',ylim=c(0,1),col='darkblue',main='Accept. Rate',xlab='Iteration',ylab='%/100')
      abline(h=seq(0,1,1/10),lty='dotted')
      abline(v=burns,lty='dotdash')
      abline(h=0.4,lty='solid',col='red',lwd=1.2)
      abline(h=0.2,lty='solid',col='red',lwd=1.2)
      box()
    }
    ret=list(par.matrix=t(par.matrix),acceptance.rate=acc,elapsed.time=tme,model.info=model.inf,failed.chain=failed.chain,covvec=covvec,prop.matrix=t(prop.matrix),errors = errs)
    class(ret) = 'GQD.mcmc'
    return(ret)
  }