Example 11 with VariableHash
use of cbit.vcell.math.VariableHash in project vcell by virtualcell.
the class MathTestingUtilities method constructOdesForSensitivity.
//
// This method is used to solve for sensitivity of variables to a given parameter.
// The mathDescription and the sensitivity parameter are passed as arguments.
// New variables and ODEs are constructed according to the rule listed below and are added to the mathDescription.
// The method returns the modified mathDescription.
//
public static MathDescription constructOdesForSensitivity(MathDescription mathDesc, Constant sensParam) throws ExpressionException, MathException, MappingException {
//
// For each ODE :
//
// dX/dt = F(X, P)
//
// (where P is the sensitivity parameter)
//
// we create two other ODEs :
//
// dX_1/dt = F(X_1, P_1) and
//
// dX_2/dt = F(X_2, P_2)
//
// where P_1 = P + epsilon, and
// P_2 = P - epsilon.
//
// We keep the initial conditions for both the new ODEs to be the same,
// i.e., X_1_init = X_2_init.
//
// Then, solving for X_1 & X_2, sensitivity of X wrt P can be computed as :
//
// dX = (X_1 - X_2)
// -- ----------- .
// dP (P_1 - P_2)
//
//
// REMOVE PRINTS AFTER CHECKING !!!
System.out.println(" \n\n------------ Old Math Description -----------------");
System.out.println(mathDesc.getVCML_database());
if (mathDesc.getGeometry().getDimension() > 0) {
throw new RuntimeException("Suppport for Spatial systems not yet implemented.");
}
VariableHash varHash = new VariableHash();
Enumeration<Variable> enumVar = mathDesc.getVariables();
while (enumVar.hasMoreElements()) {
varHash.addVariable(enumVar.nextElement());
}
//
// Get 2 values of senstivity parameter (P + epsilon) & (P - epsilon)
//
Constant epsilon = new Constant("epsilon", new Expression(sensParam.getConstantValue() * 1e-3));
Constant sensParam1 = new Constant(sensParam.getName() + "_1", new Expression(sensParam.getConstantValue() + epsilon.getConstantValue()));
Constant sensParam2 = new Constant(sensParam.getName() + "_2", new Expression(sensParam.getConstantValue() - epsilon.getConstantValue()));
//
// Iterate through each subdomain (only 1 in compartmental case), and each equation in the subdomain
//
Enumeration<SubDomain> subDomainEnum = mathDesc.getSubDomains();
//
// Create a vector of equations to store the 2 equations for each ODE variable in the subdomain.
// Later, add it to the equations list in the subdomain.
//
Vector<Equation> equnsVector = new Vector<Equation>();
Vector<Variable> varsVector = new Vector<Variable>();
Vector<Variable> var1s = new Vector<Variable>();
Vector<Variable> var2s = new Vector<Variable>();
while (subDomainEnum.hasMoreElements()) {
SubDomain subDomain = subDomainEnum.nextElement();
Enumeration<Equation> equationEnum = subDomain.getEquations();
Domain domain = new Domain(subDomain);
while (equationEnum.hasMoreElements()) {
Equation equation = equationEnum.nextElement();
if (equation instanceof OdeEquation) {
OdeEquation odeEquation = (OdeEquation) equation;
// Similar to substituteWithExactSolutions, to bind and substitute functions in the ODE
Expression substitutedRateExp = substituteFunctions(odeEquation.getRateExpression(), mathDesc);
String varName = odeEquation.getVariable().getName();
VolVariable var = new VolVariable(varName, domain);
varsVector.addElement(var);
//
// Create the variable var1, and get the initExpr and rateExpr from the original ODE.
// Substitute the new vars (var1 and param1) in the old initExpr and rateExpr and create a new ODE
//
String varName1 = new String("__" + varName + "_1");
Expression initExpr1 = odeEquation.getInitialExpression();
Expression rateExpr1 = new Expression(substitutedRateExp);
rateExpr1.substituteInPlace(new Expression(varName), new Expression(varName1));
rateExpr1.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam1.getName()));
VolVariable var1 = new VolVariable(varName1, domain);
var1s.addElement(var1);
OdeEquation odeEqun1 = new OdeEquation(var1, initExpr1, rateExpr1);
equnsVector.addElement(odeEqun1);
//
// Create the variable var2, and get the initExpr and rateExpr from the original ODE.
// Substitute the new vars (var2 and param2) in the old initExpr and rateExpr and create a new ODE
//
String varName2 = new String("__" + varName + "_2");
Expression initExpr2 = odeEquation.getInitialExpression();
Expression rateExpr2 = new Expression(substitutedRateExp);
rateExpr2.substituteInPlace(new Expression(varName), new Expression(varName2));
rateExpr2.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam2.getName()));
VolVariable var2 = new VolVariable(varName2, domain);
var2s.addElement(var2);
OdeEquation odeEqun2 = new OdeEquation(var2, initExpr2, rateExpr2);
equnsVector.addElement(odeEqun2);
//
// Create a function for the sensitivity function expression (X1-X2)/(P1-P2), and save in varHash
//
Expression diffVar = Expression.add(new Expression(var1.getName()), Expression.negate(new Expression(var2.getName())));
Expression diffParam = Expression.add(new Expression(sensParam1.getName()), Expression.negate(new Expression(sensParam2.getName())));
Expression sensitivityExpr = Expression.mult(diffVar, Expression.invert(diffParam));
Function sens_Func = new Function("__sens" + varName + "_wrt_" + sensParam.getName(), sensitivityExpr, domain);
varHash.addVariable(epsilon);
varHash.addVariable(sensParam1);
varHash.addVariable(sensParam2);
varHash.addVariable(var1);
varHash.addVariable(var2);
varHash.addVariable(sens_Func);
} else {
// sensitivity not implemented for PDEs or other equation types.
throw new RuntimeException("SolverTest.constructedExactMath(): equation type " + equation.getClass().getName() + " not yet implemented");
}
}
//
// Need to substitute the new variables in the new ODEs.
// i.e., if Rate Expr for ODE_1 for variable X_1 contains variable Y, variable Z, etc.
// Rate Expr is already substituted with X_1, but it also needs substitute Y with Y_1, Z with Z_1, etc.
// So get the volume variables, from the vectors for vars, var_1s and var_2s
// Substitute the rate expressions for the newly added ODEs in equnsVector.
//
Variable[] vars = (Variable[]) BeanUtils.getArray(varsVector, Variable.class);
Variable[] var_1s = (Variable[]) BeanUtils.getArray(var1s, Variable.class);
Variable[] var_2s = (Variable[]) BeanUtils.getArray(var2s, Variable.class);
Vector<Equation> newEqunsVector = new Vector<Equation>();
for (int i = 0; i < equnsVector.size(); i++) {
Equation equn = equnsVector.elementAt(i);
Expression initEx = equn.getInitialExpression();
Expression rateEx = equn.getRateExpression();
for (int j = 0; j < vars.length; j++) {
if (equn.getVariable().getName().endsWith("_1")) {
rateEx.substituteInPlace(new Expression(vars[j].getName()), new Expression(var_1s[j].getName()));
} else if (equn.getVariable().getName().endsWith("_2")) {
rateEx.substituteInPlace(new Expression(vars[j].getName()), new Expression(var_2s[j].getName()));
}
}
OdeEquation odeEqun = new OdeEquation(equn.getVariable(), initEx, rateEx);
newEqunsVector.addElement(odeEqun);
}
//
for (int i = 0; i < newEqunsVector.size(); i++) {
mathDesc.getSubDomain(subDomain.getName()).addEquation((Equation) newEqunsVector.elementAt(i));
}
//
// FAST SYSTEM
// If the subdomain has a fast system, create a new fast system by substituting the high-low variables/parameters
// in the expressions for the fastInvariants and fastRates and adding them to the fast system.
//
Vector<FastInvariant> invarsVector = new Vector<FastInvariant>();
Vector<FastRate> ratesVector = new Vector<FastRate>();
Enumeration<FastInvariant> fastInvarsEnum = null;
Enumeration<FastRate> fastRatesEnum = null;
// Get the fast invariants and fast rates in the system.
FastSystem fastSystem = subDomain.getFastSystem();
if (fastSystem != null) {
fastInvarsEnum = fastSystem.getFastInvariants();
fastRatesEnum = fastSystem.getFastRates();
//
while (fastInvarsEnum.hasMoreElements()) {
FastInvariant fastInvar = fastInvarsEnum.nextElement();
Expression fastInvarExpr = fastInvar.getFunction();
fastInvarExpr = MathUtilities.substituteFunctions(fastInvarExpr, mathDesc);
Expression fastInvarExpr1 = new Expression(fastInvarExpr);
Expression fastInvarExpr2 = new Expression(fastInvarExpr);
for (int i = 0; i < vars.length; i++) {
fastInvarExpr1.substituteInPlace(new Expression(vars[i].getName()), new Expression(var_1s[i].getName()));
fastInvarExpr2.substituteInPlace(new Expression(vars[i].getName()), new Expression(var_2s[i].getName()));
}
fastInvarExpr1.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam1.getName()));
FastInvariant fastInvar1 = new FastInvariant(fastInvarExpr1);
invarsVector.addElement(fastInvar1);
fastInvarExpr2.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam2.getName()));
FastInvariant fastInvar2 = new FastInvariant(fastInvarExpr2);
invarsVector.addElement(fastInvar2);
}
// Add the newly created fast invariants to the existing list of fast invariants in the fast system.
for (int i = 0; i < invarsVector.size(); i++) {
FastInvariant inVar = (FastInvariant) invarsVector.elementAt(i);
fastSystem.addFastInvariant(inVar);
}
//
while (fastRatesEnum.hasMoreElements()) {
FastRate fastRate = fastRatesEnum.nextElement();
Expression fastRateExpr = fastRate.getFunction();
fastRateExpr = MathUtilities.substituteFunctions(fastRateExpr, mathDesc);
Expression fastRateExpr1 = new Expression(fastRateExpr);
Expression fastRateExpr2 = new Expression(fastRateExpr);
for (int i = 0; i < vars.length; i++) {
fastRateExpr1.substituteInPlace(new Expression(vars[i].getName()), new Expression(var_1s[i].getName()));
fastRateExpr2.substituteInPlace(new Expression(vars[i].getName()), new Expression(var_2s[i].getName()));
}
fastRateExpr1.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam1.getName()));
FastRate fastRate1 = new FastRate(fastRateExpr1);
ratesVector.addElement(fastRate1);
fastRateExpr2.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam2.getName()));
FastRate fastRate2 = new FastRate(fastRateExpr2);
ratesVector.addElement(fastRate2);
}
// Add the newly created fast rates to the existing list of fast rates in the fast system.
for (int i = 0; i < ratesVector.size(); i++) {
FastRate rate = (FastRate) ratesVector.elementAt(i);
fastSystem.addFastRate(rate);
}
}
}
// Reset all variables in mathDesc.
mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
// REMOVE PRINTS AFTER CHECKING
System.out.println(" \n\n------------ New Math Description -----------------");
System.out.println(mathDesc.getVCML_database());
return mathDesc;
}
Also used :
MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable)
InsideVariable(cbit.vcell.math.InsideVariable)
SensVariable(cbit.vcell.solver.ode.SensVariable)
FilamentVariable(cbit.vcell.math.FilamentVariable)
VolVariable(cbit.vcell.math.VolVariable)
VolumeRegionVariable(cbit.vcell.math.VolumeRegionVariable)
ReservedVariable(cbit.vcell.math.ReservedVariable)
MemVariable(cbit.vcell.math.MemVariable)
OutsideVariable(cbit.vcell.math.OutsideVariable)
FilamentRegionVariable(cbit.vcell.math.FilamentRegionVariable)
Variable(cbit.vcell.math.Variable)
VariableHash(cbit.vcell.math.VariableHash)
Constant(cbit.vcell.math.Constant)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
Function(cbit.vcell.math.Function)
Vector(java.util.Vector)
VolVariable(cbit.vcell.math.VolVariable)
OdeEquation(cbit.vcell.math.OdeEquation)
PdeEquation(cbit.vcell.math.PdeEquation)
Equation(cbit.vcell.math.Equation)
FastRate(cbit.vcell.math.FastRate)
FastInvariant(cbit.vcell.math.FastInvariant)
OdeEquation(cbit.vcell.math.OdeEquation)
Expression(cbit.vcell.parser.Expression)
FastSystem(cbit.vcell.math.FastSystem)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
Domain(cbit.vcell.math.Variable.Domain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
Example 12 with VariableHash
use of cbit.vcell.math.VariableHash in project vcell by virtualcell.
the class MatlabOdeFileCoder method write_V6_MFile.
/**
* Insert the method's description here.
* Creation date: (3/8/00 10:31:52 PM)
*/
public void write_V6_MFile(java.io.PrintWriter pw, String functionName) throws MathException, ExpressionException {
MathDescription mathDesc = simulation.getMathDescription();
if (!mathDesc.isValid()) {
throw new MathException("invalid math description\n" + mathDesc.getWarning());
}
if (mathDesc.isSpatial()) {
throw new MathException("spatial math description, cannot create ode file");
}
if (mathDesc.hasFastSystems()) {
throw new MathException("math description contains algebraic constraints, cannot create .m file");
}
//
// print function declaration
//
pw.println("function [T,Y,yinit,param, allNames, allValues] = " + functionName + "(argTimeSpan,argYinit,argParam)");
pw.println("% [T,Y,yinit,param] = " + functionName + "(argTimeSpan,argYinit,argParam)");
pw.println("%");
pw.println("% input:");
pw.println("% argTimeSpan is a vector of start and stop times (e.g. timeSpan = [0 10.0])");
pw.println("% argYinit is a vector of initial conditions for the state variables (optional)");
pw.println("% argParam is a vector of values for the parameters (optional)");
pw.println("%");
pw.println("% output:");
pw.println("% T is the vector of times");
pw.println("% Y is the vector of state variables");
pw.println("% yinit is the initial conditions that were used");
pw.println("% param is the parameter vector that was used");
pw.println("% allNames is the output solution variable names");
pw.println("% allValues is the output solution variable values corresponding to the names");
pw.println("%");
pw.println("% example of running this file: [T,Y,yinit,param,allNames,allValues] = myMatlabFunc; <-(your main function name)");
pw.println("%");
VariableHash varHash = new VariableHash();
for (Variable var : simulationSymbolTable.getVariables()) {
varHash.addVariable(var);
}
Variable[] variables = varHash.getTopologicallyReorderedVariables();
CompartmentSubDomain subDomain = (CompartmentSubDomain) mathDesc.getSubDomains().nextElement();
//
// collect "true" constants (Constants without identifiers)
//
//
// collect "variables" (VolVariables only)
//
//
// collect "functions" (Functions and Constants with identifiers)
//
Vector<Constant> constantList = new Vector<Constant>();
Vector<VolVariable> volVarList = new Vector<VolVariable>();
Vector<Variable> functionList = new Vector<Variable>();
for (int i = 0; i < variables.length; i++) {
if (variables[i] instanceof Constant) {
Constant constant = (Constant) variables[i];
String[] symbols = constant.getExpression().getSymbols();
if (symbols == null || symbols.length == 0) {
constantList.addElement(constant);
} else {
functionList.add(constant);
}
} else if (variables[i] instanceof VolVariable) {
volVarList.addElement((VolVariable) variables[i]);
} else if (variables[i] instanceof Function) {
functionList.addElement(variables[i]);
}
}
Constant[] constants = (Constant[]) BeanUtils.getArray(constantList, Constant.class);
VolVariable[] volVars = (VolVariable[]) BeanUtils.getArray(volVarList, VolVariable.class);
Variable[] functions = (Variable[]) BeanUtils.getArray(functionList, Variable.class);
int numVars = volVarList.size() + functionList.size();
String varNamesForStringArray = "";
String varNamesForValueArray = "";
for (Variable var : volVarList) {
varNamesForStringArray = varNamesForStringArray + "'" + var.getName() + "';";
varNamesForValueArray = varNamesForValueArray + var.getName() + " ";
}
for (Variable func : functionList) {
varNamesForStringArray = varNamesForStringArray + "'" + func.getName() + "';";
varNamesForValueArray = varNamesForValueArray + func.getName() + " ";
}
pw.println("");
pw.println("%");
pw.println("% Default time span");
pw.println("%");
double beginTime = 0.0;
double endTime = simulation.getSolverTaskDescription().getTimeBounds().getEndingTime();
pw.println("timeSpan = [" + beginTime + " " + endTime + "];");
pw.println("");
pw.println("% output variable lengh and names");
pw.println("numVars = " + numVars + ";");
pw.println("allNames = {" + varNamesForStringArray + "};");
pw.println("");
pw.println("if nargin >= 1");
pw.println("\tif length(argTimeSpan) > 0");
pw.println("\t\t%");
pw.println("\t\t% TimeSpan overridden by function arguments");
pw.println("\t\t%");
pw.println("\t\ttimeSpan = argTimeSpan;");
pw.println("\tend");
pw.println("end");
pw.println("%");
pw.println("% Default Initial Conditions");
pw.println("%");
pw.println("yinit = [");
for (int j = 0; j < volVars.length; j++) {
Expression initial = subDomain.getEquation(volVars[j]).getInitialExpression();
double defaultInitialCondition = 0;
try {
initial.bindExpression(mathDesc);
defaultInitialCondition = initial.evaluateConstant();
pw.println("\t" + defaultInitialCondition + ";\t\t% yinit(" + (j + 1) + ") is the initial condition for '" + cbit.vcell.parser.SymbolUtils.getEscapedTokenMatlab(volVars[j].getName()) + "'");
} catch (ExpressionException e) {
e.printStackTrace(System.out);
pw.println("\t" + initial.infix_Matlab() + ";\t\t% yinit(" + (j + 1) + ") is the initial condition for '" + cbit.vcell.parser.SymbolUtils.getEscapedTokenMatlab(volVars[j].getName()) + "'");
// throw new RuntimeException("error evaluating initial condition for variable "+volVars[j].getName());
}
}
pw.println("];");
pw.println("if nargin >= 2");
pw.println("\tif length(argYinit) > 0");
pw.println("\t\t%");
pw.println("\t\t% initial conditions overridden by function arguments");
pw.println("\t\t%");
pw.println("\t\tyinit = argYinit;");
pw.println("\tend");
pw.println("end");
pw.println("%");
pw.println("% Default Parameters");
pw.println("% constants are only those \"Constants\" from the Math Description that are just floating point numbers (no identifiers)");
pw.println("% note: constants of the form \"A_init\" are really initial conditions and are treated in \"yinit\"");
pw.println("%");
pw.println("param = [");
int paramIndex = 0;
for (int i = 0; i < constants.length; i++) {
pw.println("\t" + constants[i].getExpression().infix_Matlab() + ";\t\t% param(" + (paramIndex + 1) + ") is '" + cbit.vcell.parser.SymbolUtils.getEscapedTokenMatlab(constants[i].getName()) + "'");
paramIndex++;
}
pw.println("];");
pw.println("if nargin >= 3");
pw.println("\tif length(argParam) > 0");
pw.println("\t\t%");
pw.println("\t\t% parameter values overridden by function arguments");
pw.println("\t\t%");
pw.println("\t\tparam = argParam;");
pw.println("\tend");
pw.println("end");
pw.println("%");
pw.println("% invoke the integrator");
pw.println("%");
pw.println("[T,Y] = ode15s(@f,timeSpan,yinit,odeset('OutputFcn',@odeplot),param,yinit);");
pw.println("");
pw.println("% get the solution");
pw.println("all = zeros(size(T), numVars);");
pw.println("for i = 1:size(T)");
pw.println("\tall(i,:) = getRow(T(i), Y(i,:), yinit, param);");
pw.println("end");
pw.println("");
pw.println("allValues = all;");
pw.println("end");
// get row data for solution
pw.println("");
pw.println("% -------------------------------------------------------");
pw.println("% get row data");
pw.println("function rowValue = getRow(t,y,y0,p)");
//
// print volVariables (in order and assign to var vector)
//
pw.println("\t% State Variables");
for (int i = 0; i < volVars.length; i++) {
pw.println("\t" + cbit.vcell.parser.SymbolUtils.getEscapedTokenMatlab(volVars[i].getName()) + " = y(" + (i + 1) + ");");
}
//
// print constants
//
pw.println("\t% Constants");
paramIndex = 0;
for (int i = 0; i < constants.length; i++) {
pw.println("\t" + cbit.vcell.parser.SymbolUtils.getEscapedTokenMatlab(constants[i].getName()) + " = p(" + (paramIndex + 1) + ");");
paramIndex++;
}
//
// print variables
//
pw.println("\t% Functions");
for (int i = 0; i < functions.length; i++) {
pw.println("\t" + cbit.vcell.parser.SymbolUtils.getEscapedTokenMatlab(functions[i].getName()) + " = " + functions[i].getExpression().infix_Matlab() + ";");
}
pw.println("");
pw.println("\trowValue = [" + varNamesForValueArray + "];");
pw.println("end");
//
// print ode-rate
//
pw.println("");
pw.println("% -------------------------------------------------------");
pw.println("% ode rate");
pw.println("function dydt = f(t,y,p,y0)");
//
// print volVariables (in order and assign to var vector)
//
pw.println("\t% State Variables");
for (int i = 0; i < volVars.length; i++) {
pw.println("\t" + cbit.vcell.parser.SymbolUtils.getEscapedTokenMatlab(volVars[i].getName()) + " = y(" + (i + 1) + ");");
}
//
// print constants
//
pw.println("\t% Constants");
paramIndex = 0;
for (int i = 0; i < constants.length; i++) {
pw.println("\t" + cbit.vcell.parser.SymbolUtils.getEscapedTokenMatlab(constants[i].getName()) + " = p(" + (paramIndex + 1) + ");");
paramIndex++;
}
//
// print variables
//
pw.println("\t% Functions");
for (int i = 0; i < functions.length; i++) {
pw.println("\t" + cbit.vcell.parser.SymbolUtils.getEscapedTokenMatlab(functions[i].getName()) + " = " + functions[i].getExpression().infix_Matlab() + ";");
}
pw.println("\t% Rates");
pw.println("\tdydt = [");
for (int i = 0; i < volVars.length; i++) {
pw.println("\t\t" + subDomain.getEquation(volVars[i]).getRateExpression().infix_Matlab() + "; % rate for " + cbit.vcell.parser.SymbolUtils.getEscapedTokenMatlab(volVars[i].getName()));
}
pw.println("\t];");
pw.println("end");
}
Also used :
Variable(cbit.vcell.math.Variable)
VolVariable(cbit.vcell.math.VolVariable)
MathDescription(cbit.vcell.math.MathDescription)
VolVariable(cbit.vcell.math.VolVariable)
VariableHash(cbit.vcell.math.VariableHash)
Constant(cbit.vcell.math.Constant)
ExpressionException(cbit.vcell.parser.ExpressionException)
Function(cbit.vcell.math.Function)
Expression(cbit.vcell.parser.Expression)
MathException(cbit.vcell.math.MathException)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
Vector(java.util.Vector)
Example 13 with VariableHash
use of cbit.vcell.math.VariableHash in project vcell by virtualcell.
the class StochMathMapping_4_8 method refreshMathDescription.
/**
* set up a math description based on current simulationContext.
*/
private void refreshMathDescription() throws MappingException, MatrixException, MathException, ExpressionException, ModelException {
// use local variable instead of using getter all the time.
SimulationContext simContext = getSimulationContext();
// local structure mapping list
StructureMapping[] structureMappings = simContext.getGeometryContext().getStructureMappings();
// We have to check if all the reactions are able to tranform to stochastic jump processes before generating the math.
String stochChkMsg = simContext.getModel().isValidForStochApp();
if (!(stochChkMsg.equals(""))) {
throw new ModelException("Problem updating math description: " + simContext.getName() + "\n" + stochChkMsg);
}
// All sizes must be set for new ODE models and ratios must be set for old ones.
simContext.checkValidity();
//
// verify that all structures are mapped to subvolumes and all subvolumes are mapped to a structure
//
Structure[] structures = simContext.getGeometryContext().getModel().getStructures();
for (int i = 0; i < structures.length; i++) {
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(structures[i]);
if (sm == null || (sm instanceof FeatureMapping && getSubVolume(((FeatureMapping) sm)) == null)) {
throw new MappingException("model structure '" + structures[i].getName() + "' not mapped to a geometry subVolume");
}
if (sm != null && (sm instanceof MembraneMapping) && ((MembraneMapping) sm).getVolumeFractionParameter() != null) {
Expression volFractExp = ((MembraneMapping) sm).getVolumeFractionParameter().getExpression();
try {
if (volFractExp != null) {
double volFract = volFractExp.evaluateConstant();
if (volFract >= 1.0) {
throw new MappingException("model structure '" + (getSimulationContext().getModel().getStructureTopology().getInsideFeature(((MembraneMapping) sm).getMembrane()).getName() + "' has volume fraction >= 1.0"));
}
}
} catch (ExpressionException e) {
e.printStackTrace(System.out);
}
}
}
SubVolume[] subVolumes = simContext.getGeometryContext().getGeometry().getGeometrySpec().getSubVolumes();
for (int i = 0; i < subVolumes.length; i++) {
if (getStructures(subVolumes[i]) == null || getStructures(subVolumes[i]).length == 0) {
throw new MappingException("geometry subVolume '" + subVolumes[i].getName() + "' not mapped from a model structure");
}
}
//
// gather only those reactionSteps that are not "excluded"
//
ReactionSpec[] reactionSpecs = simContext.getReactionContext().getReactionSpecs();
Vector<ReactionStep> rsList = new Vector<ReactionStep>();
for (int i = 0; i < reactionSpecs.length; i++) {
if (reactionSpecs[i].isExcluded() == false) {
rsList.add(reactionSpecs[i].getReactionStep());
}
}
ReactionStep[] reactionSteps = new ReactionStep[rsList.size()];
rsList.copyInto(reactionSteps);
//
for (int i = 0; i < reactionSteps.length; i++) {
Kinetics.UnresolvedParameter[] unresolvedParameters = reactionSteps[i].getKinetics().getUnresolvedParameters();
if (unresolvedParameters != null && unresolvedParameters.length > 0) {
StringBuffer buffer = new StringBuffer();
for (int j = 0; j < unresolvedParameters.length; j++) {
if (j > 0) {
buffer.append(", ");
}
buffer.append(unresolvedParameters[j].getName());
}
throw new MappingException(reactionSteps[i].getDisplayType() + " '" + reactionSteps[i].getName() + "' contains unresolved identifier(s): " + buffer);
}
}
//
// create new MathDescription (based on simContext's previous MathDescription if possible)
//
MathDescription oldMathDesc = simContext.getMathDescription();
mathDesc = null;
if (oldMathDesc != null) {
if (oldMathDesc.getVersion() != null) {
mathDesc = new MathDescription(oldMathDesc.getVersion());
} else {
mathDesc = new MathDescription(oldMathDesc.getName());
}
} else {
mathDesc = new MathDescription(simContext.getName() + "_generated");
}
//
// temporarily place all variables in a hashtable (before binding) and discarding duplicates
//
VariableHash varHash = new VariableHash();
//
// conversion factors
//
Model model = simContext.getModel();
ModelUnitSystem modelUnitSystem = model.getUnitSystem();
varHash.addVariable(new Constant(getMathSymbol(model.getKMOLE(), null), getIdentifierSubstitutions(model.getKMOLE().getExpression(), model.getKMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getN_PMOLE(), null), getIdentifierSubstitutions(model.getN_PMOLE().getExpression(), model.getN_PMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT().getExpression(), model.getFARADAY_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT_NMOLE(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT_NMOLE().getExpression(), model.getFARADAY_CONSTANT_NMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getGAS_CONSTANT(), null), getIdentifierSubstitutions(model.getGAS_CONSTANT().getExpression(), model.getGAS_CONSTANT().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getTEMPERATURE(), null), getIdentifierSubstitutions(new Expression(simContext.getTemperatureKelvin()), model.getTEMPERATURE().getUnitDefinition(), null)));
Enumeration<SpeciesContextMapping> enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = enum1.nextElement();
if (scm.getVariable() instanceof StochVolVariable) {
varHash.addVariable(scm.getVariable());
}
}
//
// add rate term for all reactions
// add current source terms for each reaction step in a membrane
//
/*for (int i = 0; i < reactionSteps.length; i++){
boolean bAllReactionParticipantsFixed = true;
ReactionParticipant rp_Array[] = reactionSteps[i].getReactionParticipants();
for (int j = 0; j < rp_Array.length; j++) {
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(rp_Array[j].getSpeciesContext());
if (!(rp_Array[j] instanceof Catalyst) && !scs.isConstant()){
bAllReactionParticipantsFixed = false; // found at least one reactionParticipant that is not fixed and needs this rate
}
}
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(reactionSteps[i].getStructure());
}---don't think it's useful, isn't it?*/
// deals with model parameters
ModelParameter[] modelParameters = simContext.getModel().getModelParameters();
for (int j = 0; j < modelParameters.length; j++) {
Expression expr = getSubstitutedExpr(modelParameters[j].getExpression(), true, false);
expr = getIdentifierSubstitutions(expr, modelParameters[j].getUnitDefinition(), null);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], null), expr));
}
// added July 2009, ElectricalStimulusParameter electric mapping tab
ElectricalStimulus[] elecStimulus = simContext.getElectricalStimuli();
if (elecStimulus.length > 0) {
throw new MappingException("Modles with electrophysiology are not supported for stochastic applications.");
}
for (int j = 0; j < structureMappings.length; j++) {
if (structureMappings[j] instanceof MembraneMapping) {
MembraneMapping memMapping = (MembraneMapping) structureMappings[j];
Parameter initialVoltageParm = memMapping.getInitialVoltageParameter();
try {
Expression exp = initialVoltageParm.getExpression();
exp.evaluateConstant();
varHash.addVariable(newFunctionOrConstant(getMathSymbol(memMapping.getMembrane().getMembraneVoltage(), memMapping), getIdentifierSubstitutions(memMapping.getInitialVoltageParameter().getExpression(), memMapping.getInitialVoltageParameter().getUnitDefinition(), memMapping)));
} catch (ExpressionException e) {
e.printStackTrace(System.out);
throw new MappingException("Membrane initial voltage: " + initialVoltageParm.getName() + " cannot be evaluated as constant.");
}
}
}
//
for (int j = 0; j < reactionSteps.length; j++) {
ReactionStep rs = reactionSteps[j];
if (simContext.getReactionContext().getReactionSpec(rs).isExcluded()) {
continue;
}
if (rs.getKinetics() instanceof LumpedKinetics) {
throw new RuntimeException("Lumped Kinetics not yet supported for Stochastic Math Generation");
}
Kinetics.KineticsParameter[] parameters = rs.getKinetics().getKineticsParameters();
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(rs.getStructure());
if (parameters != null) {
for (int i = 0; i < parameters.length; i++) {
if ((parameters[i].getRole() == Kinetics.ROLE_CurrentDensity) && (parameters[i].getExpression() == null || parameters[i].getExpression().isZero())) {
continue;
}
// don't add rate, we'll do it later when creating the jump processes
if (parameters[i].getRole() != Kinetics.ROLE_ReactionRate) {
Expression expr = getSubstitutedExpr(parameters[i].getExpression(), true, false);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[i], sm), getIdentifierSubstitutions(expr, parameters[i].getUnitDefinition(), sm)));
}
}
}
}
// the parameter "Size" is already put into mathsymbolmapping in refreshSpeciesContextMapping()
for (int i = 0; i < structureMappings.length; i++) {
StructureMapping sm = structureMappings[i];
StructureMapping.StructureMappingParameter parm = sm.getParameterFromRole(StructureMapping.ROLE_Size);
if (parm.getExpression() != null) {
try {
double value = parm.getExpression().evaluateConstant();
varHash.addVariable(new Constant(getMathSymbol(parm, sm), new Expression(value)));
} catch (ExpressionException e) {
// varHash.addVariable(new Function(getMathSymbol0(parm,sm),getIdentifierSubstitutions(parm.getExpression(),parm.getUnitDefinition(),sm)));
e.printStackTrace(System.out);
throw new MappingException("Size of structure:" + sm.getNameScope().getName() + " cannot be evaluated as constant.");
}
}
}
//
// species initial values (either function or constant)
//
SpeciesContextSpec[] speciesContextSpecs = simContext.getReactionContext().getSpeciesContextSpecs();
for (int i = 0; i < speciesContextSpecs.length; i++) {
// can be concentration or amount
SpeciesContextSpec.SpeciesContextSpecParameter initParam = null;
Expression iniExp = null;
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
if (speciesContextSpecs[i].getInitialConcentrationParameter() != null && speciesContextSpecs[i].getInitialConcentrationParameter().getExpression() != null) {
// use concentration, need to set up amount functions
initParam = speciesContextSpecs[i].getInitialConcentrationParameter();
iniExp = initParam.getExpression();
iniExp = getSubstitutedExpr(iniExp, true, !speciesContextSpecs[i].isConstant());
// now create the appropriate function or Constant for the speciesContextSpec.
varHash.addVariable(newFunctionOrConstant(getMathSymbol(initParam, sm), getIdentifierSubstitutions(iniExp, initParam.getUnitDefinition(), sm)));
// add function for initial amount
SpeciesContextSpec.SpeciesContextSpecParameter initAmountParam = speciesContextSpecs[i].getInitialCountParameter();
Expression iniAmountExp = getExpressionConcToAmt(new Expression(initParam, getNameScope()), speciesContextSpecs[i].getSpeciesContext());
// iniAmountExp.bindExpression(this);
varHash.addVariable(new Function(getMathSymbol(initAmountParam, sm), getIdentifierSubstitutions(iniAmountExp, initAmountParam.getUnitDefinition(), sm), nullDomain));
} else if (speciesContextSpecs[i].getInitialCountParameter() != null && speciesContextSpecs[i].getInitialCountParameter().getExpression() != null) {
// use amount
initParam = speciesContextSpecs[i].getInitialCountParameter();
iniExp = initParam.getExpression();
iniExp = getSubstitutedExpr(iniExp, false, !speciesContextSpecs[i].isConstant());
// now create the appropriate function or Constant for the speciesContextSpec.
varHash.addVariable(newFunctionOrConstant(getMathSymbol(initParam, sm), getIdentifierSubstitutions(iniExp, initParam.getUnitDefinition(), sm)));
}
// add spConcentration (concentration of species) to varHash as function or constant
SpeciesConcentrationParameter spConcParam = getSpeciesConcentrationParameter(speciesContextSpecs[i].getSpeciesContext());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(spConcParam, sm), getIdentifierSubstitutions(spConcParam.getExpression(), spConcParam.getUnitDefinition(), sm)));
}
//
// constant species (either function or constant)
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() instanceof Constant) {
varHash.addVariable(scm.getVariable());
}
}
//
if (simContext.getGeometryContext().getGeometry() != null) {
try {
mathDesc.setGeometry(simContext.getGeometryContext().getGeometry());
} catch (java.beans.PropertyVetoException e) {
e.printStackTrace(System.out);
throw new MappingException("failure setting geometry " + e.getMessage());
}
} else {
throw new MappingException("geometry must be defined");
}
//
// functions: species which is not a variable, but has dependency expression
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() == null && scm.getDependencyExpression() != null) {
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(scm.getSpeciesContext().getStructure());
Expression exp = scm.getDependencyExpression();
exp.bindExpression(this);
SpeciesCountParameter spCountParam = getSpeciesCountParameter(scm.getSpeciesContext());
varHash.addVariable(new Function(getMathSymbol(spCountParam, sm), getIdentifierSubstitutions(exp, spCountParam.getUnitDefinition(), sm), nullDomain));
}
}
//
// create subDomains
//
SubDomain subDomain = null;
subVolumes = simContext.getGeometryContext().getGeometry().getGeometrySpec().getSubVolumes();
for (int j = 0; j < subVolumes.length; j++) {
SubVolume subVolume = (SubVolume) subVolumes[j];
//
// get priority of subDomain
//
int priority;
Feature spatialFeature = getResolvedFeature(subVolume);
if (spatialFeature == null) {
if (simContext.getGeometryContext().getGeometry().getDimension() > 0) {
throw new MappingException("no compartment (in Physiology) is mapped to subdomain '" + subVolume.getName() + "' (in Geometry)");
} else {
priority = CompartmentSubDomain.NON_SPATIAL_PRIORITY;
}
} else {
// now does not have to match spatial feature, *BUT* needs to be unique
priority = j;
}
subDomain = new CompartmentSubDomain(subVolume.getName(), priority);
mathDesc.addSubDomain(subDomain);
}
// ReactionSpec[] reactionSpecs = simContext.getReactionContext().getReactionSpecs();---need to take a look here!
for (int i = 0; i < reactionSpecs.length; i++) {
if (reactionSpecs[i].isExcluded()) {
continue;
}
// get the reaction
ReactionStep reactionStep = reactionSpecs[i].getReactionStep();
Kinetics kinetics = reactionStep.getKinetics();
// the structure where reaction happens
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(reactionStep.getStructure());
// create symbol table for jump process based on reactionStep and structure mapping
// final ReactionStep finalRS = reactionStep;
// final StructureMapping finalSM = sm;
// SymbolTable symTable = new SymbolTable(){
// public SymbolTableEntry getEntry(String identifierString) throws ExpressionBindingException {
// SymbolTableEntry ste = finalRS.getEntry(identifierString);
// if(ste == null)
// {
// ste = finalSM.getEntry(identifierString);
// }
// return ste;
// }
// };
// Different ways to deal with simple reactions and flux reactions
// probability parameter from modelUnitSystem
VCUnitDefinition probabilityParamUnit = modelUnitSystem.getStochasticSubstanceUnit().divideBy(modelUnitSystem.getTimeUnit());
if (// simple reactions
reactionStep instanceof SimpleReaction) {
// check the reaction rate law to see if we need to decompose a reaction(reversible) into two jump processes.
// rate constants are important in calculating the probability rate.
// for Mass Action, we use KForward and KReverse,
// for General Kinetics we parse reaction rate J to see if it is in Mass Action form.
Expression forwardRate = null;
Expression reverseRate = null;
if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
forwardRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KForward).getExpression();
reverseRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KReverse).getExpression();
} else if (kinetics.getKineticsDescription().equals(KineticsDescription.General)) {
Expression rateExp = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate).getExpression();
MassActionSolver.MassActionFunction maFunc = MassActionSolver.solveMassAction(null, null, rateExp, reactionStep);
if (maFunc.getForwardRate() == null && maFunc.getReverseRate() == null) {
throw new MappingException("Cannot generate stochastic math mapping for the reaction:" + reactionStep.getName() + "\nLooking for the rate function according to the form of k1*Reactant1^Stoir1*Reactant2^Stoir2...-k2*Product1^Stoip1*Product2^Stoip2.");
} else {
if (maFunc.getForwardRate() != null) {
forwardRate = maFunc.getForwardRate();
}
if (maFunc.getReverseRate() != null) {
reverseRate = maFunc.getReverseRate();
}
}
}
/*else if (kinetics.getKineticsDescription().getName().compareTo(KineticsDescription.HMM_irreversible.getName())==0)
{
forwardRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Km).getExpression();
}
else if (kinetics.getKineticsDescription().getName().compareTo(KineticsDescription.HMM_reversible.getName())==0)
{
forwardRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KmFwd).getExpression();
reverseRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KmRev).getExpression();
}*/
boolean isForwardRatePresent = false;
boolean isReverseRatePresent = false;
if (forwardRate != null) {
isForwardRatePresent = true;
}
if (reverseRate != null) {
isReverseRatePresent = true;
}
// we process it as forward reaction
if ((isForwardRatePresent)) /*|| ((forwardRate == null) && (reverseRate == null))*/
{
// get jump process name
String jpName = TokenMangler.mangleToSName(reactionStep.getName());
// get probability
Expression exp = null;
// reactions are mass actions
exp = getProbabilityRate(reactionStep, true);
// bind symbol table before substitute identifiers in the reaction step
exp.bindExpression(this);
MathMapping_4_8.ProbabilityParameter probParm = null;
try {
probParm = addProbabilityParameter("P_" + jpName, exp, MathMapping_4_8.PARAMETER_ROLE_P, probabilityParamUnit, reactionSpecs[i]);
} catch (PropertyVetoException pve) {
pve.printStackTrace();
throw new MappingException(pve.getMessage());
}
// add probability to function or constant
varHash.addVariable(newFunctionOrConstant(getMathSymbol(probParm, sm), getIdentifierSubstitutions(exp, probabilityParamUnit, sm)));
JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probParm, sm)));
// actions
ReactionParticipant[] reacPart = reactionStep.getReactionParticipants();
for (int j = 0; j < reacPart.length; j++) {
Action action = null;
SpeciesCountParameter spCountParam = getSpeciesCountParameter(reacPart[j].getSpeciesContext());
if (reacPart[j] instanceof Reactant) {
// check if the reactant is a constant. If the species is a constant, there will be no action taken on this species
if (// not a constant
!simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
int stoi = ((Reactant) reacPart[j]).getStoichiometry();
action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression("-" + String.valueOf(stoi)));
jp.addAction(action);
}
} else if (reacPart[j] instanceof Product) {
// check if the product is a constant. If the product is a constant, there will be no action taken on this species
if (// not a constant
!simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
int stoi = ((Product) reacPart[j]).getStoichiometry();
action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression(stoi));
jp.addAction(action);
}
}
}
// add jump process to compartment subDomain
subDomain.addJumpProcess(jp);
}
if (// one more jump process for a reversible reaction
isReverseRatePresent) {
// get jump process name
String jpName = TokenMangler.mangleToSName(reactionStep.getName()) + "_reverse";
Expression exp = null;
// reactions are mass actions
exp = getProbabilityRate(reactionStep, false);
// bind symbol table before substitute identifiers in the reaction step
exp.bindExpression(this);
MathMapping_4_8.ProbabilityParameter probRevParm = null;
try {
probRevParm = addProbabilityParameter("P_" + jpName, exp, MathMapping_4_8.PARAMETER_ROLE_P_reverse, probabilityParamUnit, reactionSpecs[i]);
} catch (PropertyVetoException pve) {
pve.printStackTrace();
throw new MappingException(pve.getMessage());
}
// add probability to function or constant
varHash.addVariable(newFunctionOrConstant(getMathSymbol(probRevParm, sm), getIdentifierSubstitutions(exp, probabilityParamUnit, sm)));
JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probRevParm, sm)));
// actions
ReactionParticipant[] reacPart = reactionStep.getReactionParticipants();
for (int j = 0; j < reacPart.length; j++) {
Action action = null;
SpeciesCountParameter spCountParam = getSpeciesCountParameter(reacPart[j].getSpeciesContext());
if (reacPart[j] instanceof Reactant) {
// check if the reactant is a constant. If the species is a constant, there will be no action taken on this species
if (// not a constant
!simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
int stoi = ((Reactant) reacPart[j]).getStoichiometry();
action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression(stoi));
jp.addAction(action);
}
} else if (reacPart[j] instanceof Product) {
// check if the product is a constant. If the product is a constant, there will be no action taken on this species
if (// not a constant
!simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
int stoi = ((Product) reacPart[j]).getStoichiometry();
action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression("-" + String.valueOf(stoi)));
jp.addAction(action);
}
}
}
// add jump process to compartment subDomain
subDomain.addJumpProcess(jp);
}
// end of if(isForwardRateNonZero), if(isReverseRateNonRate)
} else if (// flux reactions
reactionStep instanceof FluxReaction) {
// we could set jump processes for general flux rate in forms of p1*Sout + p2*Sin
if (kinetics.getKineticsDescription().equals(KineticsDescription.General)) {
Expression fluxRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate).getExpression();
// we have to pass the math description para to flux solver, coz somehow math description in simulation context is not updated.
MassActionSolver.MassActionFunction fluxFunc = MassActionSolver.solveMassAction(null, null, fluxRate, (FluxReaction) reactionStep);
// create jump process for forward flux if it exists.
if (fluxFunc.getForwardRate() != null && !fluxFunc.getForwardRate().isZero()) {
// jump process name
// +"_reverse";
String jpName = TokenMangler.mangleToSName(reactionStep.getName());
// we do it here instead of fluxsolver, coz we need to use getMathSymbol0(), structuremapping...etc.
Expression rate = fluxFunc.getForwardRate();
// get species expression (depend on structure, if mem: Species/mem_Size, if vol: species*KMOLE/vol_size)
SpeciesContext scOut = fluxFunc.getReactants().get(0).getSpeciesContext();
Expression speciesFactor = null;
if (scOut.getStructure() instanceof Feature) {
Expression exp1 = new Expression(1.0 / 602.0);
Expression exp2 = new Expression(scOut.getStructure().getStructureSize(), getNameScope());
speciesFactor = Expression.div(Expression.invert(exp1), exp2);
} else {
throw new MappingException("Species involved in a flux have to be volume species.");
}
Expression speciesExp = Expression.mult(speciesFactor, new Expression(scOut, getNameScope()));
// get probability expression by adding factor to rate (rate: rate*size_mem/KMOLE)
Expression expr1 = Expression.mult(rate, speciesExp);
Expression numeratorExpr = Expression.mult(expr1, new Expression(sm.getStructure().getStructureSize(), getNameScope()));
Expression exp = new Expression(1.0 / 602.0);
Expression probExp = Expression.mult(numeratorExpr, exp);
// bind symbol table before substitute identifiers in the reaction step
probExp.bindExpression(reactionStep);
MathMapping_4_8.ProbabilityParameter probParm = null;
try {
probParm = addProbabilityParameter("P_" + jpName, probExp, MathMapping_4_8.PARAMETER_ROLE_P, probabilityParamUnit, reactionSpecs[i]);
} catch (PropertyVetoException pve) {
pve.printStackTrace();
throw new MappingException(pve.getMessage());
}
// add probability to function or constant
varHash.addVariable(newFunctionOrConstant(getMathSymbol(probParm, sm), getIdentifierSubstitutions(probExp, probabilityParamUnit, sm)));
JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probParm, sm)));
// actions
Action action = null;
SpeciesContext sc = fluxFunc.getReactants().get(0).getSpeciesContext();
if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression(-1));
jp.addAction(action);
}
sc = fluxFunc.getProducts().get(0).getSpeciesContext();
if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression(1));
jp.addAction(action);
}
subDomain.addJumpProcess(jp);
}
if (fluxFunc.getReverseRate() != null && !fluxFunc.getReverseRate().isZero()) {
// jump process name
String jpName = TokenMangler.mangleToSName(reactionStep.getName()) + "_reverse";
Expression rate = fluxFunc.getReverseRate();
// get species expression (depend on structure, if mem: Species/mem_Size, if vol: species*KMOLE/vol_size)
SpeciesContext scIn = fluxFunc.getProducts().get(0).getSpeciesContext();
Expression speciesFactor = null;
if (scIn.getStructure() instanceof Feature) {
Expression exp1 = new Expression(1.0 / 602.0);
Expression exp2 = new Expression(scIn.getStructure().getStructureSize(), getNameScope());
speciesFactor = Expression.div(Expression.invert(exp1), exp2);
} else {
throw new MappingException("Species involved in a flux have to be volume species.");
}
Expression speciesExp = Expression.mult(speciesFactor, new Expression(scIn, getNameScope()));
// get probability expression by adding factor to rate (rate: rate*size_mem/KMOLE)
Expression expr1 = Expression.mult(rate, speciesExp);
Expression numeratorExpr = Expression.mult(expr1, new Expression(sm.getStructure().getStructureSize(), getNameScope()));
Expression exp = new Expression(1.0 / 602.0);
Expression probRevExp = Expression.mult(numeratorExpr, exp);
// bind symbol table before substitute identifiers in the reaction step
probRevExp.bindExpression(reactionStep);
MathMapping_4_8.ProbabilityParameter probRevParm = null;
try {
probRevParm = addProbabilityParameter("P_" + jpName, probRevExp, MathMapping_4_8.PARAMETER_ROLE_P_reverse, probabilityParamUnit, reactionSpecs[i]);
} catch (PropertyVetoException pve) {
pve.printStackTrace();
throw new MappingException(pve.getMessage());
}
// add probability to function or constant
varHash.addVariable(newFunctionOrConstant(getMathSymbol(probRevParm, sm), getIdentifierSubstitutions(probRevExp, probabilityParamUnit, sm)));
JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probRevParm, sm)));
// actions
Action action = null;
SpeciesContext sc = fluxFunc.getReactants().get(0).getSpeciesContext();
if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression(1));
jp.addAction(action);
}
sc = fluxFunc.getProducts().get(0).getSpeciesContext();
if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression(-1));
jp.addAction(action);
}
subDomain.addJumpProcess(jp);
}
}
}
// end of if (simplereaction)...else if(fluxreaction)
}
// end of reaction step loop
//
// set Variables to MathDescription all at once with the order resolved by "VariableHash"
//
mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
// set up variable initial conditions in subDomain
SpeciesContextSpec[] scSpecs = simContext.getReactionContext().getSpeciesContextSpecs();
for (int i = 0; i < speciesContextSpecs.length; i++) {
// get stochastic variable by name
SpeciesCountParameter spCountParam = getSpeciesCountParameter(speciesContextSpecs[i].getSpeciesContext());
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
String varName = getMathSymbol(spCountParam, sm);
if (scSpecs[i].isConstant()) {
continue;
}
StochVolVariable var = (StochVolVariable) mathDesc.getVariable(varName);
// stochastic use initial number of particles
SpeciesContextSpec.SpeciesContextSpecParameter initParm = scSpecs[i].getInitialCountParameter();
// stochastic variables initial expression.
if (initParm != null) {
VarIniCondition varIni = new VarIniCount(var, new Expression(getMathSymbol(initParm, sm)));
subDomain.addVarIniCondition(varIni);
}
}
if (!mathDesc.isValid()) {
throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
}
}
Also used :
VarIniCondition(cbit.vcell.math.VarIniCondition)
MembraneMapping(cbit.vcell.mapping.MembraneMapping)
LumpedKinetics(cbit.vcell.model.LumpedKinetics)
MathDescription(cbit.vcell.math.MathDescription)
SpeciesContextMapping(cbit.vcell.mapping.SpeciesContextMapping)
Product(cbit.vcell.model.Product)
FluxReaction(cbit.vcell.model.FluxReaction)
SpeciesContext(cbit.vcell.model.SpeciesContext)
SpeciesContextSpec(cbit.vcell.mapping.SpeciesContextSpec)
Feature(cbit.vcell.model.Feature)
Reactant(cbit.vcell.model.Reactant)
ExpressionException(cbit.vcell.parser.ExpressionException)
MappingException(cbit.vcell.mapping.MappingException)
PropertyVetoException(java.beans.PropertyVetoException)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
SubVolume(cbit.vcell.geometry.SubVolume)
Vector(java.util.Vector)
ModelException(cbit.vcell.model.ModelException)
ReactionSpec(cbit.vcell.mapping.ReactionSpec)
PropertyVetoException(java.beans.PropertyVetoException)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
ReactionStep(cbit.vcell.model.ReactionStep)
Kinetics(cbit.vcell.model.Kinetics)
LumpedKinetics(cbit.vcell.model.LumpedKinetics)
ReactionParticipant(cbit.vcell.model.ReactionParticipant)
Action(cbit.vcell.math.Action)
VariableHash(cbit.vcell.math.VariableHash)
Constant(cbit.vcell.math.Constant)
StructureMapping(cbit.vcell.mapping.StructureMapping)
Function(cbit.vcell.math.Function)
FeatureMapping(cbit.vcell.mapping.FeatureMapping)
JumpProcess(cbit.vcell.math.JumpProcess)
Structure(cbit.vcell.model.Structure)
StochVolVariable(cbit.vcell.math.StochVolVariable)
ModelUnitSystem(cbit.vcell.model.ModelUnitSystem)
SimpleReaction(cbit.vcell.model.SimpleReaction)
VarIniCount(cbit.vcell.math.VarIniCount)
SimulationContext(cbit.vcell.mapping.SimulationContext)
ElectricalStimulus(cbit.vcell.mapping.ElectricalStimulus)
Expression(cbit.vcell.parser.Expression)
Model(cbit.vcell.model.Model)
ProxyParameter(cbit.vcell.model.ProxyParameter)
Parameter(cbit.vcell.model.Parameter)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
Example 14 with VariableHash
use of cbit.vcell.math.VariableHash in project vcell by virtualcell.
the class XmlReader method getMathDescription.
/**
* This method returns a MathDescription from a XML element.
* Creation date: (4/26/2001 12:11:14 PM)
* @return cbit.vcell.math.MathDescription
* @param param org.jdom.Element
* @exception cbit.vcell.xml.XmlParseException The exception description.
*/
MathDescription getMathDescription(Element param, Geometry geometry) throws XmlParseException {
MathDescription mathdes = null;
Element tempelement;
// Retrieve Metadata(Version)
Version version = getVersion(param.getChild(XMLTags.VersionTag, vcNamespace));
// Retrieve attributes
String name = unMangle(param.getAttributeValue(XMLTags.NameAttrTag));
// Create new MathDescription
if (version != null) {
mathdes = new MathDescription(version);
} else {
mathdes = new MathDescription(name);
}
try {
// this step is needed!
mathdes.setGeometry(geometry);
} catch (java.beans.PropertyVetoException e) {
e.printStackTrace();
throw new XmlParseException("a PropertyVetoException was fired when setting the Geometry to the Mathdescription in the simContext " + name, e);
}
// set attributes
try {
mathdes.setName(name);
// String annotation = param.getAttributeValue(XMLTags.AnnotationAttrTag);
// if (annotation!=null) {
// mathdes.setDescription(unMangle(annotation));
// }
// add Annotation
String annotationText = param.getChildText(XMLTags.AnnotationTag, vcNamespace);
if (annotationText != null && annotationText.length() > 0) {
mathdes.setDescription(unMangle(annotationText));
}
} catch (java.beans.PropertyVetoException e) {
e.printStackTrace();
throw new XmlParseException("A PropertyVetoException was fired when setting the name " + name + ", to a new MathDescription!", e);
}
VariableHash varHash = new VariableHash();
// Retrieve Constant
Iterator<Element> iterator = param.getChildren(XMLTags.ConstantTag, vcNamespace).iterator();
while (iterator.hasNext()) {
tempelement = (Element) iterator.next();
try {
varHash.addVariable(getConstant(tempelement));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException(e);
}
}
// Retrieve FilamentRegionVariables
iterator = param.getChildren(XMLTags.FilamentRegionVariableTag, vcNamespace).iterator();
while (iterator.hasNext()) {
tempelement = (Element) iterator.next();
try {
varHash.addVariable(getFilamentRegionVariable(tempelement));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException(e);
}
}
// Retrieve FilamentVariables
iterator = param.getChildren(XMLTags.FilamentVariableTag, vcNamespace).iterator();
while (iterator.hasNext()) {
tempelement = (Element) iterator.next();
try {
varHash.addVariable(getFilamentVariable(tempelement));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException(e);
}
}
// retrieve InsideVariables
// **** This variables are for internal USE ******
// Retrieve MembraneRegionVariable
iterator = param.getChildren(XMLTags.MembraneRegionVariableTag, vcNamespace).iterator();
while (iterator.hasNext()) {
tempelement = (Element) iterator.next();
try {
varHash.addVariable(getMembraneRegionVariable(tempelement));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException(e);
}
}
// Retrieve MembraneVariable
iterator = param.getChildren(XMLTags.MembraneVariableTag, vcNamespace).iterator();
while (iterator.hasNext()) {
tempelement = (Element) iterator.next();
try {
varHash.addVariable(getMemVariable(tempelement));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException(e);
}
}
// Retrieve PointVariable
iterator = param.getChildren(XMLTags.PointVariableTag, vcNamespace).iterator();
while (iterator.hasNext()) {
tempelement = (Element) iterator.next();
try {
varHash.addVariable(getPointVariable(tempelement));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException(e);
}
}
// retrieve OutsideVariables
// **** This variables are for internal USE ******
// Retrieve Volume Region variable
iterator = param.getChildren(XMLTags.VolumeRegionVariableTag, vcNamespace).iterator();
while (iterator.hasNext()) {
tempelement = (Element) iterator.next();
try {
varHash.addVariable(getVolumeRegionVariable(tempelement));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException(e);
}
}
// Retrieve VolumeVariable
iterator = param.getChildren(XMLTags.VolumeVariableTag, vcNamespace).iterator();
while (iterator.hasNext()) {
tempelement = (Element) iterator.next();
try {
varHash.addVariable(getVolVariable(tempelement));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException(e);
}
}
// Retrieve StochVolVariable
iterator = param.getChildren(XMLTags.StochVolVariableTag, vcNamespace).iterator();
while (iterator.hasNext()) {
tempelement = (Element) iterator.next();
try {
varHash.addVariable(getStochVolVariable(tempelement));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException(e);
}
}
// Retrieve all the Functions //This needs to be processed before all the variables are read!
iterator = param.getChildren(XMLTags.FunctionTag, vcNamespace).iterator();
while (iterator.hasNext()) {
tempelement = (Element) iterator.next();
try {
varHash.addVariable(getFunction(tempelement));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException(e);
}
}
iterator = param.getChildren(XMLTags.VolumeRandomVariableTag, vcNamespace).iterator();
while (iterator.hasNext()) {
tempelement = (Element) iterator.next();
try {
varHash.addVariable(getRandomVariable(tempelement));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException(e);
}
}
iterator = param.getChildren(XMLTags.MembraneRandomVariableTag, vcNamespace).iterator();
while (iterator.hasNext()) {
tempelement = (Element) iterator.next();
try {
varHash.addVariable(getRandomVariable(tempelement));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException(e);
}
}
iterator = param.getChildren(XMLTags.VolumeParticleVariableTag, vcNamespace).iterator();
while (iterator.hasNext()) {
tempelement = (Element) iterator.next();
try {
varHash.addVariable(getVolumeParticalVariable(tempelement));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException(e);
}
}
iterator = param.getChildren(XMLTags.MembraneParticleVariableTag, vcNamespace).iterator();
while (iterator.hasNext()) {
tempelement = (Element) iterator.next();
try {
varHash.addVariable(getMembraneParticalVariable(tempelement));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException(e);
}
}
// ParticleMolecularTypeTag getParticleMolecularTypes
// has to be done before VolumeParticleSpeciesPattern and VolumeParticleObservable
iterator = param.getChildren(XMLTags.ParticleMolecularTypeTag, vcNamespace).iterator();
while (iterator.hasNext()) {
tempelement = (Element) iterator.next();
mathdes.addParticleMolecularType(getParticleMolecularType(tempelement));
}
// VolumeParticleSpeciesPatternTag
iterator = param.getChildren(XMLTags.VolumeParticleSpeciesPatternTag, vcNamespace).iterator();
while (iterator.hasNext()) {
tempelement = (Element) iterator.next();
try {
varHash.addVariable(getVolumeParticleSpeciesPattern(tempelement, mathdes));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException(e);
}
}
// VolumeParticleObservableTag getParticleObservables
iterator = param.getChildren(XMLTags.VolumeParticleObservableTag, vcNamespace).iterator();
while (iterator.hasNext()) {
tempelement = (Element) iterator.next();
try {
varHash.addVariable(getVolumeParticleObservable(tempelement, varHash));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException(e);
}
}
//
try {
mathdes.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException("Error adding the Function variables to the MathDescription " + name, e);
} catch (ExpressionBindingException e) {
e.printStackTrace();
throw new XmlParseException("Error adding the Function variables to the MathDescription " + name, e);
}
// Retrieve CompartmentsSubdomains
iterator = param.getChildren(XMLTags.CompartmentSubDomainTag, vcNamespace).iterator();
while (iterator.hasNext()) {
tempelement = (Element) iterator.next();
try {
mathdes.addSubDomain(getCompartmentSubDomain(tempelement, mathdes));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException("Error adding a new CompartmentSubDomain to the MathDescription " + name, e);
}
}
// Retrieve MembraneSubdomains
iterator = param.getChildren(XMLTags.MembraneSubDomainTag, vcNamespace).iterator();
while (iterator.hasNext()) {
tempelement = (Element) iterator.next();
try {
mathdes.addSubDomain(getMembraneSubDomain(tempelement, mathdes));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException("Error adding a new MembraneSubDomain to the MathDescription " + name, e);
}
}
// Retrieve the FilamentSubdomain (if any)
tempelement = param.getChild(XMLTags.FilamentSubDomainTag, vcNamespace);
if (tempelement != null) {
try {
mathdes.addSubDomain(getFilamentSubDomain(tempelement, mathdes));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException("Error adding a new FilamentSubDomain to the MathDescription " + name, e);
}
}
// Retrieve the PointSubdomain (if any)
tempelement = param.getChild(XMLTags.PointSubDomainTag, vcNamespace);
if (tempelement != null) {
try {
mathdes.addSubDomain(getPointSubDomain(tempelement, mathdes));
} catch (MathException e) {
e.printStackTrace();
throw new XmlParseException("Error adding a new PointSubDomain to the MathDescription " + name, e);
}
}
iterator = param.getChildren(XMLTags.EventTag, vcNamespace).iterator();
while (iterator.hasNext()) {
tempelement = (Element) iterator.next();
Event event = getEvent(mathdes, tempelement);
try {
mathdes.addEvent(event);
} catch (MathException e) {
e.printStackTrace(System.out);
throw new XmlParseException(e);
}
}
iterator = param.getChildren(XMLTags.PostProcessingBlock, vcNamespace).iterator();
while (iterator.hasNext()) {
tempelement = (Element) iterator.next();
getPostProcessingBlock(mathdes, tempelement);
}
return mathdes;
}
Also used :
PropertyVetoException(java.beans.PropertyVetoException)
MathDescription(cbit.vcell.math.MathDescription)
Version(org.vcell.util.document.Version)
RedistributionVersion(cbit.vcell.solvers.mb.MovingBoundarySolverOptions.RedistributionVersion)
SimulationVersion(org.vcell.util.document.SimulationVersion)
VCellSoftwareVersion(org.vcell.util.document.VCellSoftwareVersion)
VariableHash(cbit.vcell.math.VariableHash)
MathException(cbit.vcell.math.MathException)
Element(org.jdom.Element)
Event(cbit.vcell.math.Event)
BioEvent(cbit.vcell.mapping.BioEvent)
ExpressionBindingException(cbit.vcell.parser.ExpressionBindingException)
Example 15 with VariableHash
use of cbit.vcell.math.VariableHash in project vcell by virtualcell.
the class XmlReader method getKinetics.
/**
* This method returns a Kinetics object from a XML Element based on the value of the kinetics type attribute.
* Creation date: (3/19/2001 4:42:04 PM)
* @return cbit.vcell.model.Kinetics
* @param param org.jdom.Element
*/
private Kinetics getKinetics(Element param, ReactionStep reaction, Model model) throws XmlParseException {
VariableHash varHash = new VariableHash();
addResevedSymbols(varHash, model);
String type = param.getAttributeValue(XMLTags.KineticsTypeAttrTag);
Kinetics newKinetics = null;
try {
if (type.equalsIgnoreCase(XMLTags.KineticsTypeGeneralKinetics)) {
// create a general kinetics
newKinetics = new GeneralKinetics(reaction);
} else if (type.equalsIgnoreCase(XMLTags.KineticsTypeGeneralCurrentKinetics)) {
// Create GeneralCurrentKinetics
newKinetics = new GeneralCurrentKinetics(reaction);
} else if (type.equalsIgnoreCase(XMLTags.KineticsTypeMassAction) && reaction instanceof SimpleReaction) {
// create a Mass Action kinetics
newKinetics = new MassActionKinetics((SimpleReaction) reaction);
} else if (type.equalsIgnoreCase(XMLTags.KineticsTypeNernst) && reaction instanceof FluxReaction) {
// create NernstKinetics
newKinetics = new NernstKinetics((FluxReaction) reaction);
} else if (type.equalsIgnoreCase(XMLTags.KineticsTypeGHK) && reaction instanceof FluxReaction) {
// create GHKKinetics
newKinetics = new GHKKinetics((FluxReaction) reaction);
} else if (type.equalsIgnoreCase(XMLTags.KineticsTypeHMM_Irr) && reaction instanceof SimpleReaction) {
// create HMM_IrrKinetics
newKinetics = new HMM_IRRKinetics((SimpleReaction) reaction);
} else if (type.equalsIgnoreCase(XMLTags.KineticsTypeHMM_Rev) && reaction instanceof SimpleReaction) {
// create HMM_RevKinetics
newKinetics = new HMM_REVKinetics((SimpleReaction) reaction);
} else if (type.equalsIgnoreCase(XMLTags.KineticsTypeGeneralTotal_oldname)) {
// create GeneralTotalKinetics
newKinetics = new GeneralLumpedKinetics(reaction);
} else if (type.equalsIgnoreCase(XMLTags.KineticsTypeGeneralLumped)) {
// create GeneralLumpedKinetics
newKinetics = new GeneralLumpedKinetics(reaction);
} else if (type.equalsIgnoreCase(XMLTags.KineticsTypeGeneralCurrentLumped)) {
// create GeneralCurrentLumpedKinetics
newKinetics = new GeneralCurrentLumpedKinetics(reaction);
} else if (type.equalsIgnoreCase(XMLTags.KineticsTypeGeneralPermeability) && reaction instanceof FluxReaction) {
// create GeneralPermeabilityKinetics
newKinetics = new GeneralPermeabilityKinetics((FluxReaction) reaction);
} else if (type.equalsIgnoreCase(XMLTags.KineticsTypeMacroscopic_Irr) && reaction instanceof SimpleReaction) {
// create Macroscopic_IRRKinetics
newKinetics = new Macroscopic_IRRKinetics((SimpleReaction) reaction);
} else if (type.equalsIgnoreCase(XMLTags.KineticsTypeMicroscopic_Irr) && reaction instanceof SimpleReaction) {
// create Microscopic_IRRKinetics
newKinetics = new Microscopic_IRRKinetics((SimpleReaction) reaction);
} else {
throw new XmlParseException("Unknown kinetics type: " + type);
}
} catch (ExpressionException e) {
e.printStackTrace();
throw new XmlParseException("Error creating the kinetics for reaction: " + reaction.getName(), e);
}
try {
// transaction begin flag ... yeah, this is a hack
newKinetics.reading(true);
// Read all of the parameters
List<Element> list = param.getChildren(XMLTags.ParameterTag, vcNamespace);
// add constants that may be used in kinetics.
// VariableHash varHash = getVariablesHash();
ArrayList<String> reserved = new ArrayList<String>();
ReservedSymbol[] reservedSymbols = reaction.getModel().getReservedSymbols();
for (ReservedSymbol rs : reservedSymbols) {
reserved.add(rs.getName());
}
try {
if (reaction.getStructure() instanceof Membrane) {
Membrane membrane = (Membrane) reaction.getStructure();
varHash.addVariable(new Constant(membrane.getMembraneVoltage().getName(), new Expression(0.0)));
reserved.add(membrane.getMembraneVoltage().getName());
}
//
// add Reactants, Products, and Catalysts (ReactionParticipants)
//
ReactionParticipant[] rp = reaction.getReactionParticipants();
for (int i = 0; i < rp.length; i++) {
varHash.addVariable(new Constant(rp[i].getName(), new Expression(0.0)));
}
} catch (MathException e) {
e.printStackTrace(System.out);
throw new XmlParseException("error reordering parameters according to dependencies: ", e);
}
//
for (Element xmlParam : list) {
String paramName = unMangle(xmlParam.getAttributeValue(XMLTags.NameAttrTag));
String role = xmlParam.getAttributeValue(XMLTags.ParamRoleAttrTag);
String paramExpStr = xmlParam.getText();
Expression paramExp = unMangleExpression(paramExpStr);
try {
if (varHash.getVariable(paramName) == null) {
varHash.addVariable(new Function(paramName, paramExp, null));
} else {
if (reserved.contains(paramName)) {
varHash.removeVariable(paramName);
varHash.addVariable(new Function(paramName, paramExp, null));
}
}
} catch (MathException e) {
e.printStackTrace(System.out);
throw new XmlParseException("error reordering parameters according to dependencies: ", e);
}
Kinetics.KineticsParameter tempParam = null;
if (!role.equals(XMLTags.ParamRoleUserDefinedTag)) {
tempParam = newKinetics.getKineticsParameterFromRole(Kinetics.getParamRoleFromDefaultDesc(role));
} else {
continue;
}
// hack for bringing in General Total kinetics without breaking.
if (tempParam == null && newKinetics instanceof GeneralLumpedKinetics) {
if (role.equals(Kinetics.GTK_AssumedCompartmentSize_oldname) || role.equals(Kinetics.GTK_ReactionRate_oldname) || role.equals(Kinetics.GTK_CurrentDensity_oldname)) {
continue;
} else if (role.equals(VCMODL.TotalRate_oldname)) {
tempParam = newKinetics.getKineticsParameterFromRole(Kinetics.ROLE_LumpedReactionRate);
}
}
// hack from bringing in chargeValence parameters without breaking
if (tempParam == null && Kinetics.getParamRoleFromDefaultDesc(role) == Kinetics.ROLE_ChargeValence) {
tempParam = newKinetics.getChargeValenceParameter();
}
if (tempParam == null) {
throw new XmlParseException("parameter with role '" + role + "' not found in kinetics type '" + type + "'");
}
//
if (!tempParam.getName().equals(paramName)) {
Kinetics.KineticsParameter multNameParam = newKinetics.getKineticsParameter(paramName);
int n = 0;
while (multNameParam != null) {
String tempName = paramName + "_" + n++;
newKinetics.renameParameter(paramName, tempName);
multNameParam = newKinetics.getKineticsParameter(tempName);
}
newKinetics.renameParameter(tempParam.getName(), paramName);
}
}
//
// create unresolved parameters for all unresolved symbols
//
String unresolvedSymbol = varHash.getFirstUnresolvedSymbol();
while (unresolvedSymbol != null) {
try {
// will turn into an UnresolvedParameter.
varHash.addVariable(new Function(unresolvedSymbol, new Expression(0.0), null));
} catch (MathException e) {
e.printStackTrace(System.out);
throw new XmlParseException(e);
}
newKinetics.addUnresolvedParameter(unresolvedSymbol);
unresolvedSymbol = varHash.getFirstUnresolvedSymbol();
}
Variable[] sortedVariables = varHash.getTopologicallyReorderedVariables();
ModelUnitSystem modelUnitSystem = reaction.getModel().getUnitSystem();
for (int i = sortedVariables.length - 1; i >= 0; i--) {
if (sortedVariables[i] instanceof Function) {
Function paramFunction = (Function) sortedVariables[i];
Element xmlParam = null;
for (int j = 0; j < list.size(); j++) {
Element tempParam = (Element) list.get(j);
if (paramFunction.getName().equals(unMangle(tempParam.getAttributeValue(XMLTags.NameAttrTag)))) {
xmlParam = tempParam;
break;
}
}
if (xmlParam == null) {
// must have been an unresolved parameter
continue;
}
String symbol = xmlParam.getAttributeValue(XMLTags.VCUnitDefinitionAttrTag);
VCUnitDefinition unit = null;
if (symbol != null) {
unit = modelUnitSystem.getInstance(symbol);
}
Kinetics.KineticsParameter tempParam = newKinetics.getKineticsParameter(paramFunction.getName());
if (tempParam == null) {
newKinetics.addUserDefinedKineticsParameter(paramFunction.getName(), paramFunction.getExpression(), unit);
} else {
newKinetics.setParameterValue(tempParam, paramFunction.getExpression());
tempParam.setUnitDefinition(unit);
}
}
}
} catch (java.beans.PropertyVetoException e) {
e.printStackTrace(System.out);
throw new XmlParseException("Exception while setting parameters for Reaction : " + reaction.getName(), e);
} catch (ExpressionException e) {
e.printStackTrace(System.out);
throw new XmlParseException("Exception while settings parameters for Reaction : " + reaction.getName(), e);
} finally {
newKinetics.reading(false);
}
return newKinetics;
}
Also used :
FilamentVariable(cbit.vcell.math.FilamentVariable)
OutsideVariable(cbit.vcell.math.OutsideVariable)
StochVolVariable(cbit.vcell.math.StochVolVariable)
RandomVariable(cbit.vcell.math.RandomVariable)
VolumeRandomVariable(cbit.vcell.math.VolumeRandomVariable)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
VolumeRegionVariable(cbit.vcell.math.VolumeRegionVariable)
InsideVariable(cbit.vcell.math.InsideVariable)
VolVariable(cbit.vcell.math.VolVariable)
MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable)
PointVariable(cbit.vcell.math.PointVariable)
MembraneRandomVariable(cbit.vcell.math.MembraneRandomVariable)
MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable)
ParticleVariable(cbit.vcell.math.ParticleVariable)
MemVariable(cbit.vcell.math.MemVariable)
FilamentRegionVariable(cbit.vcell.math.FilamentRegionVariable)
Variable(cbit.vcell.math.Variable)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
VariableHash(cbit.vcell.math.VariableHash)
HMM_IRRKinetics(cbit.vcell.model.HMM_IRRKinetics)
MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant)
Constant(cbit.vcell.math.Constant)
Element(org.jdom.Element)
ReservedSymbol(cbit.vcell.model.Model.ReservedSymbol)
NernstKinetics(cbit.vcell.model.NernstKinetics)
ArrayList(java.util.ArrayList)
FluxReaction(cbit.vcell.model.FluxReaction)
GeneralKinetics(cbit.vcell.model.GeneralKinetics)
GeneralLumpedKinetics(cbit.vcell.model.GeneralLumpedKinetics)
ExpressionException(cbit.vcell.parser.ExpressionException)
PropertyVetoException(java.beans.PropertyVetoException)
GeneralCurrentKinetics(cbit.vcell.model.GeneralCurrentKinetics)
AnnotatedFunction(cbit.vcell.solver.AnnotatedFunction)
Function(cbit.vcell.math.Function)
GeneralCurrentLumpedKinetics(cbit.vcell.model.GeneralCurrentLumpedKinetics)
Membrane(cbit.vcell.model.Membrane)
Macroscopic_IRRKinetics(cbit.vcell.model.Macroscopic_IRRKinetics)
GeneralPermeabilityKinetics(cbit.vcell.model.GeneralPermeabilityKinetics)
ModelUnitSystem(cbit.vcell.model.ModelUnitSystem)
SimpleReaction(cbit.vcell.model.SimpleReaction)
GHKKinetics(cbit.vcell.model.GHKKinetics)
HMM_REVKinetics(cbit.vcell.model.HMM_REVKinetics)
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
Expression(cbit.vcell.parser.Expression)
MathException(cbit.vcell.math.MathException)
MassActionKinetics(cbit.vcell.model.MassActionKinetics)
Macroscopic_IRRKinetics(cbit.vcell.model.Macroscopic_IRRKinetics)
Kinetics(cbit.vcell.model.Kinetics)
NernstKinetics(cbit.vcell.model.NernstKinetics)
GeneralKinetics(cbit.vcell.model.GeneralKinetics)
GeneralLumpedKinetics(cbit.vcell.model.GeneralLumpedKinetics)
GeneralPermeabilityKinetics(cbit.vcell.model.GeneralPermeabilityKinetics)
GHKKinetics(cbit.vcell.model.GHKKinetics)
MassActionKinetics(cbit.vcell.model.MassActionKinetics)
Microscopic_IRRKinetics(cbit.vcell.model.Microscopic_IRRKinetics)
GeneralCurrentKinetics(cbit.vcell.model.GeneralCurrentKinetics)
HMM_REVKinetics(cbit.vcell.model.HMM_REVKinetics)
HMM_IRRKinetics(cbit.vcell.model.HMM_IRRKinetics)
GeneralCurrentLumpedKinetics(cbit.vcell.model.GeneralCurrentLumpedKinetics)
Microscopic_IRRKinetics(cbit.vcell.model.Microscopic_IRRKinetics)
ReactionParticipant(cbit.vcell.model.ReactionParticipant)
Aggregations
VariableHash (cbit.vcell.math.VariableHash)15
Expression (cbit.vcell.parser.Expression)14
Variable (cbit.vcell.math.Variable)12
CompartmentSubDomain (cbit.vcell.math.CompartmentSubDomain)11
Constant (cbit.vcell.math.Constant)11
Function (cbit.vcell.math.Function)11
MathDescription (cbit.vcell.math.MathDescription)9
Domain (cbit.vcell.math.Variable.Domain)9
VolVariable (cbit.vcell.math.VolVariable)9
ExpressionException (cbit.vcell.parser.ExpressionException)9
PropertyVetoException (java.beans.PropertyVetoException)9
MemVariable (cbit.vcell.math.MemVariable)8
MembraneRegionVariable (cbit.vcell.math.MembraneRegionVariable)8
MathException (cbit.vcell.math.MathException)7
MembraneSubDomain (cbit.vcell.math.MembraneSubDomain)7
SubDomain (cbit.vcell.math.SubDomain)7
ModelUnitSystem (cbit.vcell.model.ModelUnitSystem)7
VCUnitDefinition (cbit.vcell.units.VCUnitDefinition)7
Vector (java.util.Vector)7
SubVolume (cbit.vcell.geometry.SubVolume)6
