# R/BiGQD.mle.R In DiffusionRgqd: Inference and Analysis for Generalized Quadratic Diffusions

``` #rm(list=ls(all=TRUE))

{
rtf=0*rtf+1
solver   =function(Xs, Xt, theta, N , delt , N2, tt  , P , alpha, lower , upper, tro  ){}
rm(list =c('solver'))
check_for_model=function()
{
txt=''
namess=c('a00','a10','a20','a01','a02','a11',
'b00','b10','b20','b01','b02','b11',
'c00','c10','c20','c01','c02','c11',
'd00','d10','d20','d01','d02','d11',
'e00','e10','e20','e01','e02','e11',
'f00','f10','f20','f01','f02','f11')
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()
a00=function(t){theta[1]*theta[2]}
a10=function(t){-theta[1]}
c00=function(t){theta[3]*theta[3]}

b00=function(t){theta[4]*theta[5]}
b01=function(t){-theta[4]}
f00=function(t){theta[6]*theta[6]}

model=BiGQD.mle(X,time,10,theta =rep(0.1,6))
--------------------------------------------------------------------------------
'
check=T
}
if((sum(func.list)>0)&&(sum(func.list[-c(1:12)])==0))
{
txt='
--------------------------------------------------------------------------------
At least two diffusion coefficients have to be defined! Try for example:
--------------------------------------------------------------------------------
GQD.remove()
c00=function(t){theta[1]*theta[1]}
f00=function(t){theta[2]*theta[2]}
model=BiGQD.mle(X,time,10,theta =rep(0.1,2))
--------------------------------------------------------------------------------
'
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)
T.seq=time
Dindex = which(Dtypes==Dtype)

# A power function used to solve the syntax gap issue between R and C++
pow=function(x,p)
{
x^p
}
prod=function(a,b){a*b}

# Separate the time series into X and Y components
X1=X[,1]
X2=X[,2]
nnn=length(X1)

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 matrix!.\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 Saddlepoint, Edgeworth or Normal.\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 < {N.updates}.\n'
,'30. MCMC : Argument {N.updates} must be > 2.\n'
,'31. MCMC : length(theta)!=length(sds).\n'
,'32. Model: There has to be at least one model coefficient.\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('a00','a10','a20','a01','a02','a11',
'b00','b10','b20','b01','b02','b11',
'c00','c10','c20','c01','c02','c11',
'd00','d10','d20','d01','d02','d11',
'e00','e10','e20','e01','e02','e11',
'f00','f10','f20','f01','f02','f11')
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('a00','a10','a20','a01','a02','a11',
'b00','b10','b20','b01','b02','b11',
'c00','c10','c20','c01','c02','c11',
'd00','d10','d20','d01','d02','d11',
'e00','e10','e20','e01','e02','e11',
'f00','f10','f20','f01','f02','f11')
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.matrix(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(dim(X)[1]<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(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(dim(X)[1]!=length(time))                                  {warntrue[28]=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)
}
############################################################################
############################################################################
############################################################################

#==============================================================================
#                           Data resolution Module
#==============================================================================
#t=seq(0,100,by=1/100)
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}
}

#==============================================================================
#                           Solution TYPE Module
#==============================================================================
namess=c('a00','a10','a20','a01','a02','a11',
'b00','b10','b20','b01','b02','b11',
'c00','c10','c20','c01','c02','c11',
'd00','d10','d20','d01','d02','d11',
'e00','e10','e20','e01','e02','e11',
'f00','f10','f20','f01','f02','f11')
func.list=rep(0,36)
obs=objects(pos=1)
for(i in 1:36)
{
if(sum(obs==namess[i])){func.list[i]=1}
}

set1.1=c(3,14,15,20,21)                          # X Nonlinear
set1.2=c(11,28,29,34,35)                         # Y Nonlinear
set4=c(3,5,6,9,11,12,14:18,20:24,26:30,32:36)    # Nonlinear terms
set5=c(4:6,8,9,12,16:18,22:24,26,27,30,32,33,36) # Dependant terms

state3.0=(     any(func.list[set1.1]==1)     &     any(func.list[set1.2]==1)     &(sum(func.list[set5])==0))  # Any nonlinear in BOTH     and NO interaction term
state3.1=((sum(any(func.list[set1.1]==1))==0)&(    any(func.list[set1.2]==1)    )&(sum(func.list[set5])==0))  # Any nonlinear Y ,Normal X and NO interaction term
state3.2=(     any(func.list[set1.1]==1)     &(sum(any(func.list[set1.2]==1))==0)&(sum(func.list[set5])==0))  # Any nonlinear X ,Normal Y and NO interaction term
state4  =(any(func.list[set4]==1)&any(func.list[set5]==1))                                                    # Any nonlinear term AND interaction terms
state1  =!(state3.0|state3.1|state3.2|state4)
sol.state=c(state1,state3.0,state3.1,state3.2,state4)
sol.state=c('2nd Ord. Truncation, Bivariate Normal','4th Ord. Truncation, Independent Saddlepoints ','2nd & 4th Ord. Truncation (Indep)','4th & 2nd Ord. Trunc. (Indep)',"4th Ord. Truncation, Bivariate-Saddlepoint")[which(sol.state==T)]

#==============================================================================
#                    Check the syntax of each coefficient
#==============================================================================
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 overloaded 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])
}

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

tro=c(5,8,8,8,14)[which(c(state1,state3.0,state3.1,state3.2,state4)==T)]

if(tro==5)
{
txtA='
#include <math.h>
#include <Rcpp.h>
#define pi           3.14159265358979323846  /* pi */
using namespace arma;
using namespace Rcpp;
using namespace R;

// [[Rcpp::export]]
vec prod(vec a,vec b)
{
return(a%b);
}
mat f(mat a,vec theta,vec t,int N2)
{
mat atemp(N2,5);'

secmom = 3
}
if(tro==8)
{
txtA='
#include <math.h>
#include <Rcpp.h>
#define pi           3.14159265358979323846  /* pi */
using namespace arma;
using namespace Rcpp;
using namespace R;

// [[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);'

secmom = 5
}

if(tro==14)
{
txtA='
#include <math.h>
#include <Rcpp.h>
#define pi           3.14159265358979323846  /* pi */
using namespace arma;
using namespace Rcpp;
using namespace R;

// [[Rcpp::export]]
vec prod(vec a,vec b)
{
return(a%b);
}
mat f(mat a,vec theta,vec t,int N2)
{
mat atemp(N2,14);'
secmom = 5
}

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

// [[Rcpp::export]]
List  solver(vec Xs,vec Ys,vec Xt,vec Yt,vec theta,int N,double delt,int N2,vec tt,mat starts,int tro,int secmom)
{
mat resss(N2,3);
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;
x0.col(secmom-1)=Ys;
vec d=tt;
for (int i = 1; i < N; 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)
{
txtB= '
return atemp;
}

// [[Rcpp::export]]
List  solver(vec Xs,vec Ys,vec Xt,vec Yt,vec theta,int N, mat delt,int N2,vec d,mat starts,int tro,int secmom)
{
mat resss(N2,3);
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;
x0.col(secmom-1)=Ys;
for (int i = 1; i < N; 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(RK.order==10)
{
if(homo.res)
{
txtB= '
return atemp;
}

// [[Rcpp::export]]
List  solver(vec Xs,vec Ys,vec Xt,vec Yt,vec theta,int N,double delt,int N2,vec tt,mat starts,int tro,int secmom)
{
mat resss(N2,3);
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);
mat x0(N2,tro);
mat x1(N2,tro);
mat x2(N2,tro);
mat x3(N2,tro);
mat x4(N2,tro);
mat x5(N2,tro);
mat x6(N2,tro);
mat x7(N2,tro);
mat x8(N2,tro);
mat x9(N2,tro);
mat x10(N2,tro);
mat x11(N2,tro);
mat x12(N2,tro);
mat x13(N2,tro);
mat x14(N2,tro);
mat x15(N2,tro);
mat x16(N2,tro);
x0.fill(0);
x0.col(0)=Xs;
x0.col(secmom-1)=Ys;
vec d=tt;
for (int i = 1; i < N; i++)
{

fx0=f(x0,theta,d,N2);
x1=x0+delt*(0.1*fx0);
fx1=f(x1,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt,N2);
x2=x0+delt*(-0.915176561375291*fx0+1.45453440217827*fx1);
fx2=f(x2,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt,N2);
x3=x0+delt*(0.202259190301118*fx0+0.606777570903354*fx2);
fx3=f(x3,theta,d+0.809036761204472681298727796821953655285910174551513523258250*delt,N2);
x4=x0+delt*(0.184024714708644*fx0+0.197966831227192*fx2-0.0729547847313633*fx3);
fx4=f(x4,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt,N2);
x5=x0+delt*(0.0879007340206681*fx0+0.410459702520261*fx3+0.482713753678866*fx4);
fx5=f(x5,theta,d+0.981074190219795268254879548310562080489056746118724882027805*delt,N2);
x6=x0+delt*(0.085970050490246*fx0+0.330885963040722*fx3+0.48966295730945*fx4-0.0731856375070851*fx5);
fx6=f(x6,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt,N2);
x7=x0+delt*(0.120930449125334*fx0+0.260124675758296*fx4+0.0325402621549091*fx5-0.0595780211817361*fx6);
fx7=f(x7,theta,d+0.354017365856802376329264185948796742115824053807373968324184*delt,N2);
x8=x0+delt*(0.110854379580391*fx0-0.0605761488255006*fx5+0.321763705601778*fx6+0.510485725608063*fx7);
fx8=f(x8,theta,d+0.882527661964732346425501486979669075182867844268052119663791*delt,N2);
x9=x0+delt*(0.112054414752879*fx0-0.144942775902866*fx5-0.333269719096257*fx6+0.49926922955688*fx7+0.509504608929686*fx8);
fx9=f(x9,theta,d+0.642615758240322548157075497020439535959501736363212695909875*delt,N2);
x10=x0+delt*(0.113976783964186*fx0-0.0768813364203357*fx5+0.239527360324391*fx6+0.397774662368095*fx7+0.0107558956873607*fx8-0.327769124164019*fx9);
fx10=f(x10,theta,d+0.357384241759677451842924502979560464040498263636787304090125*delt,N2);
x11=x0+delt*(0.0798314528280196*fx0-0.0520329686800603*fx5-0.0576954146168549*fx6+0.194781915712104*fx7+0.145384923188325*fx8-0.0782942710351671*fx9-0.114503299361099*fx10);
fx11=f(x11,theta,d+0.117472338035267653574498513020330924817132155731947880336209*delt,N2);
x12=x0+delt*(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);
fx12=f(x12,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt,N2);
x13=x0+delt*(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);
fx13=f(x13,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt,N2);
x14=x0+delt*(-0.915176561375291*fx0+1.45453440217827*fx1+0*fx2+0*fx3-0.777333643644968*fx4+0*fx5-0.0910895662155176*fx6+0.0910895662155176*fx12+0.777333643644968*fx13);
fx14=f(x14,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt,N2);
x15=x0+delt*(0.1*fx0-0.157178665799771*fx2+0.157178665799771*fx14);
fx15=f(x15,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt,N2);
x16=x0+delt*(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);
fx16=f(x16,theta,d+delt,N2);

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)*delt;
d=d+delt;
}
'

}
if(!homo.res)
{
txtB= '
return atemp;
}

// [[Rcpp::export]]
List  solver(vec Xs,vec Ys,vec Xt,vec Yt,vec theta,int N, mat delt,int N2,vec d,mat starts,int tro,int secmom)
{
mat resss(N2,3);
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);
mat x0(N2,tro);
mat x1(N2,tro);
mat x2(N2,tro);
mat x3(N2,tro);
mat x4(N2,tro);
mat x5(N2,tro);
mat x6(N2,tro);
mat x7(N2,tro);
mat x8(N2,tro);
mat x9(N2,tro);
mat x10(N2,tro);
mat x11(N2,tro);
mat x12(N2,tro);
mat x13(N2,tro);
mat x14(N2,tro);
mat x15(N2,tro);
mat x16(N2,tro);
x0.fill(0);
x0.col(0)=Xs;
x0.col(secmom-1)=Ys;
for (int i = 1; i < N; i++)
{
fx0=f(x0,theta,d,N2);
x1=x0+delt%(0.1*fx0);
fx1=f(x1,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt.col(0),N2);
x2=x0+delt%(-0.915176561375291*fx0+1.45453440217827*fx1);
fx2=f(x2,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt.col(0),N2);
x3=x0+delt%(0.202259190301118*fx0+0.606777570903354*fx2);
fx3=f(x3,theta,d+0.809036761204472681298727796821953655285910174551513523258250*delt.col(0),N2);
x4=x0+delt%(0.184024714708644*fx0+0.197966831227192*fx2-0.0729547847313633*fx3);
fx4=f(x4,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt.col(0),N2);
x5=x0+delt%(0.0879007340206681*fx0+0.410459702520261*fx3+0.482713753678866*fx4);
fx5=f(x5,theta,d+0.981074190219795268254879548310562080489056746118724882027805*delt.col(0),N2);
x6=x0+delt%(0.085970050490246*fx0+0.330885963040722*fx3+0.48966295730945*fx4-0.0731856375070851*fx5);
fx6=f(x6,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt.col(0),N2);
x7=x0+delt%(0.120930449125334*fx0+0.260124675758296*fx4+0.0325402621549091*fx5-0.0595780211817361*fx6);
fx7=f(x7,theta,d+0.354017365856802376329264185948796742115824053807373968324184*delt.col(0),N2);
x8=x0+delt%(0.110854379580391*fx0-0.0605761488255006*fx5+0.321763705601778*fx6+0.510485725608063*fx7);
fx8=f(x8,theta,d+0.882527661964732346425501486979669075182867844268052119663791*delt.col(0),N2);
x9=x0+delt%(0.112054414752879*fx0-0.144942775902866*fx5-0.333269719096257*fx6+0.49926922955688*fx7+0.509504608929686*fx8);
fx9=f(x9,theta,d+0.642615758240322548157075497020439535959501736363212695909875*delt.col(0),N2);
x10=x0+delt%(0.113976783964186*fx0-0.0768813364203357*fx5+0.239527360324391*fx6+0.397774662368095*fx7+0.0107558956873607*fx8-0.327769124164019*fx9);
fx10=f(x10,theta,d+0.357384241759677451842924502979560464040498263636787304090125*delt.col(0),N2);
x11=x0+delt%(0.0798314528280196*fx0-0.0520329686800603*fx5-0.0576954146168549*fx6+0.194781915712104*fx7+0.145384923188325*fx8-0.0782942710351671*fx9-0.114503299361099*fx10);
fx11=f(x11,theta,d+0.117472338035267653574498513020330924817132155731947880336209*delt.col(0),N2);
x12=x0+delt%(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);
fx12=f(x12,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt.col(0),N2);
x13=x0+delt%(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);
fx13=f(x13,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt.col(0),N2);
x14=x0+delt%(-0.915176561375291*fx0+1.45453440217827*fx1+0*fx2+0*fx3-0.777333643644968*fx4+0*fx5-0.0910895662155176*fx6+0.0910895662155176*fx12+0.777333643644968*fx13);
fx14=f(x14,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt.col(0),N2);
x15=x0+delt%(0.1*fx0-0.157178665799771*fx2+0.157178665799771*fx14);
fx15=f(x15,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt.col(0),N2);
x16=x0+delt%(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);
fx16=f(x16,theta,d+delt.col(0),N2);

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)%delt;
d=d+delt.col(0);
}
'
}
}

if(state1)
{
txtC=
'
vec det=(x0.col(1)%x0.col(3)-x0.col(4)%x0.col(4));
resss.col(0)=-log(2*3.141592653589793)-0.5*log(abs(det))-0.5*((Xt-x0.col(0))%(Xt-x0.col(0))%x0.col(3)/det-(Xt-x0.col(0))%(Yt-x0.col(2))%x0.col(4)/det-(Xt-x0.col(0))%(Yt-x0.col(2))%x0.col(4)/det+(Yt-x0.col(2))%(Yt-x0.col(2))%x0.col(1)/det);
List ret;
ret["like"] = resss;
ret["max"] = whch;
return(ret);
}'
}
if(state3.0)
{
txtC=
'
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);

p=(1.0/3.0) *(3*(x0.col(7)/6.0)%x0.col(5) - pow(x0.col(6)/2.0,2))/pow(x0.col(7)/6.0,2);
q=(1.0/27.0)*(27*pow(x0.col(7)/6.0,2)%(x0.col(4)-Yt) - 9*(x0.col(7)/6.0)%(x0.col(6)/2.0)%x0.col(5) + 2*pow(x0.col(6)/2.0,3))/pow(x0.col(7)/6.0,3);
chk=pow(q,2)/4.0 + pow(p,3)/27.0;
th=-(x0.col(6)/2.0)/(3*(x0.col(7)/6.0))+pow(-q/2.0+sqrt(chk),(1.0/3.0))-pow(q/2.0+sqrt(chk),(1.0/3.0));

K =x0.col(4)%th+(x0.col(5)%th%th)/2.0+(x0.col(6)%th%th%th)/6.0 +(x0.col(7)%th%th%th%th)/24.0;
K1=x0.col(4)   +(x0.col(5)%th)       +(x0.col(6)%th%th)/2.0    +(x0.col(7)%th%th%th)/6.0;
K2=x0.col(5)   +(x0.col(6)%th)       +(x0.col(7)%th%th)/2.0;
val=val-0.5*log(2*3.141592653589793*K2)+(K-th%K1);
resss.col(0)=val;
List ret;
ret["like"] = resss;
ret["max"] = whch;
return(ret);;
}'
if(Dtype=='Normal')
{
txtC=
'
vec val=   -0.5*log(2*3.141592653589793*x0.col(5))-0.5*pow(Yt-x0.col(4),2)/x0.col(5);
val=val-0.5*log(2*3.141592653589793*x0.col(1))-0.5*pow(Xt-x0.col(0),2)/x0.col(1);
resss.col(0)=val;
List ret;
ret["like"] = resss;
ret["max"] = whch;
return(ret);
}'
}
}
if(state3.1)
{
txtC=
'
vec p=(1.0/3.0) *(3*(x0.col(7)/6.0)%x0.col(5) - pow(x0.col(6)/2.0,2))/pow(x0.col(7)/6.0,2);
vec q=(1.0/27.0)*(27*pow(x0.col(7)/6.0,2)%(x0.col(4)-Yt) - 9*(x0.col(7)/6.0)%(x0.col(6)/2.0)%x0.col(5) + 2*pow(x0.col(6)/2.0,3))/pow(x0.col(7)/6.0,3);
vec chk=pow(q,2)/4.0 + pow(p,3)/27.0;
vec th=-(x0.col(6)/2.0)/(3*(x0.col(7)/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(4)%th+(x0.col(5)%th%th)/2.0+(x0.col(6)%th%th%th)/6.0 +(x0.col(7)%th%th%th%th)/24.0;
vec K1=x0.col(4)   +(x0.col(5)%th)       +(x0.col(6)%th%th)/2.0    +(x0.col(7)%th%th%th)/6.0;
vec K2=x0.col(5)   +(x0.col(6)%th)       +(x0.col(7)%th%th)/2.0;
vec val=-0.5*log(2*3.141592653589793*K2)+(K-th%K1);

val=val-0.5*log(2*3.141592653589793*x0.col(1))-0.5*pow(Xt-x0.col(0),2)/x0.col(1);
resss.col(0)=val;
List ret;
ret["like"] = resss;
ret["max"] = whch;
return(ret);
}'
if(Dtype=='Normal')
{

txtC='
vec val=   -0.5*log(2*3.141592653589793*x0.col(5))-0.5*pow(Yt-x0.col(4),2)/x0.col(5);
val=val-0.5*log(2*3.141592653589793*x0.col(1))-0.5*pow(Xt-x0.col(0),2)/x0.col(1);
resss.col(0)=val;
List ret;
ret["like"] = resss;
ret["max"] = whch;
return(ret);
}'
}
}
if(state3.2)
{
txtC='
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);

val=val-0.5*log(2*3.141592653589793*x0.col(5))-(0.5*pow(Yt-x0.col(4),2)/x0.col(5));
resss.col(0)=val;
List ret;
ret["like"] = resss;
ret["max"] = whch;
return(ret);
}'
if(Dtype=='Normal')
{

txtC='
vec val=   -0.5*log(2*3.141592653589793*x0.col(5))-0.5*pow(Yt-x0.col(4),2)/x0.col(5);
val=val-0.5*log(2*3.141592653589793*x0.col(1))-0.5*pow(Xt-x0.col(0),2)/x0.col(1);
resss.col(0)=val;
List ret;
ret["like"] = resss;
ret["max"] = whch;
return(ret);
}'
}
}
if(state4)
{
{
txtC='
vec a(N2);
vec b(N2);
vec abser(N2);
abser=0.1+abser;
a.ones();
b.ones();
a=starts.col(0);
b=starts.col(1);
vec gg(N2);
vec hh(N2);
vec gg1(N2);
vec hh1(N2);
vec gg2(N2);
vec hh2(N2);
vec ar(N2);
vec br(N2);
vec anew(N2);
vec bnew(N2);
int ind=0;
while((max(abser)>0.001)&&(ind<5000))
{
gg=x0.col(0)+x0.col(1)%a+(1.0/2.0)*x0.col(2)%a%a+(1.0/6.0)*x0.col(3)%a%a%a +x0.col(8)%b +(1.0/2.0)*x0.col(9)%b%b+x0.col(10)%a%b+(1.0/6.0)*b%b%b%x0.col(12)+(1.0/2.0)*a%a%b%x0.col(13)+(1.0/2.0)*a%b%b%x0.col(11)-Xt;
hh=x0.col(4)+x0.col(5)%b+(1.0/2.0)*x0.col(6)%b%b+(1.0/6.0)*x0.col(7)%b%b%b +x0.col(8)%a +x0.col(9)%a%b+(1.0/2.0)*x0.col(10)%a%a+(1.0/2.0)*a%b%b%x0.col(12)+(1.0/6.0)*a%a%a%x0.col(13)+(1.0/2.0)*a%a%b%x0.col(11)-Yt;

gg1=x0.col(1)+x0.col(2)%a+(1.0/2.0)*x0.col(3)%a%a+x0.col(10)%b+a%b%x0.col(13)+(1.0/2.0)*b%b%x0.col(11);
gg2=x0.col(8) +x0.col(9)%b+x0.col(10)%a+(1.0/2.0)*b%b%x0.col(12)+(1.0/2.0)*a%a%x0.col(13)+a%b%x0.col(11);

hh1=x0.col(8) +x0.col(9)%b+x0.col(10)%a+(1.0/2.0)*b%b%x0.col(12)+(1.0/2.0)*a%a%x0.col(13)+a%b%x0.col(11);
hh2=x0.col(5)+x0.col(6)%b+(1.0/2.0)*x0.col(7)%b%b +x0.col(9)%a+a%b%x0.col(12)+(1.0/2.0)*a%a%x0.col(11);

anew  =a+(hh%gg2-gg%hh2)/(gg1%hh2-gg2%hh1);
bnew  =b-(hh%gg1-gg%hh1)/(gg1%hh2-gg2%hh1);
abser =(pow(anew-a,2)+pow(bnew-b,2));
a=anew;
b=bnew;
ind=ind+1;
}

resss.col(0)=-log(2*3.141592653589793)-0.5*log(gg1%hh2-gg2%hh1)+(x0.col(0)%a+x0.col(4)%b+(1.0/2.0)*x0.col(1)%a%a+(1.0/2.0)*x0.col(5)%b%b+(1.0/6.0)*x0.col(2)%a%a%a+(1.0/6.0)*x0.col(6)%b%b%b+(1.0/24.0)*x0.col(3)%a%a%a%a+(1.0/24.0)*x0.col(7)%b%b%b%b+x0.col(8)%a%b+(1.0/2.0)*x0.col(9)%a%b%b+(1.0/2.0)*x0.col(10)%a%a%b+(1.0/6.0)*a%b%b%b%x0.col(12)+(1.0/6.0)*a%a%a%b%x0.col(13)+(1.0/4.0)*a%a%b%b%x0.col(11)-a%Xt-b%Yt);
resss.col(1)=a;
resss.col(2)=b;
List ret;
ret["like"] = resss;
ret["max"] = whch;
return(ret);
}'
}
if(Dtype=='Edgeworth')
{
txtC='
vec x = (Xt-x0.col(0))/sqrt(x0.col(1));
vec y = (Yt-x0.col(4))/sqrt(x0.col(5));

vec rho30 =x0.col(2)/(pow(x0.col(1),1.5));
vec rho40 =x0.col(3)/(pow(x0.col(1),2));
vec rho03 =x0.col(6)/(pow(x0.col(5),1.5));
vec rho04 =x0.col(7)/(pow(x0.col(5),2));
vec rho11 =x0.col(8)/(sqrt(x0.col(1))%sqrt(x0.col(5)));
vec rho12 =x0.col(9)/(sqrt(x0.col(1))%x0.col(5));
vec rho21 =x0.col(10)/(x0.col(1)%sqrt(x0.col(5)));
vec rho22 =x0.col(11)/(x0.col(1)%x0.col(5));
vec rho13 =x0.col(12)/(sqrt(x0.col(1))%pow(x0.col(5),1.5));
vec rho31 =x0.col(13)/(pow(x0.col(1),1.5)%sqrt(x0.col(5)));

vec Hx1 = x;
vec Hx2 = pow(x,2)-1;
vec Hx3 = pow(x,3)-3*x;
vec Hx4 = pow(x,4)-6*pow(x,2)+3;
vec Hx5 = pow(x,5)-10*pow(x,3)+15*x;
vec Hx6 = pow(x,6)-15*pow(x,4)+45*pow(x,2)-15;

vec Hy1 = y;
vec Hy2 = pow(y,2)-1;
vec Hy3 = pow(y,3)-3*y;
vec Hy4 = pow(y,4)-6*pow(y,2)+3;
vec Hy5 = pow(y,5)-10*pow(y,3)+15*y;
vec Hy6 = pow(y,6)-15*pow(y,4)+45*pow(y,2)-15;

vec A=Hx3%rho30+3*Hx2%Hy1%rho21+3*Hx1%Hy2%rho12+Hy3%rho03;
vec B=Hx4%rho40+4*Hx3%Hy1%rho31+6*Hx2%Hy2%rho22+4*Hx1%Hy3%rho13+Hy4%rho04;
vec C=Hx6%pow(rho30,2)+6*Hx5%Hy1%rho21%rho30+6*Hx4%Hy2%rho12%rho30+2*Hx3%Hy3%rho03%rho30+9*Hx4%Hy2%pow(rho21,2)+18*Hx3%Hy3%rho12%rho21+6*Hx2%Hy4%rho03%rho21+9*Hx2%Hy4%pow(rho12,2)+6*Hx1%Hy5%rho03%rho12+Hy6%pow(rho03,2);

resss.col(0)=-log(2*pi*sqrt(x0.col(1)%x0.col(5))%sqrt(1 - pow(rho11,2)))+(-(pow((Xt-x0.col(0)),2)/x0.col(1) -2*rho11%(Xt-x0.col(0))%(Yt-x0.col(4))/sqrt(x0.col(1)%x0.col(5))+pow((Yt-x0.col(4)),2)/x0.col(5))/(2*(1-pow(rho11,2))))+log(abs(1+0.1666667*A+0.04166667*B+0.01388889*C));

List ret;
ret["like"] = resss;
ret["max"] = whch;
return(ret);
}'
}
if(Dtype=='Normal')
{

txtC='
vec det=(x0.col(1)%x0.col(5)-x0.col(8)%x0.col(8));
resss.col(0)=-log(2*3.141592653589793)-0.5*log(abs(det))-0.5*((Xt-x0.col(0))%(Xt-x0.col(0))%x0.col(5)/det-(Xt-x0.col(0))%(Yt-x0.col(4))%x0.col(8)/det-(Xt-x0.col(0))%(Yt-x0.col(4))%x0.col(8)/det+(Yt-x0.col(4))%(Yt-x0.col(4))%x0.col(1)/det);
List ret;
ret["like"] = resss;
ret["max"] = whch;
return(ret);
}'
}
}

#==============================================================================
#                   Generate TYPE of Solution
#==============================================================================

if(state1)
{
# DATA RESOLUTION -------------------------------------------------------------
if(!homo.res)
{
delt=cbind(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)
}
}
# REFERENCE MATRIX ------------------------------------------------------------
MAT=rbind(
c('','','(+1*a.col(0))','(+1*a.col(2))','','','','','',''),
c('','','(+2*a.col(1))','(+2*a.col(4))','','','','','',''),
c('','','','','(+1*a.col(0))','(+1*a.col(2))','','','',''),
c('','','','','(+2*a.col(4))','(+2*a.col(3))','','','',''),
c('','','(+1*a.col(4))','(+1*a.col(3))','(+1*a.col(1))','(+1*a.col(4))','','','',''))
MAT[1,1]='(1+0*a.col(0))'        #ONES?
MAT[3,2]='(1+0*a.col(0))'
MAT[2,7]='(1+0*a.col(0))'
MAT[4,10]='(1+0*a.col(0))'
MAT[5,8]='(1+0*a.col(0))'
MAT[5,9]='(1+0*a.col(0))'

# HOMOGENEITY -----------------------------------------------------------------
namess2=c('a00','b00','a10','a01','b10','b01','c00','d00','e00','f00')
func.list2=rep(0,10)
obs=objects(pos=1)
for(i in 1:10)
{
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}

#func.list.timehomo[c(1:2,7:10)]=1    # Always set these to 1

#BUILD ODE --------------------------------------------------------------------
dims=rep('(',5)
for(i in 1:5)
{
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=='()')]='(1+0*a.col(0))'}

for(i in 1:5)
{
dims[i]=paste0(paste0(paste0('   atemp.col(',i-1,')='),dims[i]),';')
}

# WRIGHT AND SOURCE -----------------------------------------------------------
txt.full=paste(txtA,'\n',dims[1],'\n',dims[2],'\n',dims[3],'\n',dims[4],'\n',dims[5],txtB,txtC)
type.sol ="                  GENERALIZED LINEAR DIFFUSON"
}
if(any(c(state3.0,state3.1,state3.2)))
{
# DATA RESOLUTION -------------------------------------------------------------
if(!homo.res)
{
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 ------------------------------------------------------------
MAT=rbind(
c("()","(+1*a.col(0))","(+1*a.col(0)%a.col(0)+1*a.col(1))","(+1*a.col(4))","(+1*a.col(4)%a.col(4)+1*a.col(5))","(+1*a.col(8)+1*a.col(0)%a.col(4))","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()"
),c("()","(+2*a.col(1))","(+2*a.col(0)%a.col(1)+2*a.col(1)%a.col(0)+2*a.col(2))","(+2*a.col(8))","(+2*a.col(4)%a.col(8)+2*a.col(8)%a.col(4)+2*a.col(9))","(+2*a.col(10)+2*a.col(0)%a.col(8)+2*a.col(1)%a.col(4))","()","()","()","()","()","()","()","(+1*a.col(0))","(+1*a.col(1)+1*a.col(0)%a.col(0))","(+1*a.col(4))","(+1*a.col(5)+1*a.col(4)%a.col(4))","(+1*a.col(8)+1*a.col(4)%a.col(0))","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()"
),c("()","(+3*a.col(2))","(+3*a.col(0)%a.col(2)+6*a.col(1)%a.col(1)+3*a.col(2)%a.col(0)+3*a.col(3))","(+3*a.col(10))","(+3*a.col(4)%a.col(10)+6*a.col(8)%a.col(8)+3*a.col(10)%a.col(4)+3*a.col(11))","(+3*a.col(13)+3*a.col(0)%a.col(10)+6*a.col(1)%a.col(8)+3*a.col(2)%a.col(4))","()","()","()","()","()","()","()","(+3*a.col(1))","(+3*a.col(2)+3*a.col(0)%a.col(1)+3*a.col(1)%a.col(0))","(+3*a.col(8))","(+3*a.col(9)+3*a.col(4)%a.col(8)+3*a.col(8)%a.col(4))","(+3*a.col(10)+3*a.col(4)%a.col(1)+3*a.col(8)%a.col(0))","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()"
),c("()","(+4*a.col(3))","(+4*a.col(0)%a.col(3)+12*a.col(1)%a.col(2)+12*a.col(2)%a.col(1)+4*a.col(3)%a.col(0))","(+4*a.col(13))","(+4*a.col(4)%a.col(13)+12*a.col(8)%a.col(10)+12*a.col(10)%a.col(8)+4*a.col(13)%a.col(4))","(+4*a.col(0)%a.col(13)+12*a.col(1)%a.col(10)+12*a.col(2)%a.col(8)+4*a.col(3)%a.col(4))","()","()","()","()","()","()","()","(+6*a.col(2))","(+6*a.col(3)+6*a.col(0)%a.col(2)+12*a.col(1)%a.col(1)+6*a.col(2)%a.col(0))","(+6*a.col(10))","(+6*a.col(11)+6*a.col(4)%a.col(10)+12*a.col(8)%a.col(8)+6*a.col(10)%a.col(4))","(+6*a.col(13)+6*a.col(4)%a.col(2)+12*a.col(8)%a.col(1)+6*a.col(10)%a.col(0))","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()"
),c("()","()","()","()","()","()","()","(+1*a.col(0))","(+1*a.col(0)%a.col(0)+1*a.col(1))","(+1*a.col(4))","(+1*a.col(4)%a.col(4)+1*a.col(5))","(+1*a.col(8)+1*a.col(4)%a.col(0))","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()"
),c("()","()","()","()","()","()","()","(+2*a.col(8))","(+2*a.col(0)%a.col(8)+2*a.col(8)%a.col(0)+2*a.col(10))","(+2*a.col(5))","(+2*a.col(4)%a.col(5)+2*a.col(5)%a.col(4)+2*a.col(6))","(+2*a.col(9)+2*a.col(4)%a.col(8)+2*a.col(5)%a.col(0))","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","(+1*a.col(0))","(+1*a.col(1)+1*a.col(0)%a.col(0))","(+1*a.col(4))","(+1*a.col(5)+1*a.col(4)%a.col(4))","(+1*a.col(8)+1*a.col(4)%a.col(0))"
),c("()","()","()","()","()","()","()","(+3*a.col(9))","(+3*a.col(0)%a.col(9)+6*a.col(8)%a.col(8)+3*a.col(9)%a.col(0)+3*a.col(11))","(+3*a.col(6))","(+3*a.col(4)%a.col(6)+6*a.col(5)%a.col(5)+3*a.col(6)%a.col(4)+3*a.col(7))","(+3*a.col(12)+3*a.col(4)%a.col(9)+6*a.col(5)%a.col(8)+3*a.col(6)%a.col(0))","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","(+3*a.col(8))","(+3*a.col(10)+3*a.col(0)%a.col(8)+3*a.col(8)%a.col(0))","(+3*a.col(5))","(+3*a.col(6)+3*a.col(4)%a.col(5)+3*a.col(5)%a.col(4))","(+3*a.col(9)+3*a.col(4)%a.col(8)+3*a.col(5)%a.col(0))"
),c("()","()","()","()","()","()","()","(+4*a.col(12))","(+4*a.col(0)%a.col(12)+12*a.col(8)%a.col(9)+12*a.col(9)%a.col(8)+4*a.col(12)%a.col(0))","(+4*a.col(7))","(+4*a.col(4)%a.col(7)+12*a.col(5)%a.col(6)+12*a.col(6)%a.col(5)+4*a.col(7)%a.col(4))","(+4*a.col(4)%a.col(12)+12*a.col(5)%a.col(9)+12*a.col(6)%a.col(8)+4*a.col(7)%a.col(0))","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","(+6*a.col(9))","(+6*a.col(11)+6*a.col(0)%a.col(9)+12*a.col(8)%a.col(8)+6*a.col(9)%a.col(0))","(+6*a.col(6))","(+6*a.col(7)+6*a.col(4)%a.col(6)+12*a.col(5)%a.col(5)+6*a.col(6)%a.col(4))","(+6*a.col(12)+6*a.col(4)%a.col(9)+12*a.col(5)%a.col(8)+6*a.col(6)%a.col(0))"
),c("()","(+1*a.col(8))","(+1*a.col(0)%a.col(8)+1*a.col(8)%a.col(0)+1*a.col(10))","(+1*a.col(5))","(+1*a.col(4)%a.col(5)+1*a.col(5)%a.col(4)+1*a.col(6))","(+1*a.col(9)+1*a.col(0)%a.col(5)+1*a.col(8)%a.col(4))","()","(+1*a.col(1))","(+1*a.col(0)%a.col(1)+1*a.col(1)%a.col(0)+1*a.col(2))","(+1*a.col(8))","(+1*a.col(4)%a.col(8)+1*a.col(8)%a.col(4)+1*a.col(9))","(+1*a.col(10)+1*a.col(4)%a.col(1)+1*a.col(8)%a.col(0))","()","()","()","()","()","()","()","(+0.5*a.col(0))","(+0.5*a.col(1)+0.5*a.col(0)%a.col(0))","(+0.5*a.col(4))","(+0.5*a.col(5)+0.5*a.col(4)%a.col(4))","(+0.5*a.col(8)+0.5*a.col(4)%a.col(0))","()","(+0.5*a.col(0))","(+0.5*a.col(1)+0.5*a.col(0)%a.col(0))","(+0.5*a.col(4))","(+0.5*a.col(5)+0.5*a.col(4)%a.col(4))","(+0.5*a.col(8)+0.5*a.col(4)%a.col(0))","()","()","()","()","()","()"
),c("()","(+1*a.col(9))","(+1*a.col(0)%a.col(9)+2*a.col(8)%a.col(8)+1*a.col(9)%a.col(0)+1*a.col(11))","(+1*a.col(6))","(+1*a.col(4)%a.col(6)+2*a.col(5)%a.col(5)+1*a.col(6)%a.col(4)+1*a.col(7))","(+1*a.col(12)+1*a.col(0)%a.col(6)+2*a.col(8)%a.col(5)+1*a.col(9)%a.col(4))","()","(+2*a.col(10))","(+2*a.col(0)%a.col(10)+2*a.col(8)%a.col(1)+2*a.col(1)%a.col(8)+2*a.col(10)%a.col(0)+2*a.col(13))","(+2*a.col(9))","(+2*a.col(4)%a.col(9)+2*a.col(5)%a.col(8)+2*a.col(8)%a.col(5)+2*a.col(9)%a.col(4)+2*a.col(12))","(+2*a.col(11)+2*a.col(4)%a.col(10)+2*a.col(5)%a.col(1)+2*a.col(8)%a.col(8)+2*a.col(9)%a.col(0))","()","()","()","()","()","()","()","(+1*a.col(8))","(+1*a.col(10)+1*a.col(0)%a.col(8)+1*a.col(8)%a.col(0))","(+1*a.col(5))","(+1*a.col(6)+1*a.col(4)%a.col(5)+1*a.col(5)%a.col(4))","(+1*a.col(9)+1*a.col(4)%a.col(8)+1*a.col(5)%a.col(0))","()","(+1*a.col(8))","(+1*a.col(10)+1*a.col(0)%a.col(8)+1*a.col(8)%a.col(0))","(+1*a.col(5))","(+1*a.col(6)+1*a.col(4)%a.col(5)+1*a.col(5)%a.col(4))","(+1*a.col(9)+1*a.col(4)%a.col(8)+1*a.col(5)%a.col(0))","()","(+1*a.col(1))","(+1*a.col(2)+1*a.col(0)%a.col(1)+1*a.col(1)%a.col(0))","(+1*a.col(8))","(+1*a.col(9)+1*a.col(4)%a.col(8)+1*a.col(8)%a.col(4))","(+1*a.col(10)+1*a.col(4)%a.col(1)+1*a.col(8)%a.col(0))"
),c("()","(+2*a.col(10))","(+2*a.col(0)%a.col(10)+2*a.col(8)%a.col(1)+2*a.col(1)%a.col(8)+2*a.col(10)%a.col(0)+2*a.col(13))","(+2*a.col(9))","(+2*a.col(4)%a.col(9)+2*a.col(5)%a.col(8)+2*a.col(8)%a.col(5)+2*a.col(9)%a.col(4)+2*a.col(12))","(+2*a.col(11)+2*a.col(0)%a.col(9)+2*a.col(8)%a.col(8)+2*a.col(1)%a.col(5)+2*a.col(10)%a.col(4))","()","(+1*a.col(2))","(+1*a.col(0)%a.col(2)+2*a.col(1)%a.col(1)+1*a.col(2)%a.col(0)+1*a.col(3))","(+1*a.col(10))","(+1*a.col(4)%a.col(10)+2*a.col(8)%a.col(8)+1*a.col(10)%a.col(4)+1*a.col(11))","(+1*a.col(13)+1*a.col(4)%a.col(2)+2*a.col(8)%a.col(1)+1*a.col(10)%a.col(0))","()","(+1*a.col(8))","(+1*a.col(10)+1*a.col(0)%a.col(8)+1*a.col(8)%a.col(0))","(+1*a.col(5))","(+1*a.col(6)+1*a.col(4)%a.col(5)+1*a.col(5)%a.col(4))","(+1*a.col(9)+1*a.col(4)%a.col(8)+1*a.col(5)%a.col(0))","()","(+1*a.col(1))","(+1*a.col(2)+1*a.col(0)%a.col(1)+1*a.col(1)%a.col(0))","(+1*a.col(8))","(+1*a.col(9)+1*a.col(4)%a.col(8)+1*a.col(8)%a.col(4))","(+1*a.col(10)+1*a.col(4)%a.col(1)+1*a.col(8)%a.col(0))","()","(+1*a.col(1))","(+1*a.col(2)+1*a.col(0)%a.col(1)+1*a.col(1)%a.col(0))","(+1*a.col(8))","(+1*a.col(9)+1*a.col(4)%a.col(8)+1*a.col(8)%a.col(4))","(+1*a.col(10)+1*a.col(4)%a.col(1)+1*a.col(8)%a.col(0))","()","()","()","()","()","()"
),c("()","(+2*a.col(11))","(+2*a.col(0)%a.col(11)+4*a.col(8)%a.col(10)+2*a.col(9)%a.col(1)+2*a.col(1)%a.col(9)+4*a.col(10)%a.col(8)+2*a.col(11)%a.col(0))","(+2*a.col(12))","(+2*a.col(4)%a.col(12)+4*a.col(5)%a.col(9)+2*a.col(6)%a.col(8)+2*a.col(8)%a.col(6)+4*a.col(9)%a.col(5)+2*a.col(12)%a.col(4))","(+2*a.col(0)%a.col(12)+4*a.col(8)%a.col(9)+2*a.col(9)%a.col(8)+2*a.col(1)%a.col(6)+4*a.col(10)%a.col(5)+2*a.col(11)%a.col(4))","()","(+2*a.col(13))","(+2*a.col(0)%a.col(13)+2*a.col(8)%a.col(2)+4*a.col(1)%a.col(10)+4*a.col(10)%a.col(1)+2*a.col(2)%a.col(8)+2*a.col(13)%a.col(0))","(+2*a.col(11))","(+2*a.col(4)%a.col(11)+2*a.col(5)%a.col(10)+4*a.col(8)%a.col(9)+4*a.col(9)%a.col(8)+2*a.col(10)%a.col(5)+2*a.col(11)%a.col(4))","(+2*a.col(4)%a.col(13)+2*a.col(5)%a.col(2)+4*a.col(8)%a.col(10)+4*a.col(9)%a.col(1)+2*a.col(10)%a.col(8)+2*a.col(11)%a.col(0))","()","(+1*a.col(9))","(+1*a.col(11)+1*a.col(0)%a.col(9)+2*a.col(8)%a.col(8)+1*a.col(9)%a.col(0))","(+1*a.col(6))","(+1*a.col(7)+1*a.col(4)%a.col(6)+2*a.col(5)%a.col(5)+1*a.col(6)%a.col(4))","(+1*a.col(12)+1*a.col(4)%a.col(9)+2*a.col(5)%a.col(8)+1*a.col(6)%a.col(0))","()","(+2*a.col(10))","(+2*a.col(13)+2*a.col(0)%a.col(10)+2*a.col(8)%a.col(1)+2*a.col(1)%a.col(8)+2*a.col(10)%a.col(0))","(+2*a.col(9))","(+2*a.col(12)+2*a.col(4)%a.col(9)+2*a.col(5)%a.col(8)+2*a.col(8)%a.col(5)+2*a.col(9)%a.col(4))","(+2*a.col(11)+2*a.col(4)%a.col(10)+2*a.col(5)%a.col(1)+2*a.col(8)%a.col(8)+2*a.col(9)%a.col(0))","()","(+2*a.col(10))","(+2*a.col(13)+2*a.col(0)%a.col(10)+2*a.col(8)%a.col(1)+2*a.col(1)%a.col(8)+2*a.col(10)%a.col(0))","(+2*a.col(9))","(+2*a.col(12)+2*a.col(4)%a.col(9)+2*a.col(5)%a.col(8)+2*a.col(8)%a.col(5)+2*a.col(9)%a.col(4))","(+2*a.col(11)+2*a.col(4)%a.col(10)+2*a.col(5)%a.col(1)+2*a.col(8)%a.col(8)+2*a.col(9)%a.col(0))","()","(+1*a.col(2))","(+1*a.col(3)+1*a.col(0)%a.col(2)+2*a.col(1)%a.col(1)+1*a.col(2)%a.col(0))","(+1*a.col(10))","(+1*a.col(11)+1*a.col(4)%a.col(10)+2*a.col(8)%a.col(8)+1*a.col(10)%a.col(4))","(+1*a.col(13)+1*a.col(4)%a.col(2)+2*a.col(8)%a.col(1)+1*a.col(10)%a.col(0))"
),c("()","(+1*a.col(12))","(+1*a.col(0)%a.col(12)+3*a.col(8)%a.col(9)+3*a.col(9)%a.col(8)+1*a.col(12)%a.col(0))","(+1*a.col(7))","(+1*a.col(4)%a.col(7)+3*a.col(5)%a.col(6)+3*a.col(6)%a.col(5)+1*a.col(7)%a.col(4))","(+1*a.col(0)%a.col(7)+3*a.col(8)%a.col(6)+3*a.col(9)%a.col(5)+1*a.col(12)%a.col(4))","()","(+3*a.col(11))","(+3*a.col(0)%a.col(11)+6*a.col(8)%a.col(10)+3*a.col(9)%a.col(1)+3*a.col(1)%a.col(9)+6*a.col(10)%a.col(8)+3*a.col(11)%a.col(0))","(+3*a.col(12))","(+3*a.col(4)%a.col(12)+6*a.col(5)%a.col(9)+3*a.col(6)%a.col(8)+3*a.col(8)%a.col(6)+6*a.col(9)%a.col(5)+3*a.col(12)%a.col(4))","(+3*a.col(4)%a.col(11)+6*a.col(5)%a.col(10)+3*a.col(6)%a.col(1)+3*a.col(8)%a.col(9)+6*a.col(9)%a.col(8)+3*a.col(12)%a.col(0))","()","()","()","()","()","()","()","(+1.5*a.col(9))","(+1.5*a.col(11)+1.5*a.col(0)%a.col(9)+3*a.col(8)%a.col(8)+1.5*a.col(9)%a.col(0))","(+1.5*a.col(6))","(+1.5*a.col(7)+1.5*a.col(4)%a.col(6)+3*a.col(5)%a.col(5)+1.5*a.col(6)%a.col(4))","(+1.5*a.col(12)+1.5*a.col(4)%a.col(9)+3*a.col(5)%a.col(8)+1.5*a.col(6)%a.col(0))","()","(+1.5*a.col(9))","(+1.5*a.col(11)+1.5*a.col(0)%a.col(9)+3*a.col(8)%a.col(8)+1.5*a.col(9)%a.col(0))","(+1.5*a.col(6))","(+1.5*a.col(7)+1.5*a.col(4)%a.col(6)+3*a.col(5)%a.col(5)+1.5*a.col(6)%a.col(4))","(+1.5*a.col(12)+1.5*a.col(4)%a.col(9)+3*a.col(5)%a.col(8)+1.5*a.col(6)%a.col(0))","()","(+3*a.col(10))","(+3*a.col(13)+3*a.col(0)%a.col(10)+3*a.col(8)%a.col(1)+3*a.col(1)%a.col(8)+3*a.col(10)%a.col(0))","(+3*a.col(9))","(+3*a.col(12)+3*a.col(4)%a.col(9)+3*a.col(5)%a.col(8)+3*a.col(8)%a.col(5)+3*a.col(9)%a.col(4))","(+3*a.col(11)+3*a.col(4)%a.col(10)+3*a.col(5)%a.col(1)+3*a.col(8)%a.col(8)+3*a.col(9)%a.col(0))"
),c("()","(+3*a.col(13))","(+3*a.col(0)%a.col(13)+3*a.col(8)%a.col(2)+6*a.col(1)%a.col(10)+6*a.col(10)%a.col(1)+3*a.col(2)%a.col(8)+3*a.col(13)%a.col(0))","(+3*a.col(11))","(+3*a.col(4)%a.col(11)+3*a.col(5)%a.col(10)+6*a.col(8)%a.col(9)+6*a.col(9)%a.col(8)+3*a.col(10)%a.col(5)+3*a.col(11)%a.col(4))","(+3*a.col(0)%a.col(11)+3*a.col(8)%a.col(10)+6*a.col(1)%a.col(9)+6*a.col(10)%a.col(8)+3*a.col(2)%a.col(5)+3*a.col(13)%a.col(4))","()","(+1*a.col(3))","(+1*a.col(0)%a.col(3)+3*a.col(1)%a.col(2)+3*a.col(2)%a.col(1)+1*a.col(3)%a.col(0))","(+1*a.col(13))","(+1*a.col(4)%a.col(13)+3*a.col(8)%a.col(10)+3*a.col(10)%a.col(8)+1*a.col(13)%a.col(4))","(+1*a.col(4)%a.col(3)+3*a.col(8)%a.col(2)+3*a.col(10)%a.col(1)+1*a.col(13)%a.col(0))","()","(+3*a.col(10))","(+3*a.col(13)+3*a.col(0)%a.col(10)+3*a.col(8)%a.col(1)+3*a.col(1)%a.col(8)+3*a.col(10)%a.col(0))","(+3*a.col(9))","(+3*a.col(12)+3*a.col(4)%a.col(9)+3*a.col(5)%a.col(8)+3*a.col(8)%a.col(5)+3*a.col(9)%a.col(4))","(+3*a.col(11)+3*a.col(4)%a.col(10)+3*a.col(5)%a.col(1)+3*a.col(8)%a.col(8)+3*a.col(9)%a.col(0))","()","(+1.5*a.col(2))","(+1.5*a.col(3)+1.5*a.col(0)%a.col(2)+3*a.col(1)%a.col(1)+1.5*a.col(2)%a.col(0))","(+1.5*a.col(10))","(+1.5*a.col(11)+1.5*a.col(4)%a.col(10)+3*a.col(8)%a.col(8)+1.5*a.col(10)%a.col(4))","(+1.5*a.col(13)+1.5*a.col(4)%a.col(2)+3*a.col(8)%a.col(1)+1.5*a.col(10)%a.col(0))","()","(+1.5*a.col(2))","(+1.5*a.col(3)+1.5*a.col(0)%a.col(2)+3*a.col(1)%a.col(1)+1.5*a.col(2)%a.col(0))","(+1.5*a.col(10))","(+1.5*a.col(11)+1.5*a.col(4)%a.col(10)+3*a.col(8)%a.col(8)+1.5*a.col(10)%a.col(4))","(+1.5*a.col(13)+1.5*a.col(4)%a.col(2)+3*a.col(8)%a.col(1)+1.5*a.col(10)%a.col(0))","()","()","()","()","()","()"
))

MAT[1,1]='(1+0*a.col(0))'           # ONES(N2)
MAT[5,7]='(1+0*a.col(0))'
MAT[2,13]='(1+0*a.col(0))'
MAT[9,19]='(1+0*a.col(0))'
MAT[9,25]='(1+0*a.col(0))'
MAT[6,31]='(1+0*a.col(0))'

# HOMOGENEITY -----------------------------------------------------------------
func.list.timehomo=func.list*0
for(i in which(func.list==1))
{
# which expressions vary over time
result=eval(body(namess[i]))
func.list.timehomo[i]=2-(sum(diff(result)==0)==(length(result)-1))
}
if(any(func.list.timehomo==2)){homo=F}
#func.list.timehomo[c(1+(0:5)*6)]=1    # Always set these to 1

# BUILD ODE -------------------------------------------------------------------
dims=rep('(',8)
for(i in 1:8)
{
for(j in which(func.list==1))
{
if(MAT[i,j]!="()")
{
dims[i]=paste0(dims[i],'+(',body(namess[j])[2],')',c('*','%')[func.list.timehomo[j]],MAT[i,j])
}
}
dims[i]=paste0(dims[i],')')
}
if(any(dims=='()')){dims[which(dims=='()')]='(1+0*a.col(0))'}
for(i in 1:8)
{
dims[i]=paste0(paste0(paste0('   atemp.col(',i-1,')='),dims[i]),';')
}

# WRIGHT AND SOURCE -----------------------------------------------------------
txt.full=paste(txtA,'\n',dims[1],'\n',dims[2],'\n',dims[3],'\n',dims[4],'\n',dims[5],'\n',dims[6],'\n',dims[7],'\n',dims[8],txtB,txtC)

}

if(state4)
{
# DATA RESOLUTION -------------------------------------------------------------
if(!homo.res)
{
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,
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,diffs,diffs,diffs,diffs,diffs,diffs)
}
}
# REFERENCE MATRIX ------------------------------------------------------------
MAT=rbind(
c("()","(+1*a.col(0))","(+1*a.col(0)%a.col(0)+1*a.col(1))","(+1*a.col(4))","(+1*a.col(4)%a.col(4)+1*a.col(5))","(+1*a.col(8)+1*a.col(0)%a.col(4))","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()"
),c("()","(+2*a.col(1))","(+2*a.col(0)%a.col(1)+2*a.col(1)%a.col(0)+2*a.col(2))","(+2*a.col(8))","(+2*a.col(4)%a.col(8)+2*a.col(8)%a.col(4)+2*a.col(9))","(+2*a.col(10)+2*a.col(0)%a.col(8)+2*a.col(1)%a.col(4))","()","()","()","()","()","()","()","(+1*a.col(0))","(+1*a.col(1)+1*a.col(0)%a.col(0))","(+1*a.col(4))","(+1*a.col(5)+1*a.col(4)%a.col(4))","(+1*a.col(8)+1*a.col(4)%a.col(0))","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()"
),c("()","(+3*a.col(2))","(+3*a.col(0)%a.col(2)+6*a.col(1)%a.col(1)+3*a.col(2)%a.col(0)+3*a.col(3))","(+3*a.col(10))","(+3*a.col(4)%a.col(10)+6*a.col(8)%a.col(8)+3*a.col(10)%a.col(4)+3*a.col(11))","(+3*a.col(13)+3*a.col(0)%a.col(10)+6*a.col(1)%a.col(8)+3*a.col(2)%a.col(4))","()","()","()","()","()","()","()","(+3*a.col(1))","(+3*a.col(2)+3*a.col(0)%a.col(1)+3*a.col(1)%a.col(0))","(+3*a.col(8))","(+3*a.col(9)+3*a.col(4)%a.col(8)+3*a.col(8)%a.col(4))","(+3*a.col(10)+3*a.col(4)%a.col(1)+3*a.col(8)%a.col(0))","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()"
),c("()","(+4*a.col(3))","(+4*a.col(0)%a.col(3)+12*a.col(1)%a.col(2)+12*a.col(2)%a.col(1)+4*a.col(3)%a.col(0))","(+4*a.col(13))","(+4*a.col(4)%a.col(13)+12*a.col(8)%a.col(10)+12*a.col(10)%a.col(8)+4*a.col(13)%a.col(4))","(+4*a.col(0)%a.col(13)+12*a.col(1)%a.col(10)+12*a.col(2)%a.col(8)+4*a.col(3)%a.col(4))","()","()","()","()","()","()","()","(+6*a.col(2))","(+6*a.col(3)+6*a.col(0)%a.col(2)+12*a.col(1)%a.col(1)+6*a.col(2)%a.col(0))","(+6*a.col(10))","(+6*a.col(11)+6*a.col(4)%a.col(10)+12*a.col(8)%a.col(8)+6*a.col(10)%a.col(4))","(+6*a.col(13)+6*a.col(4)%a.col(2)+12*a.col(8)%a.col(1)+6*a.col(10)%a.col(0))","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()"
),c("()","()","()","()","()","()","()","(+1*a.col(0))","(+1*a.col(0)%a.col(0)+1*a.col(1))","(+1*a.col(4))","(+1*a.col(4)%a.col(4)+1*a.col(5))","(+1*a.col(8)+1*a.col(4)%a.col(0))","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()"
),c("()","()","()","()","()","()","()","(+2*a.col(8))","(+2*a.col(0)%a.col(8)+2*a.col(8)%a.col(0)+2*a.col(10))","(+2*a.col(5))","(+2*a.col(4)%a.col(5)+2*a.col(5)%a.col(4)+2*a.col(6))","(+2*a.col(9)+2*a.col(4)%a.col(8)+2*a.col(5)%a.col(0))","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","(+1*a.col(0))","(+1*a.col(1)+1*a.col(0)%a.col(0))","(+1*a.col(4))","(+1*a.col(5)+1*a.col(4)%a.col(4))","(+1*a.col(8)+1*a.col(4)%a.col(0))"
),c("()","()","()","()","()","()","()","(+3*a.col(9))","(+3*a.col(0)%a.col(9)+6*a.col(8)%a.col(8)+3*a.col(9)%a.col(0)+3*a.col(11))","(+3*a.col(6))","(+3*a.col(4)%a.col(6)+6*a.col(5)%a.col(5)+3*a.col(6)%a.col(4)+3*a.col(7))","(+3*a.col(12)+3*a.col(4)%a.col(9)+6*a.col(5)%a.col(8)+3*a.col(6)%a.col(0))","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","(+3*a.col(8))","(+3*a.col(10)+3*a.col(0)%a.col(8)+3*a.col(8)%a.col(0))","(+3*a.col(5))","(+3*a.col(6)+3*a.col(4)%a.col(5)+3*a.col(5)%a.col(4))","(+3*a.col(9)+3*a.col(4)%a.col(8)+3*a.col(5)%a.col(0))"
),c("()","()","()","()","()","()","()","(+4*a.col(12))","(+4*a.col(0)%a.col(12)+12*a.col(8)%a.col(9)+12*a.col(9)%a.col(8)+4*a.col(12)%a.col(0))","(+4*a.col(7))","(+4*a.col(4)%a.col(7)+12*a.col(5)%a.col(6)+12*a.col(6)%a.col(5)+4*a.col(7)%a.col(4))","(+4*a.col(4)%a.col(12)+12*a.col(5)%a.col(9)+12*a.col(6)%a.col(8)+4*a.col(7)%a.col(0))","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","()","(+6*a.col(9))","(+6*a.col(11)+6*a.col(0)%a.col(9)+12*a.col(8)%a.col(8)+6*a.col(9)%a.col(0))","(+6*a.col(6))","(+6*a.col(7)+6*a.col(4)%a.col(6)+12*a.col(5)%a.col(5)+6*a.col(6)%a.col(4))","(+6*a.col(12)+6*a.col(4)%a.col(9)+12*a.col(5)%a.col(8)+6*a.col(6)%a.col(0))"
),c("()","(+1*a.col(8))","(+1*a.col(0)%a.col(8)+1*a.col(8)%a.col(0)+1*a.col(10))","(+1*a.col(5))","(+1*a.col(4)%a.col(5)+1*a.col(5)%a.col(4)+1*a.col(6))","(+1*a.col(9)+1*a.col(0)%a.col(5)+1*a.col(8)%a.col(4))","()","(+1*a.col(1))","(+1*a.col(0)%a.col(1)+1*a.col(1)%a.col(0)+1*a.col(2))","(+1*a.col(8))","(+1*a.col(4)%a.col(8)+1*a.col(8)%a.col(4)+1*a.col(9))","(+1*a.col(10)+1*a.col(4)%a.col(1)+1*a.col(8)%a.col(0))","()","()","()","()","()","()","()","(+0.5*a.col(0))","(+0.5*a.col(1)+0.5*a.col(0)%a.col(0))","(+0.5*a.col(4))","(+0.5*a.col(5)+0.5*a.col(4)%a.col(4))","(+0.5*a.col(8)+0.5*a.col(4)%a.col(0))","()","(+0.5*a.col(0))","(+0.5*a.col(1)+0.5*a.col(0)%a.col(0))","(+0.5*a.col(4))","(+0.5*a.col(5)+0.5*a.col(4)%a.col(4))","(+0.5*a.col(8)+0.5*a.col(4)%a.col(0))","()","()","()","()","()","()"
),c("()","(+1*a.col(9))","(+1*a.col(0)%a.col(9)+2*a.col(8)%a.col(8)+1*a.col(9)%a.col(0)+1*a.col(11))","(+1*a.col(6))","(+1*a.col(4)%a.col(6)+2*a.col(5)%a.col(5)+1*a.col(6)%a.col(4)+1*a.col(7))","(+1*a.col(12)+1*a.col(0)%a.col(6)+2*a.col(8)%a.col(5)+1*a.col(9)%a.col(4))","()","(+2*a.col(10))","(+2*a.col(0)%a.col(10)+2*a.col(8)%a.col(1)+2*a.col(1)%a.col(8)+2*a.col(10)%a.col(0)+2*a.col(13))","(+2*a.col(9))","(+2*a.col(4)%a.col(9)+2*a.col(5)%a.col(8)+2*a.col(8)%a.col(5)+2*a.col(9)%a.col(4)+2*a.col(12))","(+2*a.col(11)+2*a.col(4)%a.col(10)+2*a.col(5)%a.col(1)+2*a.col(8)%a.col(8)+2*a.col(9)%a.col(0))","()","()","()","()","()","()","()","(+1*a.col(8))","(+1*a.col(10)+1*a.col(0)%a.col(8)+1*a.col(8)%a.col(0))","(+1*a.col(5))","(+1*a.col(6)+1*a.col(4)%a.col(5)+1*a.col(5)%a.col(4))","(+1*a.col(9)+1*a.col(4)%a.col(8)+1*a.col(5)%a.col(0))","()","(+1*a.col(8))","(+1*a.col(10)+1*a.col(0)%a.col(8)+1*a.col(8)%a.col(0))","(+1*a.col(5))","(+1*a.col(6)+1*a.col(4)%a.col(5)+1*a.col(5)%a.col(4))","(+1*a.col(9)+1*a.col(4)%a.col(8)+1*a.col(5)%a.col(0))","()","(+1*a.col(1))","(+1*a.col(2)+1*a.col(0)%a.col(1)+1*a.col(1)%a.col(0))","(+1*a.col(8))","(+1*a.col(9)+1*a.col(4)%a.col(8)+1*a.col(8)%a.col(4))","(+1*a.col(10)+1*a.col(4)%a.col(1)+1*a.col(8)%a.col(0))"
),c("()","(+2*a.col(10))","(+2*a.col(0)%a.col(10)+2*a.col(8)%a.col(1)+2*a.col(1)%a.col(8)+2*a.col(10)%a.col(0)+2*a.col(13))","(+2*a.col(9))","(+2*a.col(4)%a.col(9)+2*a.col(5)%a.col(8)+2*a.col(8)%a.col(5)+2*a.col(9)%a.col(4)+2*a.col(12))","(+2*a.col(11)+2*a.col(0)%a.col(9)+2*a.col(8)%a.col(8)+2*a.col(1)%a.col(5)+2*a.col(10)%a.col(4))","()","(+1*a.col(2))","(+1*a.col(0)%a.col(2)+2*a.col(1)%a.col(1)+1*a.col(2)%a.col(0)+1*a.col(3))","(+1*a.col(10))","(+1*a.col(4)%a.col(10)+2*a.col(8)%a.col(8)+1*a.col(10)%a.col(4)+1*a.col(11))","(+1*a.col(13)+1*a.col(4)%a.col(2)+2*a.col(8)%a.col(1)+1*a.col(10)%a.col(0))","()","(+1*a.col(8))","(+1*a.col(10)+1*a.col(0)%a.col(8)+1*a.col(8)%a.col(0))","(+1*a.col(5))","(+1*a.col(6)+1*a.col(4)%a.col(5)+1*a.col(5)%a.col(4))","(+1*a.col(9)+1*a.col(4)%a.col(8)+1*a.col(5)%a.col(0))","()","(+1*a.col(1))","(+1*a.col(2)+1*a.col(0)%a.col(1)+1*a.col(1)%a.col(0))","(+1*a.col(8))","(+1*a.col(9)+1*a.col(4)%a.col(8)+1*a.col(8)%a.col(4))","(+1*a.col(10)+1*a.col(4)%a.col(1)+1*a.col(8)%a.col(0))","()","(+1*a.col(1))","(+1*a.col(2)+1*a.col(0)%a.col(1)+1*a.col(1)%a.col(0))","(+1*a.col(8))","(+1*a.col(9)+1*a.col(4)%a.col(8)+1*a.col(8)%a.col(4))","(+1*a.col(10)+1*a.col(4)%a.col(1)+1*a.col(8)%a.col(0))","()","()","()","()","()","()"
),c("()","(+2*a.col(11))","(+2*a.col(0)%a.col(11)+4*a.col(8)%a.col(10)+2*a.col(9)%a.col(1)+2*a.col(1)%a.col(9)+4*a.col(10)%a.col(8)+2*a.col(11)%a.col(0))","(+2*a.col(12))","(+2*a.col(4)%a.col(12)+4*a.col(5)%a.col(9)+2*a.col(6)%a.col(8)+2*a.col(8)%a.col(6)+4*a.col(9)%a.col(5)+2*a.col(12)%a.col(4))","(+2*a.col(0)%a.col(12)+4*a.col(8)%a.col(9)+2*a.col(9)%a.col(8)+2*a.col(1)%a.col(6)+4*a.col(10)%a.col(5)+2*a.col(11)%a.col(4))","()","(+2*a.col(13))","(+2*a.col(0)%a.col(13)+2*a.col(8)%a.col(2)+4*a.col(1)%a.col(10)+4*a.col(10)%a.col(1)+2*a.col(2)%a.col(8)+2*a.col(13)%a.col(0))","(+2*a.col(11))","(+2*a.col(4)%a.col(11)+2*a.col(5)%a.col(10)+4*a.col(8)%a.col(9)+4*a.col(9)%a.col(8)+2*a.col(10)%a.col(5)+2*a.col(11)%a.col(4))","(+2*a.col(4)%a.col(13)+2*a.col(5)%a.col(2)+4*a.col(8)%a.col(10)+4*a.col(9)%a.col(1)+2*a.col(10)%a.col(8)+2*a.col(11)%a.col(0))","()","(+1*a.col(9))","(+1*a.col(11)+1*a.col(0)%a.col(9)+2*a.col(8)%a.col(8)+1*a.col(9)%a.col(0))","(+1*a.col(6))","(+1*a.col(7)+1*a.col(4)%a.col(6)+2*a.col(5)%a.col(5)+1*a.col(6)%a.col(4))","(+1*a.col(12)+1*a.col(4)%a.col(9)+2*a.col(5)%a.col(8)+1*a.col(6)%a.col(0))","()","(+2*a.col(10))","(+2*a.col(13)+2*a.col(0)%a.col(10)+2*a.col(8)%a.col(1)+2*a.col(1)%a.col(8)+2*a.col(10)%a.col(0))","(+2*a.col(9))","(+2*a.col(12)+2*a.col(4)%a.col(9)+2*a.col(5)%a.col(8)+2*a.col(8)%a.col(5)+2*a.col(9)%a.col(4))","(+2*a.col(11)+2*a.col(4)%a.col(10)+2*a.col(5)%a.col(1)+2*a.col(8)%a.col(8)+2*a.col(9)%a.col(0))","()","(+2*a.col(10))","(+2*a.col(13)+2*a.col(0)%a.col(10)+2*a.col(8)%a.col(1)+2*a.col(1)%a.col(8)+2*a.col(10)%a.col(0))","(+2*a.col(9))","(+2*a.col(12)+2*a.col(4)%a.col(9)+2*a.col(5)%a.col(8)+2*a.col(8)%a.col(5)+2*a.col(9)%a.col(4))","(+2*a.col(11)+2*a.col(4)%a.col(10)+2*a.col(5)%a.col(1)+2*a.col(8)%a.col(8)+2*a.col(9)%a.col(0))","()","(+1*a.col(2))","(+1*a.col(3)+1*a.col(0)%a.col(2)+2*a.col(1)%a.col(1)+1*a.col(2)%a.col(0))","(+1*a.col(10))","(+1*a.col(11)+1*a.col(4)%a.col(10)+2*a.col(8)%a.col(8)+1*a.col(10)%a.col(4))","(+1*a.col(13)+1*a.col(4)%a.col(2)+2*a.col(8)%a.col(1)+1*a.col(10)%a.col(0))"
),c("()","(+1*a.col(12))","(+1*a.col(0)%a.col(12)+3*a.col(8)%a.col(9)+3*a.col(9)%a.col(8)+1*a.col(12)%a.col(0))","(+1*a.col(7))","(+1*a.col(4)%a.col(7)+3*a.col(5)%a.col(6)+3*a.col(6)%a.col(5)+1*a.col(7)%a.col(4))","(+1*a.col(0)%a.col(7)+3*a.col(8)%a.col(6)+3*a.col(9)%a.col(5)+1*a.col(12)%a.col(4))","()","(+3*a.col(11))","(+3*a.col(0)%a.col(11)+6*a.col(8)%a.col(10)+3*a.col(9)%a.col(1)+3*a.col(1)%a.col(9)+6*a.col(10)%a.col(8)+3*a.col(11)%a.col(0))","(+3*a.col(12))","(+3*a.col(4)%a.col(12)+6*a.col(5)%a.col(9)+3*a.col(6)%a.col(8)+3*a.col(8)%a.col(6)+6*a.col(9)%a.col(5)+3*a.col(12)%a.col(4))","(+3*a.col(4)%a.col(11)+6*a.col(5)%a.col(10)+3*a.col(6)%a.col(1)+3*a.col(8)%a.col(9)+6*a.col(9)%a.col(8)+3*a.col(12)%a.col(0))","()","()","()","()","()","()","()","(+1.5*a.col(9))","(+1.5*a.col(11)+1.5*a.col(0)%a.col(9)+3*a.col(8)%a.col(8)+1.5*a.col(9)%a.col(0))","(+1.5*a.col(6))","(+1.5*a.col(7)+1.5*a.col(4)%a.col(6)+3*a.col(5)%a.col(5)+1.5*a.col(6)%a.col(4))","(+1.5*a.col(12)+1.5*a.col(4)%a.col(9)+3*a.col(5)%a.col(8)+1.5*a.col(6)%a.col(0))","()","(+1.5*a.col(9))","(+1.5*a.col(11)+1.5*a.col(0)%a.col(9)+3*a.col(8)%a.col(8)+1.5*a.col(9)%a.col(0))","(+1.5*a.col(6))","(+1.5*a.col(7)+1.5*a.col(4)%a.col(6)+3*a.col(5)%a.col(5)+1.5*a.col(6)%a.col(4))","(+1.5*a.col(12)+1.5*a.col(4)%a.col(9)+3*a.col(5)%a.col(8)+1.5*a.col(6)%a.col(0))","()","(+3*a.col(10))","(+3*a.col(13)+3*a.col(0)%a.col(10)+3*a.col(8)%a.col(1)+3*a.col(1)%a.col(8)+3*a.col(10)%a.col(0))","(+3*a.col(9))","(+3*a.col(12)+3*a.col(4)%a.col(9)+3*a.col(5)%a.col(8)+3*a.col(8)%a.col(5)+3*a.col(9)%a.col(4))","(+3*a.col(11)+3*a.col(4)%a.col(10)+3*a.col(5)%a.col(1)+3*a.col(8)%a.col(8)+3*a.col(9)%a.col(0))"
),c("()","(+3*a.col(13))","(+3*a.col(0)%a.col(13)+3*a.col(8)%a.col(2)+6*a.col(1)%a.col(10)+6*a.col(10)%a.col(1)+3*a.col(2)%a.col(8)+3*a.col(13)%a.col(0))","(+3*a.col(11))","(+3*a.col(4)%a.col(11)+3*a.col(5)%a.col(10)+6*a.col(8)%a.col(9)+6*a.col(9)%a.col(8)+3*a.col(10)%a.col(5)+3*a.col(11)%a.col(4))","(+3*a.col(0)%a.col(11)+3*a.col(8)%a.col(10)+6*a.col(1)%a.col(9)+6*a.col(10)%a.col(8)+3*a.col(2)%a.col(5)+3*a.col(13)%a.col(4))","()","(+1*a.col(3))","(+1*a.col(0)%a.col(3)+3*a.col(1)%a.col(2)+3*a.col(2)%a.col(1)+1*a.col(3)%a.col(0))","(+1*a.col(13))","(+1*a.col(4)%a.col(13)+3*a.col(8)%a.col(10)+3*a.col(10)%a.col(8)+1*a.col(13)%a.col(4))","(+1*a.col(4)%a.col(3)+3*a.col(8)%a.col(2)+3*a.col(10)%a.col(1)+1*a.col(13)%a.col(0))","()","(+3*a.col(10))","(+3*a.col(13)+3*a.col(0)%a.col(10)+3*a.col(8)%a.col(1)+3*a.col(1)%a.col(8)+3*a.col(10)%a.col(0))","(+3*a.col(9))","(+3*a.col(12)+3*a.col(4)%a.col(9)+3*a.col(5)%a.col(8)+3*a.col(8)%a.col(5)+3*a.col(9)%a.col(4))","(+3*a.col(11)+3*a.col(4)%a.col(10)+3*a.col(5)%a.col(1)+3*a.col(8)%a.col(8)+3*a.col(9)%a.col(0))","()","(+1.5*a.col(2))","(+1.5*a.col(3)+1.5*a.col(0)%a.col(2)+3*a.col(1)%a.col(1)+1.5*a.col(2)%a.col(0))","(+1.5*a.col(10))","(+1.5*a.col(11)+1.5*a.col(4)%a.col(10)+3*a.col(8)%a.col(8)+1.5*a.col(10)%a.col(4))","(+1.5*a.col(13)+1.5*a.col(4)%a.col(2)+3*a.col(8)%a.col(1)+1.5*a.col(10)%a.col(0))","()","(+1.5*a.col(2))","(+1.5*a.col(3)+1.5*a.col(0)%a.col(2)+3*a.col(1)%a.col(1)+1.5*a.col(2)%a.col(0))","(+1.5*a.col(10))","(+1.5*a.col(11)+1.5*a.col(4)%a.col(10)+3*a.col(8)%a.col(8)+1.5*a.col(10)%a.col(4))","(+1.5*a.col(13)+1.5*a.col(4)%a.col(2)+3*a.col(8)%a.col(1)+1.5*a.col(10)%a.col(0))","()","()","()","()","()","()"
))

MAT[1,1]='(1+0*a.col(0))'           # ONES(N2)
MAT[5,7]='(1+0*a.col(0))'
MAT[2,13]='(1+0*a.col(0))'
MAT[9,19]='(1+0*a.col(0))'
MAT[9,25]='(1+0*a.col(0))'
MAT[6,31]='(1+0*a.col(0))'

# HOMOGENEITY -----------------------------------------------------------------

func.list.timehomo=func.list*0
for(i in which(func.list==1))
{
# which expressions vary over time
result=eval(body(namess[i]))
func.list.timehomo[i]=2-(sum(diff(result)==0)==(length(result)-1))
}
if(any(func.list.timehomo==2)){homo=F}
#func.list.timehomo[c(1+(0:5)*6)]=1    # Always set these to 1

# BUILD ODE -------------------------------------------------------------------
dims=rep('(',14)
for(i in 1:14)
{
for(j in which(func.list==1))
{
if(MAT[i,j]!="()")
{
dims[i]=paste0(dims[i],'+(',body(namess[j])[2],')',c('*','%')[func.list.timehomo[j]],MAT[i,j])
}
}
dims[i]=paste0(dims[i],')')
}
if(any(dims=='()')){dims[which(dims=='()')]='(1+0*a.col(0))'}
for(i in 1:14)
{
dims[i]=paste0(paste0(paste0('        atemp.col(',i-1,')='),dims[i]),';')
}

# WRIGHT AND SOURCE -----------------------------------------------------------
txt.full=paste(txtA,'\n',dims[1],'\n',dims[2],'\n',dims[3],'\n',dims[4],'\n',dims[5],'\n',dims[6],'\n',dims[7],'\n',dims[8],'\n',dims[9],'\n',dims[10],'\n',dims[11],'\n',dims[12],'\n',dims[13],'\n',dims[14],txtB,txtC)
}

#library(Rcpp)
if(wrt)
{
write(txt.full,'BiGQD.mle.cpp')
}
cat(stre, " \r")
flush.console()
sourceCpp(code=txt.full)
cat('                                     ','\r')

#==============================================================================
#                           Interface Module
#==============================================================================
trim <- function (x) gsub("([[:space:]])", "", x)

namess4=c('a00','a10','a20','a01','a02','a11',
'b00','b10','b20','b01','b02','b11',
'c00','c10','c20','c01','c02','c11',
'd00','d10','d20','d01','d02','d11',
'e00','e10','e20','e01','e02','e11',
'f00','f10','f20','f01','f02','f11')
indnames =rep(0,36)
for(i in 1:36)
{
if(sum(obs==namess4[i]))
{
indnames[i]=TRUE
namess4[i]=paste0(namess4[i],' : ',trim(body(namess4[i])[2]))
}
}
namess4=matrix(namess4,length(namess4),1)
buffer0=c('================================================================')
buffer1=c('----------------------------------------------------------------')
buffer2=c('................................................................')
buffer3=c('...   ...   ...   ...   ...   ...   ...   ...   ...   ...   ... ')
buffer4=c('_____________________ Drift Coefficients _______________________')
buffer5=c('___________________ Diffusion Coefficients _____________________')
buffer6=c('_____________________ Prior Distributions ______________________')
buffer7=c('__________________________ Model Info __________________________')

#Info1=data.frame(rbind(buffer0,matrix(names4[1:6,],6,1),buffer3,matrix(names4[1:6+6,],6,1),buffer3,matrix(names4[1:6+12,],6,1),buffer3,matrix(names4[1:6+18,],6,1),buffer3,matrix(names4[1:6+24,],6,1),buffer3,matrix(names4[1:6+30,],6,1),buffer2,prior.list,buffer),check.names=F)
#colnames(Info1)=type.sol
#buffer00=data.frame(buffer0,check.names=F)
#colnames(buffer00)=''
#print(buffer00 ,row.names = FALSE,right=F)
#print(Info1,row.names = FALSE,right=F)
# print(indnames)
Info=c(buffer0,type.sol,buffer0,buffer4,
namess4[1:6][which(indnames[1:6]==T)],
buffer3,
namess4[7:12][which(indnames[7:12]==T)],
buffer5,
namess4[13:18][which(indnames[13:18]==T)],
buffer3,
namess4[19:24][which(indnames[19:24]==T)],
buffer3,
namess4[25:30][which(indnames[25:30]==T)],
buffer3,
namess4[31:36][which(indnames[31:36]==T)])
Info=data.frame(matrix(Info,length(Info),1))
colnames(Info)=''
if(print.output)
{
print(Info,row.names = FALSE,right=F)
}
############################################################################
############################################################################
############################################################################

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")
}
if(is.null(method))
{
}
if(is.null(control))
{
control=list(fnscale=-1)
}
if(!is.null(control))
{
if(any(names(control)=='fnscale'))
{
if(control\$fnscale>=0)
{
stop(' ==============================================================================
Incorrect input: control\$fnscale must be negative!
==============================================================================')
}
}
if(!any(names(control)=='fnscale'))
{
control=c(control,fnscale=-1)
}
if(!any(names(control)=='maxit'))
{
control=c(control,maxit=5000)
}
}
strts =matrix(0,nnn-1,2)
strts[,1]=rtf[1]
strts[,2]=rtf[2]
like=function(pars)
{
rs=solver(X1[-nnn],X2[-nnn],X1[-1],X2[-1],c(0,pars),mesh,delt,nnn-1,T.seq[-nnn],strts,tro,secmom)
sum(rs\$like[-excl,1])
}

tme=Sys.time()
result=optim(theta,like,control = control,method=method, hessian = T)
tme=tme.eval(tme)

actual.p=length(theta)

model.inf=list(elapsed.time=tme,time.homogeneous=c('Yes','No')[2-homo],p=actual.p,N=length(X[,1])-length(excl)+1,Tag=Tag)
Info2=c( buffer7,
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)),
buffer1)
Info2=data.frame(matrix(Info2,length(Info2),1))
colnames(Info2)=''
if(print.output)
{
print(Info2,row.names = FALSE,right=F)
}
ret=list(opt=result,elapsed.time=tme,model.info=model.inf)
class(ret)='GQD.mle'
return(ret)
}
```

## Try the DiffusionRgqd package in your browser

Any scripts or data that you put into this service are public.

DiffusionRgqd documentation built on May 2, 2019, 3:26 a.m.