(*Timothy I. Matis, Ph.D. and Ivan G. Guardiola, Department of Industrial Engineering, Texas Tech University, April 2006*)
BeginPackage["MomentClosure`",{"Combinatorica`"}];
Print["The functions MomentClosureSystem and MomentClosurePlots are defined in this package for the transient analysis of the cumulants (moments) of a non-linear stochastic compartmental model with Markov transitions using moment closure with cumulant neglect. In particular, MomentClosureSystem returns a closed system of approximating ordinary differential equations for the cumulants, and MomentClosurePlots returns the cumulant approximations over a finite time horizon. Further usage information may be obtained through ?MomentClosureSystem and/or ?MomentClosurePlots."]
MomentClosurePlots::usage="The function MomentClosurePlots returns closure approximations to the cumulants of a numerically defined nonlinear stochastic compartmental model with Markov transitions over a specified time horizon. Inputs to this function include the intensity functions of the system, initial conditions for the cumulants, the neglect level, the time horizon, cumulants for which output is desired, and the plot range of those cumulants. Hence, the function is entered as, MomentClosurePlots[Intensities,InitialConditions,NeglectLevel,TimeRange,PlotFunctions,PlotRanges], which may be expressed in general syntax as,
MomentClosurePlots[{{{b__},f(x__)},...},{K[c__][0]==Value_,...},m,t,{K[c__][t],...},{{lower_,upper_},...}]. As an example, the following syntax would be entered to approximate the mean and variance of a single compartment birth death model with a birth rate of 1 and a death rate of .1*x^2 between time (0,10) under a neglect level of 3 assuming that all cumulant are intially set to zero, MomentClosurePlots[{{{1},1},{{-1},.1*x^2}},{},3,10,{K[1][t],K[2][t]},{}]. Please note that the function MomentClosurePlots does not check for the stability of the moment closure differential equations apriori, which may result in error messages during evaluation if the equations for the defined system are unstable. In addition, specific usage information for the input parameters Intensities, InitialConditions, NeglectLevel, TimeRange, PlotFunctions, and PlotRanges may be obtained through ?ParameterName."
MomentClosureSystem::usage="The function MomentClosureSystem returns a closed set of ordinary differential equations that describes the closure based cumulants of a nonlinear stochastic compartmental model with Markov transitions through the neglect level. Inputs to this function include the intensity functions of the system and the desired level of cumulant neglect. Hence, the function is entered as MomentClosureSystem[Intensities_,NeglectLevel_], which may be expressed in general syntax as MomentClosureSystem[{{{b__},f(x__)},...},m]. As an example, the following syntax would be entered to generate the moment closure differential equations for a single compartment system under a neglect level of 3 with a birth rate of 1 and a death rate of .1*x^2, MomentClosureSystem[{{{1},1},{{-1},.1*Subscript[x,1]^2}},3]. Moment closure differential equations may also be generated for symbolically defined systems by subsituting a greek letter for the numerical values. As a continuation of the example, the same system defined by a birth rate of \[Lambda] and a death rate of \[Mu]*x^2 would be entered as MomentClosureSystem[{{{1},\[Lambda]},{{-1},\[Mu]*Subscript[x,1]^2}},3]. Please note that specific usage information for the input parameters Intensities and NeglectLevel may be attained through ?ParameterName."
Intensities::usage="This parameter is a list of two-tuples that describes the possible changes in the state of the system with the corresponding intensity functions. For example, a two-compartment system with immigration intensity \[Lambda], transiton intensity \[Mu]*Subscript[x,1]^2, and death intensity \[Gamma]*(Subscript[x,2])^3 would be entered as {{{1,0},\[Lambda]},{{-1,1},\[Mu]*Subscript[x,1]^2},{{0,-1},\[Gamma]*Subscript[x,2]^3}}. Note that for symbolically defined systems, all symbols must be in greek letters. In addition, please note that all variables must be subscripted with the compartment number, and there is a limit of 25 compartments. "
InitialConditions::usage="This parameter is a list of those cumulants that intially have a non-zero value. The default is the empty list which assumes all cumulants are intially zero. For example, if there were 2 entities in the first compartment of a two-compartment system at time zero, the InitialConditions list would be entered as {K[1,0][0]==2}. As another example, the InitialConditions list for a two-compartment system that is empty a time zero would be {}."
NeglectLevel::usage="This parameter specifies the moment closure cumulant neglect level. All cumulants (moments) greater than this will be set to zero, thereby closing the differential equations. The neglect level must be an integer value greater than 1."
TimeRange::usage="This parameter specifies the time horizon over which cumulant approximations are obtained. In particular, it specifies the endpoint of the finite interval (0,t)."
PlotFunctions::usage="This parameter is a list that denotes those cumulant approximations whose graphical output will be displayed. For example, the following syntax for PlotFunctions would yield a graphical plot of the moment closure approximations for the means, variances, and covariance of a two-compartment system, {K[1,0][t],K[0,1][t],K[2,0][t],K[0,2][t],K[1,1][t]}."
PlotRanges::usage="This parameter is a list of one-to-one correspondence with PlotFunctions that denotes the range of the y-axis for which the cumulant approximations will be displayed. For the example of a two-compartment system given under the help for PlotFunctions, the following syntax for PlotRanges would yield graphs whose y-axis is ranges from (0,5) for the means, (10,20) for the variances, and (-5,5) for the covariance, {{0,5},{0,5},{10,20},{10,20},{-5,5}}. The default is the empty list for which the range is automatically specified by the Mathematica program."
LHSPDE::usage="This module builds the left hand side of the partial differential equation. It returns a general form of the cumulant generating function of correct dimensionality for the system and the left hand side of the partial differential equation of the cumulant generating function of the system. "
EQUATE::usage="This module builds the right hand side of the partial differential equation of the cumulant generating function and sets it equal to the left hand side generated by the module LHSPDE. All partial derivatives and exponentials of the partial differential equation are expanded to a polynomial form, and the coefficients of powered combinations of the dummy variables theta are matched. The module returns the equated differential equations."
INITIAL::usage="This function sets the differential equations obtained from EQUATE into the list format necessary for NDSolve and appends initial conditions to this list. The module returns a list that is ready to be input into the NDSolve function."
NUMERICALSOLUTION::usage="This module passes the list obtained from the module INITIAL to the NDSolve command, upon which the command is executed through the user-specified time horizon. The module returns the interpolating functions of the cumulants obtained from the numerical solution of the ordinary differential equations."
LHSPDE[Nodes_,TRUNC_]:=
Module[{AA,AB,AC,AD,AE,AF,AG,AH,AI,AJ,AK,AL,AM,AN,AO,AP,AQ,AR,AS,AT,AU,AV,
ORLI,MGFList,bar,CGFlist,IVV,Cgf, MGF,LHSpde,orderList1},
AA=Table[Table[0,{Nodes}],{Nodes}];AB=MapThread[K,Transpose[AA]];
Clear[AA];
AD=AB/.{K[b__]\[Rule]K[b][t]};AE=Take[First[AD]];
Evaluate[AE]==0;AG={a,b,c,d,e,f,g,h,r,s,t,u,v,w,x,y,z};
AH=Take[AG,Nodes];AI=Table[0,{Nodes}];AJ=Table[TRUNC,{Nodes}];
AK=Transpose[{AH,AI,AJ}];AM=MapThread[K,AH,0];
AN=Table[Subscript[\[Theta],i],{i,Nodes}];AO=AN^AH;AP=Times@@AO;
AQ=Factorial[AH];AR=Times@@AQ;
AS=AP/AR;AT=AM/.{K[z__]\[Rule]K[z][t]};AV=AT*AS;ORLI={};
For[i=1,iorderList1[[i]]],{i,1,Length[orderList1]}];
butt=Table[Thread[AN->Sign[orderList1[[i]]]],{i,1,Length[orderList1]}];
Y12={};
butt2=
Table[Table[
If[Sign[orderList1[[j,i]]]\[Equal]1,a,butt[[j,i]]],{i,1,
Length[orderList1[[1]]]}],{j,1,Length[orderList1]}];
For[i=1,i