Nothing
R/GQD.TIpassage.R
In DiffusionRgqd: Inference and Analysis for Generalized Quadratic Diffusions
GQD.TIpassage=
function(Xs,B,s,t,delt,theta=c(0),IEQ.type='Buonocore',wrt=FALSE)
{
fpt =function(res,res2,N,delt){}
ieq =function(res1,res2,N,delt){}
solver =function(Xs, Xt, theta, N , delt , N2, tt , P , alpha, lower , upper, tro ){}
rm(list =c('fpt','ieq','solver'))
Tmax=t
# Warning Module
if(delt>=0.1){warning("Large delt may result in poor approximations.");}
if(delt>=1) {stop("delt must be <1.");}
if(Xs>B) {stop("Xs must be < B. If the fpt-problem is a down-cross problem, consider the transform Yt = -Xt.");}
pow=function(x,p)
{
x^p
}
prod=function(a,b){a*b}
if(IEQ.type=='Volterra')
{
#stop("Volterra solution offline.");
#===============================================================================
# State 1
#===============================================================================
txtA.state1=
'
#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")]]
mat f(mat a,vec theta,vec t,int N2)
{
mat atemp(N2,2);'
txtB.state1.homo= '
return atemp;
}
// [[Rcpp::export]]
mat solver(vec Xs,vec Xt,vec theta,int N,double delt,int N2,vec tt)
{
mat fx0(N2,2);
mat fx1(N2,2);
mat fx2(N2,2);
mat fx3(N2,2);
mat fx4(N2,2);
mat fx5(N2,2);
mat fx6(N2,2);
mat fx7(N2,2);
mat fx8(N2,2);
mat fx9(N2,2);
mat fx10(N2,2);
mat fx11(N2,2);
mat fx12(N2,2);
mat fx13(N2,2);
mat fx14(N2,2);
mat fx15(N2,2);
mat fx16(N2,2);
mat x0(N2,2);
mat x1(N2,2);
mat x2(N2,2);
mat x3(N2,2);
mat x4(N2,2);
mat x5(N2,2);
mat x6(N2,2);
mat x7(N2,2);
mat x8(N2,2);
mat x9(N2,2);
mat x10(N2,2);
mat x11(N2,2);
mat x12(N2,2);
mat x13(N2,2);
mat x14(N2,2);
mat x15(N2,2);
mat x16(N2,2);
x0.fill(0);
x0.col(0)=Xs;
mat res(N2,N);
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;
res.col(i) =exp(-0.5*log(2*3.141592653589793*x0.col(1))-0.5*((Xt-x0.col(0))%(Xt-x0.col(0))/x0.col(1)));
}
return(res);
}'
#===============================================================================
# State 2
#===============================================================================
txtA.state2=
'
#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")]]
mat f(mat a,vec theta,vec t,int N2)
{
mat atemp(N2,4);'
txtB.state2.homo= '
return atemp;
}
// [[Rcpp::export]]
mat solver(vec Xs,vec Xt,vec theta,int N,double delt,int N2,vec tt)
{
mat fx0(N2,4);
mat fx1(N2,4);
mat fx2(N2,4);
mat fx3(N2,4);
mat fx4(N2,4);
mat fx5(N2,4);
mat fx6(N2,4);
mat fx7(N2,4);
mat fx8(N2,4);
mat fx9(N2,4);
mat fx10(N2,4);
mat fx11(N2,4);
mat fx12(N2,4);
mat fx13(N2,4);
mat fx14(N2,4);
mat fx15(N2,4);
mat fx16(N2,4);
mat x0(N2,4);
mat x1(N2,4);
mat x2(N2,4);
mat x3(N2,4);
mat x4(N2,4);
mat x5(N2,4);
mat x6(N2,4);
mat x7(N2,4);
mat x8(N2,4);
mat x9(N2,4);
mat x10(N2,4);
mat x11(N2,4);
mat x12(N2,4);
mat x13(N2,4);
mat x14(N2,4);
mat x15(N2,4);
mat x16(N2,4);
x0.fill(0);
x0.col(0)=Xs;
mat res(N2,N);
vec d=tt;
vec p(N2);
vec q(N2);
vec chk(N2);
vec th(N2);
vec K(N2);
vec K1(N2);
vec K2(N2);
vec val(N2);
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;
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);
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);
chk=pow(q,2)/4.0 + pow(p,3)/27.0;
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));
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;
K1 =x0.col(0) +(x0.col(1)%th) +(x0.col(2)%th%th)/2.0 +(x0.col(3)%th%th%th)/6.0;
K2 =x0.col(1) +(x0.col(2)%th) +(x0.col(3)%th%th)/2.0;
val=-0.5*log(2*3.141592653589793*K2)+(K-th%K1);
res.col(i)=exp(val);
}
return(res);
}'
code2=
'
// [[Rcpp::export]]
vec fpt(mat res1,vec res2,int N,double delt)
{
vec g(N);
g.fill(0);
double sm;
for (int i = 1; i < N; i++)
{
sm=0;
for (int j = 1; j < i; j++)
{
sm=sm+res1(j,i-j)*g(j);
}
g(i) = (res2(i)/delt -sm)/res1(i,1);
}
return g;
}
'
#==============================================================================
# Data resolution Module
#==============================================================================
Tmax=t
t=seq(0,100,by=1/100)
homo=T
homo.res=T
#==============================================================================
# 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)
state2=!state1
sol.state=c(state1,state2)
sol.state=c('2nd Ord. Truncation Normal Dist.','4th Ord. Truncation Saddlepoint Appr.')[which(sol.state==T)]
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)/N.mesh,diff(T.seq)/N.mesh)
#
#}
# REFERENCE MATRIX ------------------------------------------------------------
MAT=rbind(
c('1','a.col(0)' ,''),
c('','2*a.col(1)','1'),
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]))
#print(diff(result)[1])
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,3)]=1 # Always set these to 1
#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'}
for(i in 1:2)
{
dims[i]=paste0(paste0(paste0(' atemp.col(',i-1,')='),dims[i]),';')
}
# WRIGHT AND SOURCE -----------------------------------------------------------
if(homo.res)
{
txt.full=paste(txtA.state1,'\n',dims[1],'\n',dims[2],txtB.state1.homo,code2)
}
#if(!homo.res)
#{
# txt.full=paste(txtA.state1,'\n',dims[1],'\n',dims[2],txtB.state1.NONhomo,txtC.state1)
#}
#write(txt.full,'FPT_modelbuild.cpp')
type.sol =" Gneralized Ornstein-Uhlenbeck "#paste0("GENERALIZED LINEAR DIFFUSON",c('',": TIME INHOMOGNEOUS")[2-homo],c(')','+INDEPENDENT')[1+state1])
}
if(state2)
{
# DATA RESOLUTION -------------------------------------------------------------
#if(!homo.res)
#{
# delt=cbind(diff(T.seq)/N.mesh,diff(T.seq)/N.mesh,diff(T.seq)/N.mesh,diff(T.seq)/N.mesh)
#
#}
# REFERENCE MATRIX ------------------------------------------------------------
MAT=rbind(
c('1','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','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))'))
# 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]))
#print(diff(result)[1])
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,4)]=1 # Always set these to 1
#BUILD ODE --------------------------------------------------------------------
dims=rep('(',4)
for(i in 1:4)
{
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'}
for(i in 1:4)
{
dims[i]=paste0(paste0(paste0(' atemp.col(',i-1,')='),dims[i]),';')
}
# WRIGHT AND SOURCE -----------------------------------------------------------
if(homo.res)
{
txt.full=paste(txtA.state2,'\n',dims[1],'\n',dims[2],'\n',dims[3],'\n',dims[4],txtB.state2.homo,code2)
}
#if(!homo.res)
#{
# txt.full=paste(txtA.state2,'\n',dims[1],'\n',dims[2],'\n',dims[3],'\n',dims[4],txtB.state2.NONhomo,txtC.state2)
#}
#write(txt.full,'FPT_modelbuild.cpp')
type.sol =" Generalized Quadratic Diffusion (GQD) "#paste0("GENERALIZED LINEAR DIFFUSON",c('',": TIME INHOMOGNEOUS")[2-homo],c(')','+INDEPENDENT')[1+state1])
}
#library(Rcpp)
#library(RcppArmadillo)
if(wrt)
{
write(txt.full,'FPT_modelbuild.cpp')
}
stre="Compiling C++ code. Please wait."
cat(stre, " \r")
sourceCpp(code=txt.full)
cat(' ','\r')
N=(Tmax-s)/delt
tt=seq(s,Tmax-delt,delt)
res=solver(rep(B,N),rep(B,N),c(0,theta),N,delt,N,tt)
res2=solver(rep(Xs,1),rep(B,1),c(0,theta),N,delt,1,0)
y=fpt(res,res2,N,delt)
#return(list(density=y,time=tt))
#print('Done.')
}
################################################################################
if(IEQ.type=='Buonocore')
{
#===============================================================================
# State 1
#===============================================================================
#_______________________________________________________________________________
txt1.1A= '
#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);'
#_______________________________________________________________________________
txt1.2A= '
return atemp;
}
'
#_______________________________________________________________________________
txt2A='
// [[Rcpp::export]]
vec saddle(mat xx,double Bt)
{
vec val=exp(-0.5*log(2*3.141592653589793*xx.col(1))-0.5*(Bt-xx.col(0))%(Bt-xx.col(0))/xx.col(1));
return(val);
}
// [[Rcpp::export]]
vec saddle2(mat xx,double Bt)
{
vec val=-exp(-0.5*log(2*3.141592653589793*xx.col(1))-0.5*(Bt-xx.col(0))%(Bt-xx.col(0))/xx.col(1))%(Bt-xx.col(0))/xx.col(1);
return(val);
}
// [[Rcpp::export]]
vec pcurr(mat xx,double Bt,vec theta,vec t)
{
vec val ='
#_______________________________________________________________________________
txt3A=';
return(val);
}
// [[Rcpp::export]]
mat solver(vec Xs,double Bt,vec theta,int N,double delt,int N2,vec tt)
{
mat fx0(N2,2);
mat fx1(N2,2);
mat fx2(N2,2);
mat fx3(N2,2);
mat fx4(N2,2);
mat fx5(N2,2);
mat fx6(N2,2);
mat fx7(N2,2);
mat fx8(N2,2);
mat fx9(N2,2);
mat fx10(N2,2);
mat fx11(N2,2);
mat fx12(N2,2);
mat fx13(N2,2);
mat fx14(N2,2);
mat fx15(N2,2);
mat fx16(N2,2);
mat x0(N2,2);
mat x1(N2,2);
mat x2(N2,2);
mat x3(N2,2);
mat x4(N2,2);
mat x5(N2,2);
mat x6(N2,2);
mat x7(N2,2);
mat x8(N2,2);
mat x9(N2,2);
mat x10(N2,2);
mat x11(N2,2);
mat x12(N2,2);
mat x13(N2,2);
mat x14(N2,2);
mat x15(N2,2);
mat x16(N2,2);
x0.fill(0);
x0.col(0)=Xs;
mat res(N2,N);
res.fill(0);
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;
res.col(i) = pcurr(x0,Bt,theta,d);
}
return(res);
}
// [[Rcpp::export]]
vec ieq(mat res1,vec res2,int N,double delt)
{
vec v(N);
v(0)=0;
double smtemp=0;
for (int i = 1; i < N; i++)
{
smtemp=0;
for (int j = 1; j < i; j++)
{
smtemp=smtemp+v(j)*res1(j,i-j);
}
v(i) = (2.0*res2(i)-2.0*smtemp*delt);
}
return(v);
}
'
#===============================================================================
# State 2
#===============================================================================
#_______________________________________________________________________________
txt1.1= '
#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);'
#_______________________________________________________________________________
txt1.2= '
return atemp;
}
'
#_______________________________________________________________________________
txt2='
// [[Rcpp::export]]
vec saddle(mat xx,double Bt)
{
vec p=(1.0/3.0) *(3*(xx.col(3)/6.0)%xx.col(1) - pow(xx.col(2)/2.0,2))/pow(xx.col(3)/6.0,2);
vec q=(1.0/27.0)*(27*pow(xx.col(3)/6.0,2)%(xx.col(0)-Bt) - 9*(xx.col(3)/6.0)%(xx.col(2)/2.0)%xx.col(1) + 2*pow(xx.col(2)/2.0,3))/pow(xx.col(3)/6.0,3);
vec chk=pow(q,2)/4.0 + pow(p,3)/27.0;
vec th=-(xx.col(2)/2.0)/(3*(xx.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 =xx.col(0)%th+(xx.col(1)%th%th)/2.0+(xx.col(2)%th%th%th)/6.0 +(xx.col(3)%th%th%th%th)/24.0;
vec K1=xx.col(0) +(xx.col(1)%th) +(xx.col(2)%th%th)/2.0 +(xx.col(3)%th%th%th)/6.0;
vec K2=xx.col(1) +(xx.col(2)%th) +(xx.col(3)%th%th)/2.0;
vec val=exp(-0.5*log(2*3.141592653589793*K2)+(K-th%K1));
return(val);
}
// [[Rcpp::export]]
vec saddle2(mat xx,double Bt)
{
vec p=(1.0/3.0) *(3*(xx.col(3)/6.0)%xx.col(1) - pow(xx.col(2)/2.0,2))/pow(xx.col(3)/6.0,2);
vec q=(1.0/27.0)*(27*pow(xx.col(3)/6.0,2)%(xx.col(0)-Bt) - 9*(xx.col(3)/6.0)%(xx.col(2)/2.0)%xx.col(1) + 2*pow(xx.col(2)/2.0,3))/pow(xx.col(3)/6.0,3);
vec chk=pow(q,2)/4.0 + pow(p,3)/27.0;
vec th=-(xx.col(2)/2.0)/(3*(xx.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 thdash=1.0/(xx.col(1)+th%xx.col(2)+0.5*th%th%xx.col(3));
vec K =xx.col(0)%th+(xx.col(1)%th%th)/2.0+(xx.col(2)%th%th%th)/6.0 +(xx.col(3)%th%th%th%th)/24.0;
vec K1=xx.col(0) +(xx.col(1)%th) +(xx.col(2)%th%th)/2.0 +(xx.col(3)%th%th%th)/6.0;
vec K2=xx.col(1) +(xx.col(2)%th) +(xx.col(3)%th%th)/2.0;
vec gg = 1.0/sqrt(2*3.141592653589793*(K2));
vec ggdash = -3.141592653589793*pow(2*3.141592653589793*(K2),-3.0/2.0)%(xx.col(2)%thdash+xx.col(3)%thdash%th);
vec hh = exp(K-th%K1);
vec hhdash = exp(K-th%K1)%(-th%thdash%xx.col(1)-th%th%thdash%xx.col(2)-0.5*(th%th%th)%thdash%xx.col(3));
vec val=(ggdash%hh+gg%hhdash);
return(val);
}
// [[Rcpp::export]]
vec pcurr(mat xx,double Bt,vec theta,vec t)
{
vec val ='
#_______________________________________________________________________________
txt3=';
return(val);
}
// [[Rcpp::export]]
mat solver(vec Xs,double Bt,vec theta,int N,double delt,int N2,vec tt)
{
mat fx0(N2,4);
mat fx1(N2,4);
mat fx2(N2,4);
mat fx3(N2,4);
mat fx4(N2,4);
mat fx5(N2,4);
mat fx6(N2,4);
mat fx7(N2,4);
mat fx8(N2,4);
mat fx9(N2,4);
mat fx10(N2,4);
mat fx11(N2,4);
mat fx12(N2,4);
mat fx13(N2,4);
mat fx14(N2,4);
mat fx15(N2,4);
mat fx16(N2,4);
mat x0(N2,4);
mat x1(N2,4);
mat x2(N2,4);
mat x3(N2,4);
mat x4(N2,4);
mat x5(N2,4);
mat x6(N2,4);
mat x7(N2,4);
mat x8(N2,4);
mat x9(N2,4);
mat x10(N2,4);
mat x11(N2,4);
mat x12(N2,4);
mat x13(N2,4);
mat x14(N2,4);
mat x15(N2,4);
mat x16(N2,4);
x0.fill(0);
x0.col(0)=Xs;
mat res(N2,N);
res.fill(0);
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;
res.col(i) = pcurr(x0,Bt,theta,d);
}
return(res);
}
// [[Rcpp::export]]
vec ieq(mat res1,vec res2,int N,double delt)
{
vec v(N);
v(0)=0;
double smtemp=0;
for (int i = 1; i < N; i++)
{
smtemp=0;
for (int j = 1; j < i; j++)
{
smtemp=smtemp+v(j)*res1(j,i-j);
}
v(i) = (2.0*res2(i)-2.0*smtemp*delt);
}
return(v);
}
'
# BUILD PROB. CURRENT ---------------------------------------------------------
t=seq(0,100,by=1/100)
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}
}
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))
}
T1=''
T2=''
T3=''
if(func.list[1]==1)
{
T1=paste0(T1,'+(',body(namess[1])[2],')')
}
if(func.list[2]==1)
{
T1=paste0(T1,'+(',body(namess[2])[2],')*Bt')
}
if(func.list[3]==1)
{
T1=paste0(T1,'+(',body(namess[3])[2],')*Bt*Bt')
}
if(func.list[4]==1)
{
T3=paste0(T3,'+(',body(namess[4])[2],')')
}
if(func.list[5]==1)
{
T2=paste0(T2,'+(',body(namess[5])[2],')')
T3=paste0(T3,'+(',body(namess[5])[2],')*Bt')
}
if(func.list[6]==1)
{
T2=paste0(T2,'+2*(',body(namess[6])[2],')*Bt')
T3=paste0(T3,'+(',body(namess[6])[2],')*Bt*Bt')
}
T1=paste0('(',T1,')')
T2=paste0('(',T2,')')
T3=paste0('(',T3,')')
T1.inhom=F
T2.inhom=F
T3.inhom=F
if(any(func.list.timehomo[c(1,2,3)]==2)){T1.inhom=T}
if(any(func.list.timehomo[c(5,6)]==2)) {T2.inhom=T}
if(any(func.list.timehomo[c(4,5,6)]==2)){T3.inhom=T}
text=paste0(txt2,'0.5*',T1,c('*','%')[T1.inhom+1],'saddle(xx,Bt)-','0.75*',T2,c('*','%')[T2.inhom+1],'saddle(xx,Bt)-','0.5*',T3,c('*','%')[T3.inhom+1],'saddle2(xx,Bt)',txt3)
func.list.timehomo[c(1,4)]=1 # Always set these to 1
# REFERENCE MATRIX ------------------------------------------------------------
MAT=rbind(
c('1','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','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))'))
#BUILD ODE --------------------------------------------------------------------
dims=rep('(',4)
for(i in 1:4)
{
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=='()')]='0'}
for(i in 1:4)
{
dims[i]=paste0(paste0(paste0(' atemp.col(',i-1,')='),dims[i]),';')
}
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])
}
# WRIGHT AND SOURCE -----------------------------------------------------------
if(sum(func.list[c(3,5,6)]==0)==3)
{
text=paste0(txt2A,'0.5*',T1,c('*','%')[T1.inhom+1],'saddle(xx,Bt)-','0.5*',T3,c('*','%')[T3.inhom+1],'saddle2(xx,Bt)',txt3A)
txt.full=paste(txt1.1A,'\n',dims[1],'\n',dims[2],txt1.2A)
txt.full=paste(txt.full,text)
}
else{txt.full=paste(txt1.1,'\n',dims[1],'\n',dims[2],'\n',dims[3],'\n',dims[4],txt1.2)
txt.full=paste(txt.full,text)}
#write(txt.full,'result.cpp')
#library(Rcpp)
#library(RcppArmadillo)
if(wrt)
{
write(txt.full,'FPT_modelbuild.cpp')
}
stre="Compiling C++ code. Please wait."
cat(stre, " \r")
flush.console()
sourceCpp(code=txt.full)
cat(' ','\r')
N=(Tmax-s)/delt
tt=seq(s,Tmax-delt,delt)
y=ieq(solver(rep(B,N),B,c(0,theta),N,delt,N,tt),t(solver(rep(Xs,1),B,c(0,theta),N,delt,1,tt[1])),N,delt)
}
xcx = list(density=y,time=tt)
med.surv = function(x)
{
dens = x$density
time = x$time
dt = diff(time)[1]
cdf = cumsum(dens*dt)
if(any(cdf==0.5)){return(time[which(cdf==0.5)])}
wh.min = max(which(cdf < 0.5))
wh.pls = min(which(cdf > 0.5))
res = time[wh.min]+(0.5-cdf[wh.min])*(time[wh.pls]-time[wh.min])/(cdf[wh.pls]-cdf[wh.min])
return(res)
}
return(list(density=y,time=tt,prob.coverage = round(sum(y*delt),3),med.survival = med.surv(xcx)))
}
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.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GQD.TIpassage=
function(Xs,B,s,t,delt,theta=c(0),IEQ.type='Buonocore',wrt=FALSE)
{
fpt =function(res,res2,N,delt){}
ieq =function(res1,res2,N,delt){}
solver =function(Xs, Xt, theta, N , delt , N2, tt , P , alpha, lower , upper, tro ){}
rm(list =c('fpt','ieq','solver'))
Tmax=t
# Warning Module
if(delt>=0.1){warning("Large delt may result in poor approximations.");}
if(delt>=1) {stop("delt must be <1.");}
if(Xs>B) {stop("Xs must be < B. If the fpt-problem is a down-cross problem, consider the transform Yt = -Xt.");}
pow=function(x,p)
{
x^p
}
prod=function(a,b){a*b}
if(IEQ.type=='Volterra')
{
#stop("Volterra solution offline.");
#===============================================================================
# State 1
#===============================================================================
txtA.state1=
'
#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")]]
mat f(mat a,vec theta,vec t,int N2)
{
mat atemp(N2,2);'
txtB.state1.homo= '
return atemp;
}
// [[Rcpp::export]]
mat solver(vec Xs,vec Xt,vec theta,int N,double delt,int N2,vec tt)
{
mat fx0(N2,2);
mat fx1(N2,2);
mat fx2(N2,2);
mat fx3(N2,2);
mat fx4(N2,2);
mat fx5(N2,2);
mat fx6(N2,2);
mat fx7(N2,2);
mat fx8(N2,2);
mat fx9(N2,2);
mat fx10(N2,2);
mat fx11(N2,2);
mat fx12(N2,2);
mat fx13(N2,2);
mat fx14(N2,2);
mat fx15(N2,2);
mat fx16(N2,2);
mat x0(N2,2);
mat x1(N2,2);
mat x2(N2,2);
mat x3(N2,2);
mat x4(N2,2);
mat x5(N2,2);
mat x6(N2,2);
mat x7(N2,2);
mat x8(N2,2);
mat x9(N2,2);
mat x10(N2,2);
mat x11(N2,2);
mat x12(N2,2);
mat x13(N2,2);
mat x14(N2,2);
mat x15(N2,2);
mat x16(N2,2);
x0.fill(0);
x0.col(0)=Xs;
mat res(N2,N);
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;
res.col(i) =exp(-0.5*log(2*3.141592653589793*x0.col(1))-0.5*((Xt-x0.col(0))%(Xt-x0.col(0))/x0.col(1)));
}
return(res);
}'
#===============================================================================
# State 2
#===============================================================================
txtA.state2=
'
#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")]]
mat f(mat a,vec theta,vec t,int N2)
{
mat atemp(N2,4);'
txtB.state2.homo= '
return atemp;
}
// [[Rcpp::export]]
mat solver(vec Xs,vec Xt,vec theta,int N,double delt,int N2,vec tt)
{
mat fx0(N2,4);
mat fx1(N2,4);
mat fx2(N2,4);
mat fx3(N2,4);
mat fx4(N2,4);
mat fx5(N2,4);
mat fx6(N2,4);
mat fx7(N2,4);
mat fx8(N2,4);
mat fx9(N2,4);
mat fx10(N2,4);
mat fx11(N2,4);
mat fx12(N2,4);
mat fx13(N2,4);
mat fx14(N2,4);
mat fx15(N2,4);
mat fx16(N2,4);
mat x0(N2,4);
mat x1(N2,4);
mat x2(N2,4);
mat x3(N2,4);
mat x4(N2,4);
mat x5(N2,4);
mat x6(N2,4);
mat x7(N2,4);
mat x8(N2,4);
mat x9(N2,4);
mat x10(N2,4);
mat x11(N2,4);
mat x12(N2,4);
mat x13(N2,4);
mat x14(N2,4);
mat x15(N2,4);
mat x16(N2,4);
x0.fill(0);
x0.col(0)=Xs;
mat res(N2,N);
vec d=tt;
vec p(N2);
vec q(N2);
vec chk(N2);
vec th(N2);
vec K(N2);
vec K1(N2);
vec K2(N2);
vec val(N2);
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;
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);
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);
chk=pow(q,2)/4.0 + pow(p,3)/27.0;
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));
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;
K1 =x0.col(0) +(x0.col(1)%th) +(x0.col(2)%th%th)/2.0 +(x0.col(3)%th%th%th)/6.0;
K2 =x0.col(1) +(x0.col(2)%th) +(x0.col(3)%th%th)/2.0;
val=-0.5*log(2*3.141592653589793*K2)+(K-th%K1);
res.col(i)=exp(val);
}
return(res);
}'
code2=
'
// [[Rcpp::export]]
vec fpt(mat res1,vec res2,int N,double delt)
{
vec g(N);
g.fill(0);
double sm;
for (int i = 1; i < N; i++)
{
sm=0;
for (int j = 1; j < i; j++)
{
sm=sm+res1(j,i-j)*g(j);
}
g(i) = (res2(i)/delt -sm)/res1(i,1);
}
return g;
}
'
#==============================================================================
# Data resolution Module
#==============================================================================
Tmax=t
t=seq(0,100,by=1/100)
homo=T
homo.res=T
#==============================================================================
# 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)
state2=!state1
sol.state=c(state1,state2)
sol.state=c('2nd Ord. Truncation Normal Dist.','4th Ord. Truncation Saddlepoint Appr.')[which(sol.state==T)]
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)/N.mesh,diff(T.seq)/N.mesh)
#
#}
# REFERENCE MATRIX ------------------------------------------------------------
MAT=rbind(
c('1','a.col(0)' ,''),
c('','2*a.col(1)','1'),
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]))
#print(diff(result)[1])
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,3)]=1 # Always set these to 1
#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'}
for(i in 1:2)
{
dims[i]=paste0(paste0(paste0(' atemp.col(',i-1,')='),dims[i]),';')
}
# WRIGHT AND SOURCE -----------------------------------------------------------
if(homo.res)
{
txt.full=paste(txtA.state1,'\n',dims[1],'\n',dims[2],txtB.state1.homo,code2)
}
#if(!homo.res)
#{
# txt.full=paste(txtA.state1,'\n',dims[1],'\n',dims[2],txtB.state1.NONhomo,txtC.state1)
#}
#write(txt.full,'FPT_modelbuild.cpp')
type.sol =" Gneralized Ornstein-Uhlenbeck "#paste0("GENERALIZED LINEAR DIFFUSON",c('',": TIME INHOMOGNEOUS")[2-homo],c(')','+INDEPENDENT')[1+state1])
}
if(state2)
{
# DATA RESOLUTION -------------------------------------------------------------
#if(!homo.res)
#{
# delt=cbind(diff(T.seq)/N.mesh,diff(T.seq)/N.mesh,diff(T.seq)/N.mesh,diff(T.seq)/N.mesh)
#
#}
# REFERENCE MATRIX ------------------------------------------------------------
MAT=rbind(
c('1','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','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))'))
# 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]))
#print(diff(result)[1])
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,4)]=1 # Always set these to 1
#BUILD ODE --------------------------------------------------------------------
dims=rep('(',4)
for(i in 1:4)
{
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'}
for(i in 1:4)
{
dims[i]=paste0(paste0(paste0(' atemp.col(',i-1,')='),dims[i]),';')
}
# WRIGHT AND SOURCE -----------------------------------------------------------
if(homo.res)
{
txt.full=paste(txtA.state2,'\n',dims[1],'\n',dims[2],'\n',dims[3],'\n',dims[4],txtB.state2.homo,code2)
}
#if(!homo.res)
#{
# txt.full=paste(txtA.state2,'\n',dims[1],'\n',dims[2],'\n',dims[3],'\n',dims[4],txtB.state2.NONhomo,txtC.state2)
#}
#write(txt.full,'FPT_modelbuild.cpp')
type.sol =" Generalized Quadratic Diffusion (GQD) "#paste0("GENERALIZED LINEAR DIFFUSON",c('',": TIME INHOMOGNEOUS")[2-homo],c(')','+INDEPENDENT')[1+state1])
}
#library(Rcpp)
#library(RcppArmadillo)
if(wrt)
{
write(txt.full,'FPT_modelbuild.cpp')
}
stre="Compiling C++ code. Please wait."
cat(stre, " \r")
sourceCpp(code=txt.full)
cat(' ','\r')
N=(Tmax-s)/delt
tt=seq(s,Tmax-delt,delt)
res=solver(rep(B,N),rep(B,N),c(0,theta),N,delt,N,tt)
res2=solver(rep(Xs,1),rep(B,1),c(0,theta),N,delt,1,0)
y=fpt(res,res2,N,delt)
#return(list(density=y,time=tt))
#print('Done.')
}
################################################################################
if(IEQ.type=='Buonocore')
{
#===============================================================================
# State 1
#===============================================================================
#_______________________________________________________________________________
txt1.1A= '
#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);'
#_______________________________________________________________________________
txt1.2A= '
return atemp;
}
'
#_______________________________________________________________________________
txt2A='
// [[Rcpp::export]]
vec saddle(mat xx,double Bt)
{
vec val=exp(-0.5*log(2*3.141592653589793*xx.col(1))-0.5*(Bt-xx.col(0))%(Bt-xx.col(0))/xx.col(1));
return(val);
}
// [[Rcpp::export]]
vec saddle2(mat xx,double Bt)
{
vec val=-exp(-0.5*log(2*3.141592653589793*xx.col(1))-0.5*(Bt-xx.col(0))%(Bt-xx.col(0))/xx.col(1))%(Bt-xx.col(0))/xx.col(1);
return(val);
}
// [[Rcpp::export]]
vec pcurr(mat xx,double Bt,vec theta,vec t)
{
vec val ='
#_______________________________________________________________________________
txt3A=';
return(val);
}
// [[Rcpp::export]]
mat solver(vec Xs,double Bt,vec theta,int N,double delt,int N2,vec tt)
{
mat fx0(N2,2);
mat fx1(N2,2);
mat fx2(N2,2);
mat fx3(N2,2);
mat fx4(N2,2);
mat fx5(N2,2);
mat fx6(N2,2);
mat fx7(N2,2);
mat fx8(N2,2);
mat fx9(N2,2);
mat fx10(N2,2);
mat fx11(N2,2);
mat fx12(N2,2);
mat fx13(N2,2);
mat fx14(N2,2);
mat fx15(N2,2);
mat fx16(N2,2);
mat x0(N2,2);
mat x1(N2,2);
mat x2(N2,2);
mat x3(N2,2);
mat x4(N2,2);
mat x5(N2,2);
mat x6(N2,2);
mat x7(N2,2);
mat x8(N2,2);
mat x9(N2,2);
mat x10(N2,2);
mat x11(N2,2);
mat x12(N2,2);
mat x13(N2,2);
mat x14(N2,2);
mat x15(N2,2);
mat x16(N2,2);
x0.fill(0);
x0.col(0)=Xs;
mat res(N2,N);
res.fill(0);
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;
res.col(i) = pcurr(x0,Bt,theta,d);
}
return(res);
}
// [[Rcpp::export]]
vec ieq(mat res1,vec res2,int N,double delt)
{
vec v(N);
v(0)=0;
double smtemp=0;
for (int i = 1; i < N; i++)
{
smtemp=0;
for (int j = 1; j < i; j++)
{
smtemp=smtemp+v(j)*res1(j,i-j);
}
v(i) = (2.0*res2(i)-2.0*smtemp*delt);
}
return(v);
}
'
#===============================================================================
# State 2
#===============================================================================
#_______________________________________________________________________________
txt1.1= '
#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);'
#_______________________________________________________________________________
txt1.2= '
return atemp;
}
'
#_______________________________________________________________________________
txt2='
// [[Rcpp::export]]
vec saddle(mat xx,double Bt)
{
vec p=(1.0/3.0) *(3*(xx.col(3)/6.0)%xx.col(1) - pow(xx.col(2)/2.0,2))/pow(xx.col(3)/6.0,2);
vec q=(1.0/27.0)*(27*pow(xx.col(3)/6.0,2)%(xx.col(0)-Bt) - 9*(xx.col(3)/6.0)%(xx.col(2)/2.0)%xx.col(1) + 2*pow(xx.col(2)/2.0,3))/pow(xx.col(3)/6.0,3);
vec chk=pow(q,2)/4.0 + pow(p,3)/27.0;
vec th=-(xx.col(2)/2.0)/(3*(xx.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 =xx.col(0)%th+(xx.col(1)%th%th)/2.0+(xx.col(2)%th%th%th)/6.0 +(xx.col(3)%th%th%th%th)/24.0;
vec K1=xx.col(0) +(xx.col(1)%th) +(xx.col(2)%th%th)/2.0 +(xx.col(3)%th%th%th)/6.0;
vec K2=xx.col(1) +(xx.col(2)%th) +(xx.col(3)%th%th)/2.0;
vec val=exp(-0.5*log(2*3.141592653589793*K2)+(K-th%K1));
return(val);
}
// [[Rcpp::export]]
vec saddle2(mat xx,double Bt)
{
vec p=(1.0/3.0) *(3*(xx.col(3)/6.0)%xx.col(1) - pow(xx.col(2)/2.0,2))/pow(xx.col(3)/6.0,2);
vec q=(1.0/27.0)*(27*pow(xx.col(3)/6.0,2)%(xx.col(0)-Bt) - 9*(xx.col(3)/6.0)%(xx.col(2)/2.0)%xx.col(1) + 2*pow(xx.col(2)/2.0,3))/pow(xx.col(3)/6.0,3);
vec chk=pow(q,2)/4.0 + pow(p,3)/27.0;
vec th=-(xx.col(2)/2.0)/(3*(xx.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 thdash=1.0/(xx.col(1)+th%xx.col(2)+0.5*th%th%xx.col(3));
vec K =xx.col(0)%th+(xx.col(1)%th%th)/2.0+(xx.col(2)%th%th%th)/6.0 +(xx.col(3)%th%th%th%th)/24.0;
vec K1=xx.col(0) +(xx.col(1)%th) +(xx.col(2)%th%th)/2.0 +(xx.col(3)%th%th%th)/6.0;
vec K2=xx.col(1) +(xx.col(2)%th) +(xx.col(3)%th%th)/2.0;
vec gg = 1.0/sqrt(2*3.141592653589793*(K2));
vec ggdash = -3.141592653589793*pow(2*3.141592653589793*(K2),-3.0/2.0)%(xx.col(2)%thdash+xx.col(3)%thdash%th);
vec hh = exp(K-th%K1);
vec hhdash = exp(K-th%K1)%(-th%thdash%xx.col(1)-th%th%thdash%xx.col(2)-0.5*(th%th%th)%thdash%xx.col(3));
vec val=(ggdash%hh+gg%hhdash);
return(val);
}
// [[Rcpp::export]]
vec pcurr(mat xx,double Bt,vec theta,vec t)
{
vec val ='
#_______________________________________________________________________________
txt3=';
return(val);
}
// [[Rcpp::export]]
mat solver(vec Xs,double Bt,vec theta,int N,double delt,int N2,vec tt)
{
mat fx0(N2,4);
mat fx1(N2,4);
mat fx2(N2,4);
mat fx3(N2,4);
mat fx4(N2,4);
mat fx5(N2,4);
mat fx6(N2,4);
mat fx7(N2,4);
mat fx8(N2,4);
mat fx9(N2,4);
mat fx10(N2,4);
mat fx11(N2,4);
mat fx12(N2,4);
mat fx13(N2,4);
mat fx14(N2,4);
mat fx15(N2,4);
mat fx16(N2,4);
mat x0(N2,4);
mat x1(N2,4);
mat x2(N2,4);
mat x3(N2,4);
mat x4(N2,4);
mat x5(N2,4);
mat x6(N2,4);
mat x7(N2,4);
mat x8(N2,4);
mat x9(N2,4);
mat x10(N2,4);
mat x11(N2,4);
mat x12(N2,4);
mat x13(N2,4);
mat x14(N2,4);
mat x15(N2,4);
mat x16(N2,4);
x0.fill(0);
x0.col(0)=Xs;
mat res(N2,N);
res.fill(0);
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;
res.col(i) = pcurr(x0,Bt,theta,d);
}
return(res);
}
// [[Rcpp::export]]
vec ieq(mat res1,vec res2,int N,double delt)
{
vec v(N);
v(0)=0;
double smtemp=0;
for (int i = 1; i < N; i++)
{
smtemp=0;
for (int j = 1; j < i; j++)
{
smtemp=smtemp+v(j)*res1(j,i-j);
}
v(i) = (2.0*res2(i)-2.0*smtemp*delt);
}
return(v);
}
'
# BUILD PROB. CURRENT ---------------------------------------------------------
t=seq(0,100,by=1/100)
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}
}
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))
}
T1=''
T2=''
T3=''
if(func.list[1]==1)
{
T1=paste0(T1,'+(',body(namess[1])[2],')')
}
if(func.list[2]==1)
{
T1=paste0(T1,'+(',body(namess[2])[2],')*Bt')
}
if(func.list[3]==1)
{
T1=paste0(T1,'+(',body(namess[3])[2],')*Bt*Bt')
}
if(func.list[4]==1)
{
T3=paste0(T3,'+(',body(namess[4])[2],')')
}
if(func.list[5]==1)
{
T2=paste0(T2,'+(',body(namess[5])[2],')')
T3=paste0(T3,'+(',body(namess[5])[2],')*Bt')
}
if(func.list[6]==1)
{
T2=paste0(T2,'+2*(',body(namess[6])[2],')*Bt')
T3=paste0(T3,'+(',body(namess[6])[2],')*Bt*Bt')
}
T1=paste0('(',T1,')')
T2=paste0('(',T2,')')
T3=paste0('(',T3,')')
T1.inhom=F
T2.inhom=F
T3.inhom=F
if(any(func.list.timehomo[c(1,2,3)]==2)){T1.inhom=T}
if(any(func.list.timehomo[c(5,6)]==2)) {T2.inhom=T}
if(any(func.list.timehomo[c(4,5,6)]==2)){T3.inhom=T}
text=paste0(txt2,'0.5*',T1,c('*','%')[T1.inhom+1],'saddle(xx,Bt)-','0.75*',T2,c('*','%')[T2.inhom+1],'saddle(xx,Bt)-','0.5*',T3,c('*','%')[T3.inhom+1],'saddle2(xx,Bt)',txt3)
func.list.timehomo[c(1,4)]=1 # Always set these to 1
# REFERENCE MATRIX ------------------------------------------------------------
MAT=rbind(
c('1','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','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))'))
#BUILD ODE --------------------------------------------------------------------
dims=rep('(',4)
for(i in 1:4)
{
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=='()')]='0'}
for(i in 1:4)
{
dims[i]=paste0(paste0(paste0(' atemp.col(',i-1,')='),dims[i]),';')
}
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])
}
# WRIGHT AND SOURCE -----------------------------------------------------------
if(sum(func.list[c(3,5,6)]==0)==3)
{
text=paste0(txt2A,'0.5*',T1,c('*','%')[T1.inhom+1],'saddle(xx,Bt)-','0.5*',T3,c('*','%')[T3.inhom+1],'saddle2(xx,Bt)',txt3A)
txt.full=paste(txt1.1A,'\n',dims[1],'\n',dims[2],txt1.2A)
txt.full=paste(txt.full,text)
}
else{txt.full=paste(txt1.1,'\n',dims[1],'\n',dims[2],'\n',dims[3],'\n',dims[4],txt1.2)
txt.full=paste(txt.full,text)}
#write(txt.full,'result.cpp')
#library(Rcpp)
#library(RcppArmadillo)
if(wrt)
{
write(txt.full,'FPT_modelbuild.cpp')
}
stre="Compiling C++ code. Please wait."
cat(stre, " \r")
flush.console()
sourceCpp(code=txt.full)
cat(' ','\r')
N=(Tmax-s)/delt
tt=seq(s,Tmax-delt,delt)
y=ieq(solver(rep(B,N),B,c(0,theta),N,delt,N,tt),t(solver(rep(Xs,1),B,c(0,theta),N,delt,1,tt[1])),N,delt)
}
xcx = list(density=y,time=tt)
med.surv = function(x)
{
dens = x$density
time = x$time
dt = diff(time)[1]
cdf = cumsum(dens*dt)
if(any(cdf==0.5)){return(time[which(cdf==0.5)])}
wh.min = max(which(cdf < 0.5))
wh.pls = min(which(cdf > 0.5))
res = time[wh.min]+(0.5-cdf[wh.min])*(time[wh.pls]-time[wh.min])/(cdf[wh.pls]-cdf[wh.min])
return(res)
}
return(list(density=y,time=tt,prob.coverage = round(sum(y*delt),3),med.survival = med.surv(xcx)))
}
Try the DiffusionRgqd package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.