Example 1 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 2 with VariableHash
use of cbit.vcell.math.VariableHash in project vcell by virtualcell.
the class DiffEquMathMapping method refreshMathDescription.
/**
* This method was created in VisualAge.
*/
@SuppressWarnings("deprecation")
private void refreshMathDescription() throws MappingException, MatrixException, MathException, ExpressionException, ModelException {
// All sizes must be set for new ODE models and ratios must be set for old ones.
simContext.checkValidity();
//
// temporarily place all variables in a hashtable (before binding) and discarding duplicates (check for equality)
//
VariableHash varHash = new VariableHash();
StructureMapping[] structureMappings = simContext.getGeometryContext().getStructureMappings();
//
// 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.getGeometryClass() == null)){
// localIssueList.add(new Issue(structures[i], IssueCategory.StructureNotMapped,"In Application '" + simContext.getName() + "', model structure '"+structures[i].getName()+"' not mapped to a geometry subdomain",Issue.SEVERITY_WARNING));
// }
// }
// SubVolume subVolumes[] = simContext.getGeometryContext().getGeometry().getGeometrySpec().getSubVolumes();
// for (int i = 0; i < subVolumes.length; i++){
// Structure[] mappedStructures = simContext.getGeometryContext().getStructuresFromGeometryClass(subVolumes[i]);
// if (mappedStructures==null || mappedStructures.length==0){
// localIssueList.add(new Issue(subVolumes[i], IssueCategory.GeometryClassNotMapped,"In Application '" + simContext.getName() + "', geometry subVolume '"+subVolumes[i].getName()+"' not mapped from a model structure",Issue.SEVERITY_WARNING));
// }
// }
// deals with model parameters
HashMap<VolVariable, EventAssignmentOrRateRuleInitParameter> eventVolVarHash = new HashMap<VolVariable, EventAssignmentOrRateRuleInitParameter>();
HashMap<Variable, RateRuleRateParameter> rateRuleRateParamHash = new HashMap<Variable, RateRuleRateParameter>();
ArrayList<SymbolTableEntry> rateRuleVarTargets = new ArrayList<SymbolTableEntry>();
ArrayList<SymbolTableEntry> assignmentRuleVarTargets = new ArrayList<SymbolTableEntry>();
ArrayList<SymbolTableEntry> eventAssignTargets = new ArrayList<SymbolTableEntry>();
Model model = simContext.getModel();
ModelUnitSystem modelUnitSystem = model.getUnitSystem();
VCUnitDefinition timeUnit = modelUnitSystem.getTimeUnit();
ModelParameter[] modelParameters = model.getModelParameters();
if (simContext.getGeometry().getDimension() == 0) {
//
// global parameters from model (that presently are constants)
//
BioEvent[] bioEvents = simContext.getBioEvents();
if (bioEvents != null && bioEvents.length > 0) {
for (BioEvent be : bioEvents) {
ArrayList<EventAssignment> eventAssignments = be.getEventAssignments();
if (eventAssignments != null) {
for (EventAssignment ea : eventAssignments) {
if (!eventAssignTargets.contains(ea.getTarget())) {
eventAssignTargets.add(ea.getTarget());
}
}
}
}
}
RateRule[] rrs = simContext.getRateRules();
if (rrs != null && rrs.length > 0) {
for (RateRule rr : rrs) {
SymbolTableEntry rrVar = rr.getRateRuleVar();
if (!rateRuleVarTargets.contains(rrVar)) {
rateRuleVarTargets.add(rrVar);
}
}
}
AssignmentRule[] ars = simContext.getAssignmentRules();
if (ars != null && ars.length > 0) {
for (AssignmentRule ar : ars) {
SymbolTableEntry arVar = ar.getAssignmentRuleVar();
if (!assignmentRuleVarTargets.contains(arVar)) {
assignmentRuleVarTargets.add(arVar);
}
}
}
for (int j = 0; j < modelParameters.length; j++) {
Expression modelParamExpr = modelParameters[j].getExpression();
GeometryClass geometryClass = getDefaultGeometryClass(modelParamExpr);
VCUnitDefinition paramUnit = modelParameters[j].getUnitDefinition();
modelParamExpr = getIdentifierSubstitutions(modelParamExpr, paramUnit, geometryClass);
// if (eventAssignTargets.contains(modelParameters[j]) || rateRuleVarTargets.contains(modelParameters[j])) {
if (eventAssignTargets.contains(modelParameters[j])) {
EventAssignmentOrRateRuleInitParameter eap = null;
try {
eap = addEventAssignmentOrRateRuleInitParameter(modelParameters[j], modelParamExpr, PARAMETER_ROLE_EVENTASSIGN_OR_RATERULE_INITCONDN, paramUnit);
} catch (PropertyVetoException e) {
e.printStackTrace(System.out);
throw new MappingException(e.getMessage());
}
if (geometryClass == null) {
GeometryClass[] geometryClasses = simContext.getGeometryContext().getGeometry().getGeometryClasses();
geometryClass = geometryClasses[0];
}
Domain domain = null;
if (geometryClass != null) {
// the volume variable will look like Compartment::g0 rather than just g0
domain = new Domain(geometryClass);
}
VolVariable volVar = new VolVariable(modelParameters[j].getName(), domain);
varHash.addVariable(volVar);
eventVolVarHash.put(volVar, eap);
} else if (rateRuleVarTargets.contains(modelParameters[j])) {
// do nothing, will do elsewhere
;
} else if (assignmentRuleVarTargets.contains(modelParameters[j])) {
// do nothing, will do elsewhere
;
} else {
Variable variable = newFunctionOrConstant(getMathSymbol(modelParameters[j], geometryClass), modelParamExpr, geometryClass);
varHash.addVariable(variable);
}
}
} else {
for (int j = 0; j < modelParameters.length; j++) {
Expression modelParamExpr = modelParameters[j].getExpression();
GeometryClass geometryClass = getDefaultGeometryClass(modelParamExpr);
modelParamExpr = getIdentifierSubstitutions(modelParamExpr, modelParameters[j].getUnitDefinition(), geometryClass);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], geometryClass), modelParamExpr, geometryClass));
}
}
//
for (SimulationContextParameter scParameter : simContext.getSimulationContextParameters()) {
Expression scParameterExpression = scParameter.getExpression();
GeometryClass gc = getDefaultGeometryClass(scParameterExpression);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(scParameter, gc), getIdentifierSubstitutions(scParameter.getExpression(), scParameter.getUnitDefinition(), gc), gc));
}
//
for (DataSymbol dataSymbol : simContext.getDataContext().getDataSymbols()) {
if (dataSymbol instanceof FieldDataSymbol) {
FieldDataSymbol fieldDataSymbol = (FieldDataSymbol) dataSymbol;
GeometryClass geometryClass = null;
FieldFunctionArguments ffs = new FieldFunctionArguments(fieldDataSymbol.getExternalDataIdentifier().getName(), fieldDataSymbol.getFieldDataVarName(), new Expression(fieldDataSymbol.getFieldDataVarTime()), VariableType.getVariableTypeFromVariableTypeName(fieldDataSymbol.getFieldDataVarType()));
Expression exp = new Expression(ffs.infix());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(dataSymbol, geometryClass), getIdentifierSubstitutions(exp, dataSymbol.getUnitDefinition(), geometryClass), geometryClass));
} else {
throw new RuntimeException("In Application '" + simContext.getName() + "', dataSymbol type '" + dataSymbol.getClass().getName() + "' not yet supported for math generation");
}
}
//
// 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("In Application '" + simContext.getName() + "', " + 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");
}
//
// volume variables
//
Enumeration<SpeciesContextMapping> enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = enum1.nextElement();
if (scm.getVariable() instanceof VolVariable) {
if (!(mathDesc.getVariable(scm.getVariable().getName()) instanceof VolVariable)) {
varHash.addVariable(scm.getVariable());
}
}
}
//
// membrane variables
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() instanceof MemVariable) {
varHash.addVariable(scm.getVariable());
}
}
//
// volume region variables
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() instanceof VolumeRegionVariable) {
varHash.addVariable(scm.getVariable());
}
}
//
// membrane region variables
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() instanceof MembraneRegionVariable) {
varHash.addVariable(scm.getVariable());
}
}
//
// add compartment and membrane subdomains
//
ArrayList<CompartmentSubdomainContext> compartmentSubdomainContexts = new ArrayList<CompartmentSubdomainContext>();
ArrayList<MembraneSubdomainContext> membraneSubdomainContexts = new ArrayList<MembraneSubdomainContext>();
addSubdomains(model, compartmentSubdomainContexts, membraneSubdomainContexts);
// membrane velocities set on MembraneSubdomains later.
addSpatialProcesses(varHash, compartmentSubdomainContexts, membraneSubdomainContexts);
varHash.addVariable(new Constant(getMathSymbol(model.getPI_CONSTANT(), null), getIdentifierSubstitutions(model.getPI_CONSTANT().getExpression(), model.getPI_CONSTANT().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)));
//
// only calculate potential if at least one MembraneMapping has CalculateVoltage == true
//
boolean bCalculatePotential = false;
for (int i = 0; i < structureMappings.length; i++) {
if (structureMappings[i] instanceof MembraneMapping) {
if (((MembraneMapping) structureMappings[i]).getCalculateVoltage()) {
bCalculatePotential = true;
}
}
}
potentialMapping = new PotentialMapping(simContext, this);
if (bCalculatePotential) {
potentialMapping.computeMath();
//
// copy functions for currents and constants for capacitances
//
ElectricalDevice[] devices = potentialMapping.getElectricalDevices();
for (int j = 0; j < devices.length; j++) {
if (devices[j] instanceof MembraneElectricalDevice) {
MembraneElectricalDevice membraneElectricalDevice = (MembraneElectricalDevice) devices[j];
MembraneMapping memMapping = membraneElectricalDevice.getMembraneMapping();
Parameter specificCapacitanceParm = memMapping.getParameterFromRole(MembraneMapping.ROLE_SpecificCapacitance);
varHash.addVariable(new Constant(getMathSymbol(specificCapacitanceParm, memMapping.getGeometryClass()), getIdentifierSubstitutions(specificCapacitanceParm.getExpression(), specificCapacitanceParm.getUnitDefinition(), memMapping.getGeometryClass())));
ElectricalDevice.ElectricalDeviceParameter transmembraneCurrentParm = membraneElectricalDevice.getParameterFromRole(ElectricalDevice.ROLE_TransmembraneCurrent);
ElectricalDevice.ElectricalDeviceParameter totalCurrentParm = membraneElectricalDevice.getParameterFromRole(ElectricalDevice.ROLE_TotalCurrent);
ElectricalDevice.ElectricalDeviceParameter capacitanceParm = membraneElectricalDevice.getParameterFromRole(ElectricalDevice.ROLE_Capacitance);
GeometryClass geometryClass = membraneElectricalDevice.getMembraneMapping().getGeometryClass();
if (totalCurrentParm != null && /* totalCurrentDensityParm.getExpression()!=null && */
memMapping.getCalculateVoltage()) {
Expression totalCurrentDensityExp = (totalCurrentParm.getExpression() != null) ? (totalCurrentParm.getExpression()) : (new Expression(0.0));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrentParm, geometryClass), getIdentifierSubstitutions(totalCurrentDensityExp, totalCurrentParm.getUnitDefinition(), geometryClass), geometryClass));
}
if (transmembraneCurrentParm != null && transmembraneCurrentParm.getExpression() != null && memMapping.getCalculateVoltage()) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(transmembraneCurrentParm, geometryClass), getIdentifierSubstitutions(transmembraneCurrentParm.getExpression(), transmembraneCurrentParm.getUnitDefinition(), geometryClass), geometryClass));
}
if (capacitanceParm != null && capacitanceParm.getExpression() != null && memMapping.getCalculateVoltage()) {
StructureMappingParameter sizeParameter = membraneElectricalDevice.getMembraneMapping().getSizeParameter();
if (simContext.getGeometry().getDimension() == 0 && (sizeParameter.getExpression() == null || sizeParameter.getExpression().isZero())) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(capacitanceParm, geometryClass), getIdentifierSubstitutions(Expression.mult(memMapping.getNullSizeParameterValue(), specificCapacitanceParm.getExpression()), capacitanceParm.getUnitDefinition(), geometryClass), geometryClass));
} else {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(capacitanceParm, geometryClass), getIdentifierSubstitutions(capacitanceParm.getExpression(), capacitanceParm.getUnitDefinition(), geometryClass), geometryClass));
}
}
//
if (membraneElectricalDevice.getDependentVoltageExpression() == null) {
// is Voltage Independent?
StructureMapping.StructureMappingParameter initialVoltageParm = memMapping.getInitialVoltageParameter();
varHash.addVariable(newFunctionOrConstant(getMathSymbol(initialVoltageParm, memMapping.getGeometryClass()), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), memMapping.getGeometryClass()), memMapping.getGeometryClass()));
} else //
// membrane forced potential
//
{
varHash.addVariable(newFunctionOrConstant(getMathSymbol(memMapping.getMembrane().getMembraneVoltage(), memMapping.getGeometryClass()), getIdentifierSubstitutions(membraneElectricalDevice.getDependentVoltageExpression(), memMapping.getMembrane().getMembraneVoltage().getUnitDefinition(), memMapping.getGeometryClass()), memMapping.getGeometryClass()));
}
} else if (devices[j] instanceof CurrentClampElectricalDevice) {
CurrentClampElectricalDevice currentClampDevice = (CurrentClampElectricalDevice) devices[j];
// total current = current source (no capacitance)
Parameter totalCurrentParm = currentClampDevice.getParameterFromRole(CurrentClampElectricalDevice.ROLE_TotalCurrent);
Parameter currentParm = currentClampDevice.getParameterFromRole(CurrentClampElectricalDevice.ROLE_TransmembraneCurrent);
// Parameter dependentVoltage = currentClampDevice.getCurrentClampStimulus().getVoltageParameter();
Feature deviceElectrodeFeature = currentClampDevice.getCurrentClampStimulus().getElectrode().getFeature();
Feature groundElectrodeFeature = simContext.getGroundElectrode().getFeature();
Membrane membrane = model.getStructureTopology().getMembrane(deviceElectrodeFeature, groundElectrodeFeature);
GeometryClass geometryClass = null;
if (membrane != null) {
StructureMapping membraneStructureMapping = simContext.getGeometryContext().getStructureMapping(membrane);
geometryClass = membraneStructureMapping.getGeometryClass();
}
varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrentParm, geometryClass), getIdentifierSubstitutions(totalCurrentParm.getExpression(), totalCurrentParm.getUnitDefinition(), geometryClass), geometryClass));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(currentParm, geometryClass), getIdentifierSubstitutions(currentParm.getExpression(), currentParm.getUnitDefinition(), geometryClass), geometryClass));
// varHash.addVariable(newFunctionOrConstant(getMathSymbol(dependentVoltage,null),getIdentifierSubstitutions(currentClampDevice.getDependentVoltageExpression(),dependentVoltage.getUnitDefinition(),null)));
//
// add user-defined parameters
//
ElectricalDevice.ElectricalDeviceParameter[] parameters = currentClampDevice.getParameters();
for (int k = 0; k < parameters.length; k++) {
if (parameters[k].getExpression() != null) {
// guards against voltage parameters that are "variable".
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[k], null), getIdentifierSubstitutions(parameters[k].getExpression(), parameters[k].getUnitDefinition(), geometryClass), geometryClass));
}
}
} else if (devices[j] instanceof VoltageClampElectricalDevice) {
VoltageClampElectricalDevice voltageClampDevice = (VoltageClampElectricalDevice) devices[j];
Feature deviceElectrodeFeature = voltageClampDevice.getVoltageClampStimulus().getElectrode().getFeature();
Feature groundElectrodeFeature = simContext.getGroundElectrode().getFeature();
Membrane membrane = model.getStructureTopology().getMembrane(deviceElectrodeFeature, groundElectrodeFeature);
GeometryClass geometryClass = null;
if (membrane != null) {
StructureMapping membraneStructureMapping = simContext.getGeometryContext().getStructureMapping(membrane);
geometryClass = membraneStructureMapping.getGeometryClass();
}
// total current = current source (no capacitance)
Parameter totalCurrent = voltageClampDevice.getParameterFromRole(VoltageClampElectricalDevice.ROLE_TotalCurrent);
Parameter totalCurrentParm = voltageClampDevice.getParameterFromRole(VoltageClampElectricalDevice.ROLE_TotalCurrent);
Parameter voltageParm = voltageClampDevice.getParameterFromRole(VoltageClampElectricalDevice.ROLE_Voltage);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrent, geometryClass), getIdentifierSubstitutions(totalCurrent.getExpression(), totalCurrent.getUnitDefinition(), geometryClass), geometryClass));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrentParm, geometryClass), getIdentifierSubstitutions(totalCurrentParm.getExpression(), totalCurrentParm.getUnitDefinition(), geometryClass), geometryClass));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(voltageParm, geometryClass), getIdentifierSubstitutions(voltageParm.getExpression(), voltageParm.getUnitDefinition(), geometryClass), geometryClass));
//
// add user-defined parameters
//
ElectricalDevice.ElectricalDeviceParameter[] parameters = voltageClampDevice.getParameters();
for (int k = 0; k < parameters.length; k++) {
if (parameters[k].getRole() == ElectricalDevice.ROLE_UserDefined) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[k], geometryClass), getIdentifierSubstitutions(parameters[k].getExpression(), parameters[k].getUnitDefinition(), geometryClass), geometryClass));
}
}
}
}
} else {
//
for (int j = 0; j < structureMappings.length; j++) {
if (structureMappings[j] instanceof MembraneMapping) {
MembraneMapping memMapping = (MembraneMapping) structureMappings[j];
varHash.addVariable(newFunctionOrConstant(getMathSymbol(memMapping.getMembrane().getMembraneVoltage(), memMapping.getGeometryClass()), getIdentifierSubstitutions(memMapping.getInitialVoltageParameter().getExpression(), memMapping.getInitialVoltageParameter().getUnitDefinition(), memMapping.getGeometryClass()), memMapping.getGeometryClass()));
}
}
}
//
for (int j = 0; j < structureMappings.length; j++) {
if (structureMappings[j] instanceof MembraneMapping) {
MembraneMapping membraneMapping = (MembraneMapping) structureMappings[j];
Membrane.MembraneVoltage membraneVoltage = membraneMapping.getMembrane().getMembraneVoltage();
ElectricalDevice[] membraneDevices = potentialMapping.getElectricalDevices(membraneMapping.getMembrane());
// ElectricalDevice membraneDevice = null;
for (int i = 0; i < membraneDevices.length; i++) {
if (membraneDevices[i].hasCapacitance() && membraneDevices[i].getDependentVoltageExpression() == null) {
GeometryClass geometryClass = membraneMapping.getGeometryClass();
if (geometryClass == null) {
throw new MappingException("Application '" + getSimulationContext().getName() + "'\nGeometry->StructureMapping->(" + structureMappings[j].getStructure().getTypeName() + ")'" + structureMappings[j].getStructure().getName() + "' must be mapped to geometry domain.\n(see 'Problems' tab)");
}
Domain domain = new Domain(geometryClass);
if (membraneMapping.getCalculateVoltage() && bCalculatePotential) {
if (geometryClass instanceof SurfaceClass) {
//
if (mathDesc.getVariable(Membrane.MEMBRANE_VOLTAGE_REGION_NAME) == null) {
// varHash.addVariable(new MembraneRegionVariable(MembraneVoltage.MEMBRANE_VOLTAGE_REGION_NAME));
varHash.addVariable(new MembraneRegionVariable(getMathSymbol(membraneVoltage, geometryClass), domain));
}
} else {
//
// spatially unresolved membrane, and must solve for potential ... make VolVariable for this compartment
//
varHash.addVariable(new VolVariable(getMathSymbol(membraneVoltage, geometryClass), domain));
}
Parameter initialVoltageParm = membraneMapping.getInitialVoltageParameter();
Variable initVoltageFunction = newFunctionOrConstant(getMathSymbol(initialVoltageParm, geometryClass), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), geometryClass), geometryClass);
varHash.addVariable(initVoltageFunction);
} else {
//
// don't calculate voltage, still may need it though
//
Parameter initialVoltageParm = membraneMapping.getInitialVoltageParameter();
Variable voltageFunction = newFunctionOrConstant(getMathSymbol(membraneMapping.getMembrane().getMembraneVoltage(), geometryClass), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), geometryClass), geometryClass);
varHash.addVariable(voltageFunction);
}
}
}
}
}
//
for (int j = 0; j < reactionSteps.length; j++) {
ReactionStep rs = reactionSteps[j];
if (simContext.getReactionContext().getReactionSpec(rs).isExcluded()) {
continue;
}
Kinetics.KineticsParameter[] parameters = rs.getKinetics().getKineticsParameters();
GeometryClass geometryClass = null;
if (rs.getStructure() != null) {
geometryClass = simContext.getGeometryContext().getStructureMapping(rs.getStructure()).getGeometryClass();
}
if (parameters != null) {
for (int i = 0; i < parameters.length; i++) {
if (((parameters[i].getRole() == Kinetics.ROLE_CurrentDensity) || (parameters[i].getRole() == Kinetics.ROLE_LumpedCurrent)) && (parameters[i].getExpression() == null || parameters[i].getExpression().isZero())) {
continue;
}
String mathSymbol = getMathSymbol(parameters[i], geometryClass);
Expression expr = getIdentifierSubstitutions(parameters[i].getExpression(), parameters[i].getUnitDefinition(), geometryClass);
varHash.addVariable(newFunctionOrConstant(mathSymbol, expr, geometryClass));
}
}
}
//
// initial conditions (either function or constant) for rate rule variables that are model parameters
//
// the init variables with expressions still containing variables
Map<ModelParameter, Variable> initModelParameterHashTmp = new HashMap<>();
// here we store the init parameter of the model parameter
Map<EventAssignmentOrRateRuleInitParameter, ModelParameter> rateRuleInitToModelParamMapping = new HashMap<>();
// here we store the init parameter of the model parameter
Map<ModelParameter, EventAssignmentOrRateRuleInitParameter> modelParamTorateRuleInitMapping = new HashMap<>();
for (ModelParameter mp : modelParameters) {
// initial assignment for global parameter used as rate rule variable
RateRule rr = simContext.getRateRule(mp);
if (rr == null) {
// we only care about global parameters that are rate rule variables
continue;
}
Variable var = varHash.getVariable(mp.getName());
if (var != null) {
if (eventVolVarHash.containsKey(var)) {
System.out.println("Global Parameters that are rate rule Variables should be unmapped at this point, unless they are EventAssignments too.");
} else {
throw new MappingException("Global Parameters that are rate rule Variables should be unmapped at this point.");
}
}
Expression modelParamExpr = mp.getExpression();
if (modelParamExpr == null) {
continue;
}
GeometryClass gc = getDefaultGeometryClass(modelParamExpr);
VCUnitDefinition paramUnit = modelUnitSystem.getInstance_TBD();
if (mp.getUnitDefinition() != null && !mp.getUnitDefinition().equals(modelUnitSystem.getInstance_TBD())) {
paramUnit = mp.getUnitDefinition();
}
// TODO: is this really needed? or could I directly use modelParamExpr in addEventAssignmentOrRateRuleInitParameter()
Expression mpInitExpr = getIdentifierSubstitutions(modelParamExpr, paramUnit, gc);
EventAssignmentOrRateRuleInitParameter mpInitParam;
try {
mpInitParam = addEventAssignmentOrRateRuleInitParameter(mp, mpInitExpr, PARAMETER_ROLE_EVENTASSIGN_OR_RATERULE_INITCONDN, paramUnit);
} catch (PropertyVetoException e) {
e.printStackTrace(System.out);
throw new MappingException(e.getMessage());
}
rateRuleInitToModelParamMapping.put(mpInitParam, mp);
modelParamTorateRuleInitMapping.put(mp, mpInitParam);
}
//
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
fieldMathMappingParameters[i].getExpression().bindExpression(this);
Expression exp = getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass);
Variable var = newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), exp, geometryClass);
varHash.addVariable(var);
ModelParameter mp = rateRuleInitToModelParamMapping.get(fieldMathMappingParameters[i]);
if (mp != null) {
initModelParameterHashTmp.put(mp, var);
}
}
//
// initial conditions (either function or constant) for species variables
//
SpeciesContextSpec[] speciesContextSpecs = simContext.getReactionContext().getSpeciesContextSpecs();
for (int i = 0; i < speciesContextSpecs.length; i++) {
// add initial count if present (!= null)
SpeciesContextSpecParameter initCountParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialCount);
SpeciesContext speciesContext = speciesContextSpecs[i].getSpeciesContext();
if (initCountParm != null && initCountParm.getExpression() != null) {
Expression initCountExpr = new Expression(initCountParm.getExpression());
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContext.getStructure());
String[] symbols = initCountExpr.getSymbols();
// Check if 'initExpr' has other speciesContexts in its expression, need to replace it with 'spContext_init'
for (int j = 0; symbols != null && j < symbols.length; j++) {
// if symbol is a speciesContext, replacing it with a reference to initial condition for that speciesContext.
SpeciesContext spC = null;
SymbolTableEntry ste = initCountExpr.getSymbolBinding(symbols[j]);
if (ste instanceof SpeciesContextSpecProxyParameter) {
SpeciesContextSpecProxyParameter spspp = (SpeciesContextSpecProxyParameter) ste;
if (spspp.getTarget() instanceof SpeciesContext) {
spC = (SpeciesContext) spspp.getTarget();
SpeciesContextSpec spcspec = simContext.getReactionContext().getSpeciesContextSpec(spC);
SpeciesContextSpecParameter spCInitParm = spcspec.getParameterFromRole(SpeciesContextSpec.ROLE_InitialCount);
// need to get init condn expression, but can't get it from getMathSymbol() (mapping between bio and math), hence get it as below.
Expression scsInitExpr = new Expression(spCInitParm, getNameScope());
initCountExpr.substituteInPlace(new Expression(spC.getName()), scsInitExpr);
}
}
}
// now create the appropriate function for the current speciesContextSpec.
varHash.addVariable(newFunctionOrConstant(getMathSymbol(initCountParm, sm.getGeometryClass()), getIdentifierSubstitutions(initCountExpr, initCountParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
// add initial concentration (may be derived from initial count if necessary)
SpeciesContextSpecParameter initConcParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
if (initConcParm != null) {
Expression initConcExpr = null;
if (initConcParm.getExpression() != null) {
initConcExpr = new Expression(initConcParm.getExpression());
} else if (initCountParm != null && initCountParm.getExpression() != null) {
Expression structureSizeExpr = new Expression(speciesContext.getStructure().getStructureSize(), getNameScope());
VCUnitDefinition concUnit = initConcParm.getUnitDefinition();
VCUnitDefinition countDensityUnit = initCountParm.getUnitDefinition().divideBy(speciesContext.getStructure().getStructureSize().getUnitDefinition());
Expression unitFactor = getUnitFactor(concUnit.divideBy(countDensityUnit));
initConcExpr = Expression.mult(Expression.div(new Expression(initCountParm, getNameScope()), structureSizeExpr), unitFactor);
}
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContext.getStructure());
String[] symbols = initConcExpr.getSymbols();
// Check if 'initExpr' has other speciesContexts in its expression, need to replace it with 'spContext_init'
for (int j = 0; symbols != null && j < symbols.length; j++) {
// if symbol is a speciesContext, replacing it with a reference to initial condition for that speciesContext.
SpeciesContext spC = null;
SymbolTableEntry ste = initConcExpr.getSymbolBinding(symbols[j]);
if (ste == null) {
String msg = initConcParm.getName() == null ? "??" : initConcParm.getName();
System.out.println("Unexpected NULL symbol in the initial expression of " + msg);
} else if (ste instanceof SpeciesContextSpecProxyParameter) {
SpeciesContextSpecProxyParameter spspp = (SpeciesContextSpecProxyParameter) ste;
if (spspp.getTarget() instanceof SpeciesContext) {
spC = (SpeciesContext) spspp.getTarget();
SpeciesContextSpec spcspec = simContext.getReactionContext().getSpeciesContextSpec(spC);
SpeciesContextSpecParameter spCInitParm = spcspec.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
// if initConc param expression is null, try initCount
if (spCInitParm.getExpression() == null) {
spCInitParm = spcspec.getParameterFromRole(SpeciesContextSpec.ROLE_InitialCount);
}
// need to get init condn expression, but can't get it from getMathSymbol() (mapping between bio and math), hence get it as below.
Expression scsInitExpr = new Expression(spCInitParm, getNameScope());
initConcExpr.substituteInPlace(new Expression(spC.getName()), scsInitExpr);
}
} else if (ste instanceof ModelParameter) {
ModelParameter mpArg = (ModelParameter) ste;
System.out.println(mpArg.getName());
if (simContext.getRateRule(mpArg) == null) {
// only globals that are RateRule variables need to be replaced with their _init variable
continue;
}
EventAssignmentOrRateRuleInitParameter mpInitParam = modelParamTorateRuleInitMapping.get(mpArg);
if (mpInitParam != null) {
// we already made it, we only need to use it
Expression mpArgInitExpr = new Expression(mpInitParam, getNameScope());
initConcExpr.substituteInPlace(new Expression(ste.getName()), mpArgInitExpr);
}
} else {
String msg = ste.getName() == null ? "??" : ste.getName();
String msg2 = initConcParm.getName() == null ? "??" : initConcParm.getName();
System.out.println("Unexpected symbol type for " + msg + " in the initial expression of " + msg2);
}
}
// now create the appropriate function for the current speciesContextSpec.
varHash.addVariable(newFunctionOrConstant(getMathSymbol(initConcParm, sm.getGeometryClass()), getIdentifierSubstitutions(initConcExpr, initConcParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextMapping scm = getSpeciesContextMapping(speciesContextSpecs[i].getSpeciesContext());
SpeciesContextSpec.SpeciesContextSpecParameter diffParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_DiffusionRate);
if (diffParm != null && (scm.isPDERequired())) {
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(diffParm, sm.getGeometryClass()), getIdentifierSubstitutions(diffParm.getExpression(), diffParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter bc_xm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueXm);
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
if (bc_xm != null && (bc_xm.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_xm, sm.getGeometryClass()), getIdentifierSubstitutions(bc_xm.getExpression(), bc_xm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_xp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueXp);
if (bc_xp != null && (bc_xp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_xp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_xp.getExpression(), bc_xp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_ym = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueYm);
if (bc_ym != null && (bc_ym.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_ym, sm.getGeometryClass()), getIdentifierSubstitutions(bc_ym.getExpression(), bc_ym.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_yp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueYp);
if (bc_yp != null && (bc_yp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_yp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_yp.getExpression(), bc_yp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_zm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueZm);
if (bc_zm != null && (bc_zm.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_zm, sm.getGeometryClass()), getIdentifierSubstitutions(bc_zm.getExpression(), bc_zm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_zp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueZp);
if (bc_zp != null && (bc_zp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_zp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_zp.getExpression(), bc_zp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter advection_velX = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityX);
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
GeometryClass geometryClass = sm.getGeometryClass();
if (advection_velX != null && (advection_velX.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velX, geometryClass), getIdentifierSubstitutions(advection_velX.getExpression(), advection_velX.getUnitDefinition(), geometryClass), geometryClass));
}
SpeciesContextSpec.SpeciesContextSpecParameter advection_velY = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityY);
if (advection_velY != null && (advection_velY.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velY, geometryClass), getIdentifierSubstitutions(advection_velY.getExpression(), advection_velY.getUnitDefinition(), geometryClass), geometryClass));
}
SpeciesContextSpec.SpeciesContextSpecParameter advection_velZ = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityZ);
if (advection_velZ != null && (advection_velZ.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velZ, geometryClass), getIdentifierSubstitutions(advection_velZ.getExpression(), advection_velZ.getUnitDefinition(), geometryClass), geometryClass));
}
}
//
// 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());
}
}
//
// conversion factors
//
varHash.addVariable(new Constant(model.getKMOLE().getName(), getIdentifierSubstitutions(model.getKMOLE().getExpression(), model.getKMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(model.getN_PMOLE().getName(), getIdentifierSubstitutions(model.getN_PMOLE().getExpression(), model.getN_PMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(model.getKMILLIVOLTS().getName(), getIdentifierSubstitutions(model.getKMILLIVOLTS().getExpression(), model.getKMILLIVOLTS().getUnitDefinition(), null)));
varHash.addVariable(new Constant(model.getK_GHK().getName(), getIdentifierSubstitutions(model.getK_GHK().getExpression(), model.getK_GHK().getUnitDefinition(), null)));
//
for (int i = 0; i < structureMappings.length; i++) {
StructureMapping sm = structureMappings[i];
if (simContext.getGeometry().getDimension() == 0) {
StructureMappingParameter sizeParm = sm.getSizeParameter();
if (sizeParm != null && sizeParm.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(sizeParm, sm.getGeometryClass()), getIdentifierSubstitutions(sizeParm.getExpression(), sizeParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
} else {
if (sm instanceof MembraneMapping) {
MembraneMapping mm = (MembraneMapping) sm;
StructureMappingParameter volFrac = mm.getVolumeFractionParameter();
if (volFrac != null && volFrac.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(volFrac, sm.getGeometryClass()), getIdentifierSubstitutions(volFrac.getExpression(), volFrac.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
StructureMappingParameter surfToVol = mm.getSurfaceToVolumeParameter();
if (surfToVol != null && surfToVol.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(surfToVol, sm.getGeometryClass()), getIdentifierSubstitutions(surfToVol.getExpression(), surfToVol.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
}
} else {
Parameter parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitArea);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitVolume);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitArea);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitVolume);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
StructureMappingParameter sizeParm = sm.getSizeParameter();
if (sm.getGeometryClass() != null && sizeParm != null) {
if (simContext.getGeometry().getDimension() == 0) {
if (sizeParm.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(sizeParm, sm.getGeometryClass()), getIdentifierSubstitutions(sizeParm.getExpression(), sizeParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
} else {
String compartmentName = sm.getGeometryClass().getName();
VCUnitDefinition sizeUnit = sm.getSizeParameter().getUnitDefinition();
String sizeFunctionName = null;
if (sm instanceof MembraneMapping) {
MembraneMapping mm = (MembraneMapping) sm;
if (mm.getGeometryClass() instanceof SurfaceClass) {
sizeFunctionName = MathFunctionDefinitions.Function_regionArea_current.getFunctionName();
} else if (mm.getGeometryClass() instanceof SubVolume) {
sizeFunctionName = MathFunctionDefinitions.Function_regionVolume_current.getFunctionName();
}
} else if (sm instanceof FeatureMapping) {
sizeFunctionName = MathFunctionDefinitions.Function_regionVolume_current.getFunctionName();
} else {
throw new RuntimeException("structure mapping " + sm.getClass().getName() + " not yet supported");
}
Expression totalVolumeCorrection = sm.getStructureSizeCorrection(simContext, this);
Expression sizeFunctionExpression = Expression.function(sizeFunctionName, new Expression[] { new Expression("'" + compartmentName + "'") });
// sizeFunctionExpression.bindExpression(mathDesc);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(sizeParm, sm.getGeometryClass()), getIdentifierSubstitutions(Expression.mult(totalVolumeCorrection, sizeFunctionExpression), sizeUnit, sm.getGeometryClass()), sm.getGeometryClass()));
}
}
}
//
// functions
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() == null && scm.getDependencyExpression() != null) {
// check if speciesContext has a rateRule; then the speciesContext should not be added as a constant
if (simContext.getRateRule(scm.getSpeciesContext()) == null) {
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(scm.getSpeciesContext().getStructure());
if (sm.getGeometryClass() == null) {
Structure s = sm.getStructure();
if (s != null) {
throw new RuntimeException("unmapped structure " + s.getName());
}
throw new RuntimeException("structure mapping with no structure or mapping");
}
Variable dependentVariable = newFunctionOrConstant(getMathSymbol(scm.getSpeciesContext(), sm.getGeometryClass()), getIdentifierSubstitutions(scm.getDependencyExpression(), scm.getSpeciesContext().getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass());
dependentVariable.setDomain(new Domain(sm.getGeometryClass()));
varHash.addVariable(dependentVariable);
}
}
}
BioEvent[] bioevents = simContext.getBioEvents();
if (bioevents != null && bioevents.length > 0) {
for (BioEvent be : bioevents) {
// transform the bioEvent trigger/delay to math Event
for (LocalParameter p : be.getEventParameters()) {
if (p.getExpression() != null) {
// ex: eventName.delay and eventName.triggerFunction
String name = getMathSymbol(p, null);
Expression exp = getIdentifierSubstitutions(p.getExpression(), p.getUnitDefinition(), null);
Variable var = newFunctionOrConstant(name, exp, null);
varHash.addVariable(var);
} else if (be.getParameter(BioEventParameterType.GeneralTriggerFunction) == p) {
//
// use generated function here.
//
String name = getMathSymbol(p, null);
Expression exp = getIdentifierSubstitutions(be.generateTriggerExpression(), p.getUnitDefinition(), null);
Variable var = newFunctionOrConstant(name, exp, null);
varHash.addVariable(var);
}
}
}
}
//
// substitute init functions for event assignment variables
//
// for (Map.Entry<VolVariable,EventAssignmentOrRateRuleInitParameter> entry : eventVolVarHash.entrySet()) {
// EventAssignmentOrRateRuleInitParameter eap = entry.getValue();
//
// String argName = eap.getName();
// Expression modelParamExpr = eap.getExpression();
// GeometryClass gc = getDefaultGeometryClass(modelParamExpr);
// VCUnitDefinition paramUnit = eap.getUnitDefinition();
// Expression mpInitExpr = new Expression(modelParamExpr);
// String[] symbols = mpInitExpr.getSymbols();
// if(symbols == null || symbols.length == 0) {
// continue;
// }
// // TODO: this is still not working well
// // check if 'initExpr' has other speciesContexts or rate rule global parameter variables in its expression
// // need to replace it with 'spContext_init', modelParameter_init
// for (String symbol : symbols) {
// // if symbol is a speciesContext, replacing it with a reference to initial condition for that speciesContext.
// SymbolTableEntry ste = mpInitExpr.getSymbolBinding(symbol);
// if (ste == null) {
// System.out.println("Unexpected NULL symbol in the initial expression of " + argName);
// } else if (ste instanceof SpeciesContext) {
// SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec((SpeciesContext)ste);
// // TODO: what if initial count???
// SpeciesContextSpecParameter spCInitParm = scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
// // need to get init condn expression, but can't get it from getMathSymbol() (mapping between bio and math), hence get it as below.
// Expression scsInitExpr = new Expression(spCInitParm, getNameScope());
// mpInitExpr.substituteInPlace(new Expression(ste.getName()), scsInitExpr);
// } else if(ste instanceof ModelParameter) {
// ModelParameter mpArg = (ModelParameter)ste;
// System.out.println(mpArg.getName());
// if(simContext.getRateRule(mpArg) == null) {
// continue; // only globals that are RateRule variables need to be replaced with their _init variable
// }
// Variable mpArgVar = initModelParameterHashTmp.get(mpArg);
// if(mpArgVar != null && eventVolVarHash.get(mpArgVar) != null) {
// EventAssignmentOrRateRuleInitParameter mpInitParam = eventVolVarHash.get(mpArgVar);
// Expression mpArgInitExpr = new Expression(mpInitParam, getNameScope());
// mpInitExpr.substituteInPlace(new Expression(ste.getName()), mpArgInitExpr);
//
// }
// } else {
// String msg = ste.getName() == null ? "??" : ste.getName();
// System.out.println("Unexpected symbol type for " + msg + " in the initial expression of " + argName);
// }
// }
// varHash.removeVariable(argName);
// Expression exp = getIdentifierSubstitutions(mpInitExpr, paramUnit, gc);
// Variable varInit = newFunctionOrConstant(argName, exp, gc);
// varHash.addVariable(varInit);
// }
//
// deal with rate rules
//
// first, substitute the init functions for rate rule variables that are model parameters
// we'll need this init variable (function or constant) for the ODE Equation
//
// here we store the init variable with the final substitutions within their expressions
Map<ModelParameter, Variable> initModelParameterHash = new HashMap<>();
Map<String, SymbolTableEntry> entryMap = new HashMap<String, SymbolTableEntry>();
simContext.getEntries(entryMap);
for (Map.Entry<ModelParameter, Variable> entry : initModelParameterHashTmp.entrySet()) {
ModelParameter mp = entry.getKey();
Variable mpInitVariable = entry.getValue();
String argName = mpInitVariable.getName();
Expression modelParamExpr = mp.getExpression();
GeometryClass gc = getDefaultGeometryClass(modelParamExpr);
Expression mpInitExpr = new Expression(modelParamExpr);
String[] symbols = mpInitExpr.getSymbols();
if (symbols == null || symbols.length == 0) {
// stays as it is in variable hash
// we just move it into the initModelParameterHash
initModelParameterHash.put(mp, mpInitVariable);
continue;
}
// need to replace it with 'spContext_init', modelParameter_init
for (String symbol : symbols) {
// if symbol is a speciesContext, replacing it with a reference to initial condition for that speciesContext.
SymbolTableEntry ste = mpInitExpr.getSymbolBinding(symbol);
if (ste == null) {
System.out.println("Unexpected NULL symbol in the initial expression of " + argName);
} else if (ste instanceof SpeciesContext) {
SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec((SpeciesContext) ste);
// TODO: what if initial count???
SpeciesContextSpecParameter spCInitParm = scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
// need to get init condn expression, but can't get it from getMathSymbol() (mapping between bio and math), hence get it as below.
Expression scsInitExpr = new Expression(spCInitParm, getNameScope());
mpInitExpr.substituteInPlace(new Expression(ste.getName()), scsInitExpr);
} else if (ste instanceof ModelParameter) {
ModelParameter mpArg = (ModelParameter) ste;
System.out.println(mpArg.getName());
if (simContext.getRateRule(mpArg) == null) {
// only globals that are RateRule variables need to be replaced with their _init variable
continue;
}
EventAssignmentOrRateRuleInitParameter mpInitParam = modelParamTorateRuleInitMapping.get(mpArg);
if (mpInitParam != null) {
// we already made it, we only need to use it
Expression mpArgInitExpr = new Expression(mpInitParam, getNameScope());
mpInitExpr.substituteInPlace(new Expression(ste.getName()), mpArgInitExpr);
}
} else {
String msg = ste.getName() == null ? "??" : ste.getName();
System.out.println("Unexpected symbol type for " + msg + " in the initial expression of " + argName);
}
}
VCUnitDefinition paramUnit = modelUnitSystem.getInstance_TBD();
if (mp.getUnitDefinition() != null && !mp.getUnitDefinition().equals(modelUnitSystem.getInstance_TBD())) {
paramUnit = mp.getUnitDefinition();
}
varHash.removeVariable(mpInitVariable);
Expression exp = getIdentifierSubstitutions(mpInitExpr, paramUnit, gc);
mpInitVariable = newFunctionOrConstant(argName, exp, gc);
varHash.addVariable(mpInitVariable);
initModelParameterHash.put(mp, mpInitVariable);
}
//
// create the VolVariable for the species context used as rate rule variable
// create the Variable (function or constant) for its rate (need it for the ODE Equation)
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
// species context used as rate rule variable
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
Variable var = scm.getVariable();
Expression exp = scm.getDependencyExpression();
if (var == null && exp != null) {
RateRule rr = simContext.getRateRule(scm.getSpeciesContext());
if (rr != null && (rr.getRateRuleVar() instanceof SpeciesContext)) {
SpeciesContext sc = scm.getSpeciesContext();
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(sc.getStructure());
if (sm.getGeometryClass() == null) {
Structure s = sm.getStructure();
if (s != null) {
throw new RuntimeException("unmapped structure " + s.getName());
}
throw new RuntimeException("structure mapping with no structure or mapping");
}
String name = getMathSymbol(scm.getSpeciesContext(), sm.getGeometryClass());
Expression orig = rr.getRateRuleExpression();
Expression ex = getIdentifierSubstitutions(orig, scm.getSpeciesContext().getUnitDefinition(), sm.getGeometryClass());
GeometryClass gc = sm.getGeometryClass();
Domain domain = null;
if (gc != null) {
domain = new Domain(gc);
}
if (gc instanceof SurfaceClass) {
scm.setVariable(new MemVariable(scm.getSpeciesContext().getName(), domain));
} else {
scm.setVariable(new VolVariable(scm.getSpeciesContext().getName(), domain));
}
Variable oldVariablre = varHash.getVariable(name);
if (oldVariablre != null) {
// should always be null
varHash.removeVariable(name);
}
varHash.addVariable(scm.getVariable());
// // create the rate parameter
SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(sc);
SpeciesContextSpecParameter scsInitParam = scs.getInitialConditionParameter();
VCUnitDefinition scsInitParamUnit = scsInitParam.getUnitDefinition();
RateRuleRateParameter rateParam = null;
try {
Expression origExp = simContext.getRateRule(sc).getRateRuleExpression();
VCUnitDefinition rateUnit = modelUnitSystem.getInstance_TBD();
if (scsInitParamUnit != null && !scsInitParamUnit.equals(modelUnitSystem.getInstance_TBD())) {
rateUnit = scsInitParamUnit;
}
Expression rateExpr = getIdentifierSubstitutions(origExp, rateUnit, gc);
String argName = sc.getName() + MATH_FUNC_SUFFIX_RATERULE_RATE;
Variable param = newFunctionOrConstant(argName, rateExpr, gc);
varHash.addVariable(param);
rateParam = addRateRuleRateParameter(sc, rateExpr, PARAMETER_ROLE_RATERULE_RATE, rateUnit);
} catch (PropertyVetoException e) {
e.printStackTrace(System.out);
throw new MappingException(e.getMessage());
}
// we generate the ODE equation elsewhere (later)
rateRuleRateParamHash.put(scm.getVariable(), rateParam);
}
} else if (var != null && exp == null) {
// could be an event variable AND a rate rule variable - in which case we need a rate parameter for the event ODE equation
SpeciesContext sc = scm.getSpeciesContext();
boolean isRateRuleVar = rateRuleVarTargets.contains(sc);
boolean isEventAssignVar = eventAssignTargets.contains(sc);
if (isRateRuleVar && isEventAssignVar) {
// is both, so we make a rate parameter, like above
SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(sc);
SpeciesContextSpecParameter scsInitParam = scs.getInitialConditionParameter();
VCUnitDefinition scsInitParamUnit = scsInitParam.getUnitDefinition();
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(sc.getStructure());
GeometryClass gc = sm.getGeometryClass();
RateRuleRateParameter rateParam = null;
try {
Expression origExp = simContext.getRateRule(sc).getRateRuleExpression();
VCUnitDefinition rateUnit = modelUnitSystem.getInstance_TBD();
if (scsInitParamUnit != null && !scsInitParamUnit.equals(modelUnitSystem.getInstance_TBD())) {
rateUnit = scsInitParamUnit;
}
Expression rateExpr = getIdentifierSubstitutions(origExp, rateUnit, gc);
String argName = sc.getName() + MATH_FUNC_SUFFIX_RATERULE_RATE;
Variable param = newFunctionOrConstant(argName, rateExpr, gc);
varHash.addVariable(param);
rateParam = addRateRuleRateParameter(sc, rateExpr, PARAMETER_ROLE_RATERULE_RATE, rateUnit);
} catch (PropertyVetoException e) {
e.printStackTrace(System.out);
throw new MappingException(e.getMessage());
}
// we generate the ODE equation elsewhere (later)
rateRuleRateParamHash.put(var, rateParam);
}
}
}
//
for (ModelParameter mp : modelParameters) {
// global parameter used as rate rule variable
Variable var = varHash.getVariable(mp.getName());
RateRule rr = simContext.getRateRule(mp);
Expression modelParamExpr = mp.getExpression();
if (var == null && rr != null) {
// at this point var should be a constant
// we're under the assumption that it's non-spatial
GeometryClass[] geometryClasses = simContext.getGeometryContext().getGeometry().getGeometryClasses();
GeometryClass gc = geometryClasses[0];
// SubDomain subDomain = mathDesc.getSubDomains().nextElement();
// GeometryClass gc = getDefaultGeometryClass(modelParamExpr);
Domain domain = null;
if (gc != null) {
domain = new Domain(gc);
}
Variable variable;
if (gc instanceof SurfaceClass) {
variable = new MemVariable(mp.getName(), domain);
} else {
variable = new VolVariable(mp.getName(), domain);
}
varHash.addVariable(variable);
RateRuleRateParameter rateParam = null;
try {
Expression origExp = rr.getRateRuleExpression();
VCUnitDefinition rateUnit = modelUnitSystem.getInstance_TBD();
if (mp.getUnitDefinition() != null && !mp.getUnitDefinition().equals(modelUnitSystem.getInstance_TBD())) {
rateUnit = mp.getUnitDefinition().divideBy(timeUnit);
}
Expression rateExpr = getIdentifierSubstitutions(origExp, rateUnit, gc);
String argName = mp.getName() + MATH_FUNC_SUFFIX_RATERULE_RATE;
Variable param = newFunctionOrConstant(argName, rateExpr, gc);
varHash.addVariable(param);
rateParam = addRateRuleRateParameter(mp, rateExpr, PARAMETER_ROLE_RATERULE_RATE, rateUnit);
} catch (PropertyVetoException e) {
e.printStackTrace(System.out);
throw new MappingException(e.getMessage());
}
// no need to put it in the hash, we make the ODE Equation right here
// rateRuleRateParamHash.put(variable, rateParam);
// we know it's non-spatial
SubDomain subDomain = mathDesc.getSubDomains().nextElement();
Equation equation = null;
// TODO: replace the expression with the variable ex: "g0_protocol_init" computed above
Expression initial = new Expression(mp.getExpression());
// TODO: can it be null? should check and maybe try mp.getConstantValue() too ???
Variable mpInitVariable = initModelParameterHash.get(mp);
if (mpInitVariable != null) {
initial = new Expression(mpInitVariable.getName());
}
Expression rateExpr = new Expression(0.0);
// RateRuleRateParameter rateParam = rateRuleRateParamHash.get(variable);
if (rateParam != null) {
// ex: g0_rate
rateExpr = new Expression(getMathSymbol(rateParam, gc));
}
// ODE Equation for rate rule variable being a global parameter
equation = new OdeEquation(variable, initial, rateExpr);
subDomain.addEquation(equation);
}
}
//
// deal with assignment rules
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
// species context used as assignment rule variable
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() == null && scm.getDependencyExpression() != null) {
AssignmentRule ar = simContext.getAssignmentRule(scm.getSpeciesContext());
if (ar != null && (ar.getAssignmentRuleVar() instanceof SpeciesContext)) {
// TODO: we limit assignment rules to SpeciesContext for now
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(scm.getSpeciesContext().getStructure());
if (sm.getGeometryClass() == null) {
Structure s = sm.getStructure();
if (s != null) {
throw new RuntimeException("unmapped structure " + s.getName());
}
throw new RuntimeException("structure mapping with no structure or mapping");
}
String name = getMathSymbol(scm.getSpeciesContext(), sm.getGeometryClass());
Expression orig = ar.getAssignmentRuleExpression();
Expression ex = getIdentifierSubstitutions(orig, scm.getSpeciesContext().getUnitDefinition(), sm.getGeometryClass());
GeometryClass gc = sm.getGeometryClass();
Variable dependentVariable = newFunctionOrConstant(name, ex, gc);
dependentVariable.setDomain(new Domain(sm.getGeometryClass()));
varHash.removeVariable(name);
varHash.addVariable(dependentVariable);
}
}
}
for (ModelParameter mp : modelParameters) {
// global parameter used as assignment rule variable
Variable var = varHash.getVariable(mp.getName());
AssignmentRule ar = simContext.getAssignmentRule(mp);
Expression modelParamExpr = mp.getExpression();
if (var == null && ar != null) {
// at this point var (global parameter used as assignment rule variable) should be null
// we're under the assumption that it's non-spatial
GeometryClass[] geometryClasses = simContext.getGeometryContext().getGeometry().getGeometryClasses();
GeometryClass gc = geometryClasses[0];
SubDomain subDomain = mathDesc.getSubDomains().nextElement();
Expression origExp = ar.getAssignmentRuleExpression();
VCUnitDefinition rateUnit = modelUnitSystem.getInstance_TBD();
if (mp.getUnitDefinition() != null && !mp.getUnitDefinition().equals(modelUnitSystem.getInstance_TBD())) {
rateUnit = mp.getUnitDefinition();
}
Expression rateExpr = getIdentifierSubstitutions(origExp, rateUnit, gc);
String argName = mp.getName();
Variable param = newFunctionOrConstant(argName, rateExpr, gc);
varHash.addVariable(param);
}
}
//
// set Variables to MathDescription all at once with the order resolved by "VariableHash"
//
mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
//
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");
}
//
for (CompartmentSubdomainContext compartmentSubDomainContext : compartmentSubdomainContexts) {
SubVolume subVolume = compartmentSubDomainContext.subvolume;
CompartmentSubDomain subDomain = mathDesc.getCompartmentSubDomain(subVolume.getName());
//
// assign boundary condition types
//
StructureMapping[] mappedSMs = simContext.getGeometryContext().getStructureMappings(subVolume);
FeatureMapping mappedFM = null;
for (int i = 0; i < mappedSMs.length; i++) {
if (mappedSMs[i] instanceof FeatureMapping) {
if (mappedFM != null) {
lg.warn("WARNING:::: MathMapping.refreshMathDescription() ... assigning boundary condition types not unique");
}
mappedFM = (FeatureMapping) mappedSMs[i];
}
}
if (mappedFM != null) {
if (simContext.getGeometry().getDimension() > 0) {
subDomain.setBoundaryConditionXm(mappedFM.getBoundaryConditionTypeXm());
subDomain.setBoundaryConditionXp(mappedFM.getBoundaryConditionTypeXp());
}
if (simContext.getGeometry().getDimension() > 1) {
subDomain.setBoundaryConditionYm(mappedFM.getBoundaryConditionTypeYm());
subDomain.setBoundaryConditionYp(mappedFM.getBoundaryConditionTypeYp());
}
if (simContext.getGeometry().getDimension() > 2) {
subDomain.setBoundaryConditionZm(mappedFM.getBoundaryConditionTypeZm());
subDomain.setBoundaryConditionZp(mappedFM.getBoundaryConditionTypeZp());
}
}
//
// create equations
//
VolumeStructureAnalyzer structureAnalyzer = getVolumeStructureAnalyzer(subVolume);
Enumeration<SpeciesContextMapping> enumSCM = getSpeciesContextMappings();
while (enumSCM.hasMoreElements()) {
SpeciesContextMapping scm = enumSCM.nextElement();
SpeciesContext sc = scm.getSpeciesContext();
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(sc.getStructure());
SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(sc);
//
// if an independent volume variable, then create equation for it (if mapped to this subDomain)
//
final GeometryClass gc = sm.getGeometryClass();
if (gc == null || !gc.getName().equals(subDomain.getName())) {
continue;
}
SpeciesContextSpecParameter initConcParameter = scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
if ((scm.getVariable() instanceof VolumeRegionVariable) && scm.getDependencyExpression() == null) {
VolumeRegionVariable volumeRegionVariable = (VolumeRegionVariable) scm.getVariable();
Expression initial = getIdentifierSubstitutions(new Expression(initConcParameter, getNameScope()), initConcParameter.getUnitDefinition(), sm.getGeometryClass());
Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass());
VolumeRegionEquation volumeRegionEquation = new VolumeRegionEquation(volumeRegionVariable, initial);
volumeRegionEquation.setVolumeRateExpression(rate);
subDomain.addEquation(volumeRegionEquation);
} else if (scm.getVariable() instanceof VolVariable && scm.getDependencyExpression() == null) {
VolVariable variable = (VolVariable) scm.getVariable();
Equation equation = null;
if (sm.getGeometryClass() == subVolume) {
if (scm.isPDERequired()) {
//
// species context belongs to this subDomain
//
Expression initial = getIdentifierSubstitutions(new Expression(initConcParameter, getNameScope()), initConcParameter.getUnitDefinition(), sm.getGeometryClass());
Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass());
SpeciesContextSpecParameter diffusionParameter = scs.getDiffusionParameter();
Expression diffusion = getIdentifierSubstitutions(new Expression(diffusionParameter, getNameScope()), diffusionParameter.getUnitDefinition(), sm.getGeometryClass());
equation = new PdeEquation(variable, initial, rate, diffusion);
((PdeEquation) equation).setBoundaryXm((scs.getBoundaryXmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXmParameter(), sm.getGeometryClass())));
((PdeEquation) equation).setBoundaryXp((scs.getBoundaryXpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXpParameter(), sm.getGeometryClass())));
((PdeEquation) equation).setBoundaryYm((scs.getBoundaryYmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYmParameter(), sm.getGeometryClass())));
((PdeEquation) equation).setBoundaryYp((scs.getBoundaryYpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYpParameter(), sm.getGeometryClass())));
((PdeEquation) equation).setBoundaryZm((scs.getBoundaryZmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZmParameter(), sm.getGeometryClass())));
((PdeEquation) equation).setBoundaryZp((scs.getBoundaryZpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZpParameter(), sm.getGeometryClass())));
if (simContext.getGeometry().getDimension() >= 1) {
Expression velXExp = null;
if (scs.getVelocityXParameter().getExpression() != null) {
velXExp = new Expression(getMathSymbol(scs.getVelocityXParameter(), sm.getGeometryClass()));
} else {
SpatialQuantity[] velX_quantities = scs.getVelocityQuantities(QuantityComponent.X);
if (velX_quantities.length > 0) {
int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(subVolume).length;
if (velX_quantities.length == 1 && numRegions == 1) {
velXExp = new Expression(getMathSymbol(velX_quantities[0], sm.getGeometryClass()));
} else {
throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
}
}
}
((PdeEquation) equation).setVelocityX(velXExp);
}
if (simContext.getGeometry().getDimension() >= 2) {
Expression velYExp = null;
if (scs.getVelocityYParameter().getExpression() != null) {
velYExp = new Expression(getMathSymbol(scs.getVelocityYParameter(), sm.getGeometryClass()));
} else {
SpatialQuantity[] velY_quantities = scs.getVelocityQuantities(QuantityComponent.Y);
if (velY_quantities.length > 0) {
int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(subVolume).length;
if (velY_quantities.length == 1 && numRegions == 1) {
velYExp = new Expression(getMathSymbol(velY_quantities[0], sm.getGeometryClass()));
} else {
throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
}
}
}
((PdeEquation) equation).setVelocityY(velYExp);
}
if (simContext.getGeometry().getDimension() == 3) {
Expression velZExp = null;
if (scs.getVelocityZParameter().getExpression() != null) {
velZExp = new Expression(getMathSymbol(scs.getVelocityZParameter(), sm.getGeometryClass()));
} else {
SpatialQuantity[] velZ_quantities = scs.getVelocityQuantities(QuantityComponent.Z);
if (velZ_quantities.length > 0) {
int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(subVolume).length;
if (velZ_quantities.length == 1 && numRegions == 1) {
velZExp = new Expression(getMathSymbol(velZ_quantities[0], sm.getGeometryClass()));
} else {
throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
}
}
}
((PdeEquation) equation).setVelocityZ(velZExp);
}
subDomain.replaceEquation(equation);
} else {
//
// ODE - species context belongs to this subDomain
//
Expression initial = new Expression(getMathSymbol(initConcParameter, null));
Expression rate = (scm.getRate() == null) ? new Expression(0.0) : getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass());
//
// if it's an event assignment variable AND a rate rule variable
// we replace the event rate computed above (which should be zero) with the RateRuleParameter expression
//
RateRuleRateParameter rateParam = rateRuleRateParamHash.get(variable);
if (rateParam != null) {
rate = new Expression(getMathSymbol(rateParam, null));
}
equation = new OdeEquation(variable, initial, rate);
subDomain.replaceEquation(equation);
}
}
} else if (scm.getVariable() instanceof VolVariable && scm.getDependencyExpression() != null) {
// rate rule variables are like this
RateRule rr = simContext.getRateRule(scm.getSpeciesContext());
if (rr != null && (rr.getRateRuleVar() instanceof SpeciesContext)) {
//
// we generate rate rule ODE equation only for species variable that are NOT event assignment variable (see right above)
// for global parameters variable we do it elsewhere
//
VolVariable variable = (VolVariable) scm.getVariable();
Equation equation = null;
if (sm.getGeometryClass() == subVolume) {
Expression initial = new Expression(getMathSymbol(initConcParameter, null));
Expression rateExpr = new Expression(0.0);
RateRuleRateParameter rateParam = rateRuleRateParamHash.get(variable);
if (rateParam != null) {
rateExpr = new Expression(getMathSymbol(rateParam, null));
}
equation = new OdeEquation(variable, initial, rateExpr);
subDomain.addEquation(equation);
}
}
}
}
//
// create fast system (if neccessary)
//
SpeciesContextMapping[] fastSpeciesContextMappings = structureAnalyzer.getFastSpeciesContextMappings();
if (fastSpeciesContextMappings != null) {
FastSystem fastSystem = new FastSystem(mathDesc);
for (int i = 0; i < fastSpeciesContextMappings.length; i++) {
SpeciesContextMapping scm = fastSpeciesContextMappings[i];
if (scm.getFastInvariant() == null) {
//
// independant-fast variable, create a fastRate object
//
Expression rate = getIdentifierSubstitutions(scm.getFastRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), subVolume);
FastRate fastRate = new FastRate(rate);
fastSystem.addFastRate(fastRate);
} else {
//
// dependant-fast variable, create a fastInvariant object
//
Expression rate = getIdentifierSubstitutions(scm.getFastInvariant(), modelUnitSystem.getVolumeConcentrationUnit(), subVolume);
FastInvariant fastInvariant = new FastInvariant(rate);
fastSystem.addFastInvariant(fastInvariant);
}
}
subDomain.setFastSystem(fastSystem);
// constructor calls the 'refresh' method which constructs depemdency matrix, dependent/independent vars and pseudoconstants, etc.
// FastSystemAnalyzer fs_analyzer =
new FastSystemAnalyzer(fastSystem, mathDesc);
}
//
// create ode's for voltages to be calculated on unresolved membranes mapped to this subVolume
//
Structure[] localStructures = simContext.getGeometryContext().getStructuresFromGeometryClass(subVolume);
for (int sIndex = 0; sIndex < localStructures.length; sIndex++) {
if (localStructures[sIndex] instanceof Membrane) {
Membrane membrane = (Membrane) localStructures[sIndex];
MembraneMapping membraneMapping = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(membrane);
if ((membraneMapping.getGeometryClass() instanceof SubVolume) && membraneMapping.getCalculateVoltage()) {
MembraneElectricalDevice capacitiveDevice = potentialMapping.getCapacitiveDevice(membrane);
if (capacitiveDevice.getDependentVoltageExpression() == null) {
VolVariable vVar = (VolVariable) mathDesc.getVariable(getMathSymbol(capacitiveDevice.getVoltageSymbol(), membraneMapping.getGeometryClass()));
Expression initExp = new Expression(getMathSymbol(capacitiveDevice.getMembraneMapping().getInitialVoltageParameter(), membraneMapping.getGeometryClass()));
subDomain.addEquation(new OdeEquation(vVar, initExp, getIdentifierSubstitutions(potentialMapping.getOdeRHS(capacitiveDevice, this), membrane.getMembraneVoltage().getUnitDefinition().divideBy(timeUnit), membraneMapping.getGeometryClass())));
} else {
//
//
//
}
}
}
}
}
//
for (MembraneSubdomainContext memSubdomainContext : membraneSubdomainContexts) {
MembraneSubDomain memSubDomain = memSubdomainContext.membraneSubdomain;
SurfaceClass surfaceClass = memSubdomainContext.surfaceClass;
for (SurfaceRegionObject surfaceRegionObject : memSubdomainContext.surfaceRegionObjects) {
if (surfaceRegionObject.isQuantityCategoryEnabled(QuantityCategory.SurfaceVelocity)) {
int dim = simContext.getGeometry().getDimension();
if (dim != 2) {
throw new MappingException("Membrane Velocity only supported for 2D geometries");
}
if (simContext.getGeometry().getDimension() >= 1) {
SpatialQuantity velXQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.SurfaceVelocity, QuantityComponent.X);
Expression velXExp = new Expression(velXQuantity, simContext.getNameScope());
memSubDomain.setVelocityX(getIdentifierSubstitutions(velXExp, velXQuantity.getUnitDefinition(), surfaceClass));
}
if (simContext.getGeometry().getDimension() >= 2) {
SpatialQuantity velYQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.SurfaceVelocity, QuantityComponent.Y);
Expression velYExp = new Expression(velYQuantity, simContext.getNameScope());
memSubDomain.setVelocityY(getIdentifierSubstitutions(velYExp, velYQuantity.getUnitDefinition(), surfaceClass));
}
if (simContext.getGeometry().getDimension() == 3) {
SpatialQuantity velZQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.SurfaceVelocity, QuantityComponent.Z);
Expression velZExp = new Expression(velZQuantity, simContext.getNameScope());
// memSubDomain.setVelocityZ(getIdentifierSubstitutions(velZExp, velZQuantity.getUnitDefinition(), surfaceClass));
throw new MappingException("Membrane Velocity not supported for 2D problems");
}
}
}
//
// create equations for membrane-bound molecular species
//
MembraneStructureAnalyzer membraneStructureAnalyzer = getMembraneStructureAnalyzer(surfaceClass);
Enumeration<SpeciesContextMapping> enumSCM = getSpeciesContextMappings();
while (enumSCM.hasMoreElements()) {
SpeciesContextMapping scm = enumSCM.nextElement();
SpeciesContext sc = scm.getSpeciesContext();
SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(sc);
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(sc.getStructure());
//
if ((scm.getVariable() instanceof MembraneRegionVariable) && scm.getDependencyExpression() == null) {
MembraneRegionEquation equation = null;
MembraneRegionVariable memRegionVar = (MembraneRegionVariable) scm.getVariable();
if (sm.getGeometryClass() == surfaceClass) {
//
// species context belongs to this subDomain
//
Expression initial = new Expression(getMathSymbol(scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration), sm.getGeometryClass()));
Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass());
equation = new MembraneRegionEquation(memRegionVar, initial);
equation.setMembraneRateExpression(rate);
// equation.setUniformRateExpression(newUniformRateExpression);
memSubDomain.replaceEquation(equation);
}
} else if ((scm.getVariable() instanceof MemVariable) && scm.getDependencyExpression() == null) {
//
if (sm.getGeometryClass() == surfaceClass) {
Equation equation = null;
MemVariable variable = (MemVariable) scm.getVariable();
if (scm.isPDERequired()) {
//
// PDE
//
//
// species context belongs to this subDomain
//
Expression initial = new Expression(getMathSymbol(scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration), sm.getGeometryClass()));
Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass());
Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), sm.getGeometryClass()));
equation = new PdeEquation(variable, initial, rate, diffusion);
((PdeEquation) equation).setBoundaryXm((scs.getBoundaryXmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXmParameter(), sm.getGeometryClass())));
((PdeEquation) equation).setBoundaryXp((scs.getBoundaryXpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXpParameter(), sm.getGeometryClass())));
((PdeEquation) equation).setBoundaryYm((scs.getBoundaryYmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYmParameter(), sm.getGeometryClass())));
((PdeEquation) equation).setBoundaryYp((scs.getBoundaryYpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYpParameter(), sm.getGeometryClass())));
((PdeEquation) equation).setBoundaryZm((scs.getBoundaryZmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZmParameter(), sm.getGeometryClass())));
((PdeEquation) equation).setBoundaryZp((scs.getBoundaryZpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZpParameter(), sm.getGeometryClass())));
memSubDomain.replaceEquation(equation);
} else {
//
// ODE
//
//
// species context belongs to this subDomain
//
Expression initial = new Expression(getMathSymbol(scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration), null));
Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass());
equation = new OdeEquation(variable, initial, rate);
memSubDomain.replaceEquation(equation);
}
}
}
}
Enumeration<SpeciesContextMapping> enum_scm = getSpeciesContextMappings();
while (enum_scm.hasMoreElements()) {
SpeciesContextMapping scm = enum_scm.nextElement();
if (scm.isPDERequired() || scm.getVariable() instanceof VolumeRegionVariable) {
// Species species = scm.getSpeciesContext().getSpecies();
Variable var = scm.getVariable();
final Domain dm = var.getDomain();
if (dm != null) {
final String domainName = dm.getName();
if (sameName(domainName, memSubDomain.getInsideCompartment()) || sameName(domainName, memSubDomain.getOutsideCompartment())) {
JumpCondition jc = memSubDomain.getJumpCondition(var);
if (jc == null) {
// System.out.println("MathMapping.refreshMathDescription(), adding jump condition for diffusing variable "+var.getName()+" on membrane "+membraneStructureAnalyzer.getMembrane().getName());
if (var instanceof VolVariable) {
jc = new JumpCondition((VolVariable) var);
} else if (var instanceof VolumeRegionVariable) {
jc = new JumpCondition((VolumeRegionVariable) var);
} else {
throw new RuntimeException("unexpected Variable type " + var.getClass().getName());
}
memSubDomain.addJumpCondition(jc);
}
}
}
}
}
//
// set jump conditions for any volume variables or volume region variables that have explicitly defined fluxes
//
ResolvedFlux[] resolvedFluxes = membraneStructureAnalyzer.getResolvedFluxes();
if (resolvedFluxes != null) {
for (int i = 0; i < resolvedFluxes.length; i++) {
SpeciesContext sc = resolvedFluxes[i].getSpeciesContext();
SpeciesContextMapping scm = getSpeciesContextMapping(sc);
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure());
if (scm.getVariable() instanceof VolVariable && scm.isPDERequired()) {
VolVariable volVar = (VolVariable) scm.getVariable();
JumpCondition jc = memSubDomain.getJumpCondition(volVar);
if (jc == null) {
jc = new JumpCondition(volVar);
memSubDomain.addJumpCondition(jc);
}
Expression flux = getIdentifierSubstitutions(resolvedFluxes[i].getFluxExpression(), resolvedFluxes[i].getUnitDefinition(), membraneStructureAnalyzer.getSurfaceClass());
if (memSubDomain.getInsideCompartment().getName().equals(sm.getGeometryClass().getName())) {
jc.setInFlux(flux);
} else if (memSubDomain.getOutsideCompartment().getName().equals(sm.getGeometryClass().getName())) {
jc.setOutFlux(flux);
} else {
throw new RuntimeException("Application " + simContext.getName() + " : " + scm.getSpeciesContext().getName() + " has spatially resolved flux at membrane " + scm.getSpeciesContext().getStructure().getName() + " with a non-local flux species " + scm.getSpeciesContext().getName());
}
} else if (scm.getVariable() instanceof VolumeRegionVariable) {
VolumeRegionVariable volRegionVar = (VolumeRegionVariable) scm.getVariable();
JumpCondition jc = memSubDomain.getJumpCondition(volRegionVar);
if (jc == null) {
jc = new JumpCondition(volRegionVar);
memSubDomain.addJumpCondition(jc);
}
Expression flux = getIdentifierSubstitutions(resolvedFluxes[i].getFluxExpression(), resolvedFluxes[i].getUnitDefinition(), membraneStructureAnalyzer.getSurfaceClass());
if (memSubDomain.getInsideCompartment().getName().equals(sm.getGeometryClass().getName())) {
jc.setInFlux(flux);
} else if (memSubDomain.getOutsideCompartment().getName().equals(sm.getGeometryClass().getName())) {
jc.setOutFlux(flux);
} else {
throw new RuntimeException("Application " + simContext.getName() + " : " + scm.getSpeciesContext().getName() + " has spatially resolved flux at membrane " + scm.getSpeciesContext().getStructure().getName() + " with a non-local flux species " + scm.getSpeciesContext().getName());
}
} else {
throw new MappingException("Application " + simContext.getName() + " : " + scm.getSpeciesContext().getName() + " has spatially resolved flux at membrane " + scm.getSpeciesContext().getStructure().getName() + ", but doesn't diffuse in compartment " + scm.getSpeciesContext().getStructure().getName());
}
}
}
//
// create fast system (if neccessary)
//
SpeciesContextMapping[] fastSpeciesContextMappings = membraneStructureAnalyzer.getFastSpeciesContextMappings();
if (fastSpeciesContextMappings != null) {
FastSystem fastSystem = new FastSystem(mathDesc);
for (int i = 0; i < fastSpeciesContextMappings.length; i++) {
SpeciesContextMapping scm = fastSpeciesContextMappings[i];
if (scm.getFastInvariant() == null) {
//
// independant-fast variable, create a fastRate object
//
VCUnitDefinition rateUnit = scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit);
FastRate fastRate = new FastRate(getIdentifierSubstitutions(scm.getFastRate(), rateUnit, surfaceClass));
fastSystem.addFastRate(fastRate);
} else {
//
// dependant-fast variable, create a fastInvariant object
//
VCUnitDefinition invariantUnit = scm.getSpeciesContext().getUnitDefinition();
FastInvariant fastInvariant = new FastInvariant(getIdentifierSubstitutions(scm.getFastInvariant(), invariantUnit, surfaceClass));
fastSystem.addFastInvariant(fastInvariant);
}
}
memSubDomain.setFastSystem(fastSystem);
// constructor calls the 'refresh' method which constructs depemdency matrix, dependent/independent vars and pseudoconstants, etc.
// FastSystemAnalyzer fs_analyzer =
new FastSystemAnalyzer(fastSystem, mathDesc);
}
//
// create Membrane-region equations for potential of this resolved membrane
//
Structure[] resolvedSurfaceStructures = membraneStructureAnalyzer.getStructures();
for (int m = 0; m < resolvedSurfaceStructures.length; m++) {
if (resolvedSurfaceStructures[m] instanceof Membrane) {
Membrane membrane = (Membrane) resolvedSurfaceStructures[m];
MembraneMapping membraneMapping = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(membrane);
if (membraneMapping.getCalculateVoltage()) {
ElectricalDevice[] membraneDevices = potentialMapping.getElectricalDevices(membrane);
int numCapacitiveDevices = 0;
MembraneElectricalDevice capacitiveDevice = null;
for (int i = 0; i < membraneDevices.length; i++) {
if (membraneDevices[i] instanceof MembraneElectricalDevice) {
numCapacitiveDevices++;
capacitiveDevice = (MembraneElectricalDevice) membraneDevices[i];
}
}
if (numCapacitiveDevices != 1) {
throw new MappingException("expecting 1 capacitive electrical device on graph edge for membrane " + membrane.getName() + ", found '" + numCapacitiveDevices + "'");
}
if (mathDesc.getVariable(getMathSymbol(capacitiveDevice.getVoltageSymbol(), membraneMapping.getGeometryClass())) instanceof MembraneRegionVariable) {
MembraneRegionVariable vVar = (MembraneRegionVariable) mathDesc.getVariable(getMathSymbol(capacitiveDevice.getVoltageSymbol(), membraneMapping.getGeometryClass()));
Parameter initialVoltageParm = capacitiveDevice.getMembraneMapping().getInitialVoltageParameter();
Expression initExp = getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), capacitiveDevice.getMembraneMapping().getGeometryClass());
MembraneRegionEquation vEquation = new MembraneRegionEquation(vVar, initExp);
vEquation.setMembraneRateExpression(getIdentifierSubstitutions(potentialMapping.getOdeRHS(capacitiveDevice, this), membrane.getMembraneVoltage().getUnitDefinition().divideBy(timeUnit), capacitiveDevice.getMembraneMapping().getGeometryClass()));
memSubDomain.addEquation(vEquation);
}
}
}
}
}
// create equations for event assignment or rate rule targets that are model params/species, etc.
Set<VolVariable> hashKeySet = eventVolVarHash.keySet();
Iterator<VolVariable> volVarsIter = hashKeySet.iterator();
// working under the assumption that we are dealing with non-spatial math, hence only one compartment domain!
SubDomain subDomain = mathDesc.getSubDomains().nextElement();
while (volVarsIter.hasNext()) {
VolVariable volVar = volVarsIter.next();
EventAssignmentOrRateRuleInitParameter initParam = eventVolVarHash.get(volVar);
// check event initial condition, it shouldn't contain vars, we have to do it here, coz we want to substitute functions...etc.
Expression eapExp = MathUtilities.substituteFunctions(initParam.getExpression(), mathDesc);
if (eapExp.getSymbols() != null) {
for (String symbol : eapExp.getSymbols()) {
SymbolTableEntry ste = eapExp.getSymbolBinding(symbol);
if (ste instanceof VolVariable || ste instanceof MemVariable) {
throw new MathException("Variables are not allowed in Event assignment initial condition.\nEvent assignment target: " + volVar.getName() + " has variable (" + symbol + ") in its expression.");
}
}
}
Expression rateExpr = new Expression(0.0);
RateRuleRateParameter rateParam = rateRuleRateParamHash.get(volVar);
if (rateParam != null) {
// this is a rate rule, get its expression.
rateExpr = new Expression(getMathSymbol(rateParam, null));
}
Equation equation = new OdeEquation(volVar, new Expression(getMathSymbol(initParam, null)), rateExpr);
subDomain.addEquation(equation);
}
// events - add events to math desc for event assignments that have parameters as target variables
if (bioevents != null && bioevents.length > 0) {
for (BioEvent be : bioevents) {
// transform the bioEvent trigger/delay to math Event
LocalParameter genTriggerParam = be.getParameter(BioEventParameterType.GeneralTriggerFunction);
Expression mathTriggerExpr = getIdentifierSubstitutions(new Expression(genTriggerParam, be.getNameScope()), modelUnitSystem.getInstance_DIMENSIONLESS(), null);
Delay mathDelay = null;
LocalParameter delayParam = be.getParameter(BioEventParameterType.TriggerDelay);
if (delayParam != null && delayParam.getExpression() != null && !delayParam.getExpression().compareEqual(new Expression(0.0))) {
boolean bUseValsFromTriggerTime = be.getUseValuesFromTriggerTime();
Expression mathDelayExpr = getIdentifierSubstitutions(new Expression(delayParam, be.getNameScope()), timeUnit, null);
mathDelay = new Delay(bUseValsFromTriggerTime, mathDelayExpr);
}
// now deal with (bio)event Assignment translation to math EventAssignment
ArrayList<EventAssignment> eventAssignments = be.getEventAssignments();
ArrayList<Event.EventAssignment> mathEventAssignmentsList = new ArrayList<Event.EventAssignment>();
if (eventAssignments != null) {
for (EventAssignment ea : eventAssignments) {
SymbolTableEntry ste = simContext.getEntry(ea.getTarget().getName());
if (ste instanceof StructureSize) {
throw new RuntimeException("Event Assignment Variable for compartment size is not supported yet");
}
VCUnitDefinition eventAssignVarUnit = ste.getUnitDefinition();
Variable variable = varHash.getVariable(ste.getName());
Event.EventAssignment mathEA = new Event.EventAssignment(variable, getIdentifierSubstitutions(ea.getAssignmentExpression(), eventAssignVarUnit, null));
mathEventAssignmentsList.add(mathEA);
}
}
// use the translated trigger, delay and event assignments to create (math) event
Event mathEvent = new Event(be.getName(), mathTriggerExpr, mathDelay, mathEventAssignmentsList);
mathDesc.addEvent(mathEvent);
}
}
if (simContext.getMicroscopeMeasurement() != null && simContext.getMicroscopeMeasurement().getFluorescentSpecies().size() > 0) {
MicroscopeMeasurement measurement = simContext.getMicroscopeMeasurement();
Expression volumeConcExp = new Expression(0.0);
Expression membraneDensityExp = new Expression(0.0);
for (SpeciesContext speciesContext : measurement.getFluorescentSpecies()) {
GeometryClass geometryClass = simContext.getGeometryContext().getStructureMapping(speciesContext.getStructure()).getGeometryClass();
StructureMapping structureMapping = simContext.getGeometryContext().getStructureMapping(speciesContext.getStructure());
StructureMappingParameter unitSizeParameter = structureMapping.getUnitSizeParameter();
Expression mappedSpeciesContextExpression = Expression.mult(unitSizeParameter.getExpression(), new Expression(getMathSymbol(speciesContext, geometryClass)));
VCUnitDefinition mappedSpeciesContextUnit = unitSizeParameter.getUnitDefinition().multiplyBy(speciesContext.getUnitDefinition());
if (geometryClass instanceof SubVolume) {
// volume function
int dimension = 3;
VCUnitDefinition desiredConcUnits = model.getUnitSystem().getInstance("molecules").divideBy(model.getUnitSystem().getLengthUnit().raiseTo(new ucar.units_vcell.RationalNumber(dimension)));
Expression unitFactor = getUnitFactor(desiredConcUnits.divideBy(mappedSpeciesContextUnit));
volumeConcExp = Expression.add(volumeConcExp, Expression.mult(unitFactor, mappedSpeciesContextExpression)).flatten();
} else if (geometryClass instanceof SurfaceClass) {
// membrane function
int dimension = 2;
VCUnitDefinition desiredSurfaceDensityUnits = model.getUnitSystem().getInstance("molecules").divideBy(model.getUnitSystem().getLengthUnit().raiseTo(new ucar.units_vcell.RationalNumber(dimension)));
Expression unitFactor = getUnitFactor(desiredSurfaceDensityUnits.divideBy(mappedSpeciesContextUnit));
membraneDensityExp = Expression.add(membraneDensityExp, Expression.mult(unitFactor, mappedSpeciesContextExpression)).flatten();
} else {
throw new MathException("unsupported geometry mapping for microscopy measurement");
}
}
ConvolutionKernel kernel = measurement.getConvolutionKernel();
if (kernel instanceof ExperimentalPSF) {
if (!membraneDensityExp.isZero()) {
throw new MappingException("membrane variables and functions not yet supported for Z projection in Microcopy Measurements");
}
ExperimentalPSF psf = (ExperimentalPSF) kernel;
DataSymbol psfDataSymbol = psf.getPSFDataSymbol();
if (psfDataSymbol instanceof FieldDataSymbol) {
FieldDataSymbol fieldDataSymbol = (FieldDataSymbol) psfDataSymbol;
String fieldDataName = ((FieldDataSymbol) psfDataSymbol).getExternalDataIdentifier().getName();
Expression psfExp = Expression.function(FieldFunctionDefinition.FUNCTION_name, new Expression("'" + fieldDataName + "'"), new Expression("'" + fieldDataSymbol.getFieldDataVarName() + "'"), new Expression(fieldDataSymbol.getFieldDataVarTime()), new Expression("'" + fieldDataSymbol.getFieldDataVarType() + "'"));
varHash.addVariable(new Function("__PSF__", psfExp, null));
}
Expression convExp = Expression.function(ConvFunctionDefinition.FUNCTION_name, volumeConcExp, new Expression("__PSF__"));
varHash.addVariable(newFunctionOrConstant(measurement.getName(), convExp, null));
} else if (kernel instanceof GaussianConvolutionKernel) {
GaussianConvolutionKernel gaussianConvolutionKernel = (GaussianConvolutionKernel) kernel;
GaussianConvolutionDataGeneratorKernel mathKernel = new GaussianConvolutionDataGeneratorKernel(gaussianConvolutionKernel.getSigmaXY_um(), gaussianConvolutionKernel.getSigmaZ_um());
ConvolutionDataGenerator dataGenerator = new ConvolutionDataGenerator(measurement.getName(), mathKernel, volumeConcExp, membraneDensityExp);
mathDesc.getPostProcessingBlock().addDataGenerator(dataGenerator);
} else if (kernel instanceof ProjectionZKernel) {
if (mathDesc.getGeometry().getDimension() == 3) {
if (!membraneDensityExp.isZero()) {
throw new MappingException("membrane variables and functions not yet supported for Z projection in Microcopy Measurements");
}
ProjectionDataGenerator dataGenerator = new ProjectionDataGenerator(measurement.getName(), null, ProjectionDataGenerator.Axis.z, ProjectionDataGenerator.Operation.sum, volumeConcExp);
mathDesc.getPostProcessingBlock().addDataGenerator(dataGenerator);
} else {
throw new MappingException("Z Projection is only supported in 3D spatial applications.");
}
}
}
//
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
Variable variable = newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass());
if (mathDesc.getVariable(variable.getName()) == null) {
mathDesc.addVariable(variable);
}
}
}
if (!mathDesc.isValid()) {
System.out.println(mathDesc.getVCML_database());
throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
}
// System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string begin ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
// System.out.println(mathDesc.getVCML());
// System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string end ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
}
Also used :
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
ArrayList(java.util.ArrayList)
StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter)
SpeciesContext(cbit.vcell.model.SpeciesContext)
Feature(cbit.vcell.model.Feature)
MemVariable(cbit.vcell.math.MemVariable)
SubVolume(cbit.vcell.geometry.SubVolume)
Vector(java.util.Vector)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
FastInvariant(cbit.vcell.math.FastInvariant)
GaussianConvolutionKernel(cbit.vcell.mapping.MicroscopeMeasurement.GaussianConvolutionKernel)
PropertyVetoException(java.beans.PropertyVetoException)
FieldDataSymbol(cbit.vcell.data.FieldDataSymbol)
DataSymbol(cbit.vcell.data.DataSymbol)
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
FastSystem(cbit.vcell.math.FastSystem)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
ReactionStep(cbit.vcell.model.ReactionStep)
Map(java.util.Map)
HashMap(java.util.HashMap)
MembraneRegionEquation(cbit.vcell.math.MembraneRegionEquation)
SurfaceClass(cbit.vcell.geometry.SurfaceClass)
VariableHash(cbit.vcell.math.VariableHash)
ExperimentalPSF(cbit.vcell.mapping.MicroscopeMeasurement.ExperimentalPSF)
ConvolutionDataGenerator(cbit.vcell.math.ConvolutionDataGenerator)
GaussianConvolutionDataGeneratorKernel(cbit.vcell.math.ConvolutionDataGenerator.GaussianConvolutionDataGeneratorKernel)
Structure(cbit.vcell.model.Structure)
ModelUnitSystem(cbit.vcell.model.ModelUnitSystem)
LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter)
VoltageClampElectricalDevice(cbit.vcell.mapping.potential.VoltageClampElectricalDevice)
StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter)
Expression(cbit.vcell.parser.Expression)
Model(cbit.vcell.model.Model)
StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter)
Parameter(cbit.vcell.model.Parameter)
SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter)
SimulationContextParameter(cbit.vcell.mapping.SimulationContext.SimulationContextParameter)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter)
Event(cbit.vcell.math.Event)
ProjectionDataGenerator(cbit.vcell.math.ProjectionDataGenerator)
SurfaceRegionObject(cbit.vcell.mapping.spatial.SurfaceRegionObject)
FieldDataSymbol(cbit.vcell.data.FieldDataSymbol)
MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable)
HashMap(java.util.HashMap)
MathDescription(cbit.vcell.math.MathDescription)
MembraneElectricalDevice(cbit.vcell.mapping.potential.MembraneElectricalDevice)
Delay(cbit.vcell.math.Event.Delay)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
PointSubDomain(cbit.vcell.math.PointSubDomain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
PropertyVetoException(java.beans.PropertyVetoException)
PdeEquation(cbit.vcell.math.PdeEquation)
CurrentClampElectricalDevice(cbit.vcell.mapping.potential.CurrentClampElectricalDevice)
SpatialQuantity(cbit.vcell.mapping.spatial.SpatialObject.SpatialQuantity)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
EventAssignment(cbit.vcell.mapping.BioEvent.EventAssignment)
CurrentClampElectricalDevice(cbit.vcell.mapping.potential.CurrentClampElectricalDevice)
MembraneElectricalDevice(cbit.vcell.mapping.potential.MembraneElectricalDevice)
ElectricalDevice(cbit.vcell.mapping.potential.ElectricalDevice)
VoltageClampElectricalDevice(cbit.vcell.mapping.potential.VoltageClampElectricalDevice)
VolVariable(cbit.vcell.math.VolVariable)
StructureSize(cbit.vcell.model.Structure.StructureSize)
ProjectionZKernel(cbit.vcell.mapping.MicroscopeMeasurement.ProjectionZKernel)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
OdeEquation(cbit.vcell.math.OdeEquation)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
PointSubDomain(cbit.vcell.math.PointSubDomain)
Domain(cbit.vcell.math.Variable.Domain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
JumpCondition(cbit.vcell.math.JumpCondition)
GeometryClass(cbit.vcell.geometry.GeometryClass)
VolVariable(cbit.vcell.math.VolVariable)
MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable)
PointVariable(cbit.vcell.math.PointVariable)
VolumeRegionVariable(cbit.vcell.math.VolumeRegionVariable)
MemVariable(cbit.vcell.math.MemVariable)
Variable(cbit.vcell.math.Variable)
SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter)
Constant(cbit.vcell.math.Constant)
VolumeRegionVariable(cbit.vcell.math.VolumeRegionVariable)
Function(cbit.vcell.math.Function)
SymbolTableEntry(cbit.vcell.parser.SymbolTableEntry)
Membrane(cbit.vcell.model.Membrane)
VolumeRegionEquation(cbit.vcell.math.VolumeRegionEquation)
PotentialMapping(cbit.vcell.mapping.potential.PotentialMapping)
EventAssignment(cbit.vcell.mapping.BioEvent.EventAssignment)
FieldFunctionArguments(cbit.vcell.field.FieldFunctionArguments)
PdeEquation(cbit.vcell.math.PdeEquation)
ComputeMembraneMetricEquation(cbit.vcell.math.ComputeMembraneMetricEquation)
VolumeRegionEquation(cbit.vcell.math.VolumeRegionEquation)
OdeEquation(cbit.vcell.math.OdeEquation)
MembraneRegionEquation(cbit.vcell.math.MembraneRegionEquation)
Equation(cbit.vcell.math.Equation)
FastRate(cbit.vcell.math.FastRate)
SimulationContextParameter(cbit.vcell.mapping.SimulationContext.SimulationContextParameter)
ConvolutionKernel(cbit.vcell.mapping.MicroscopeMeasurement.ConvolutionKernel)
GaussianConvolutionKernel(cbit.vcell.mapping.MicroscopeMeasurement.GaussianConvolutionKernel)
MathException(cbit.vcell.math.MathException)
Example 3 with VariableHash
use of cbit.vcell.math.VariableHash in project vcell by virtualcell.
the class ParticleMathMapping method refreshMathDescription.
/**
* This method was created in VisualAge.
*/
private void refreshMathDescription() throws MappingException, MatrixException, MathException, ExpressionException, ModelException {
getSimulationContext().checkValidity();
if (getSimulationContext().getGeometry().getDimension() == 0) {
throw new MappingException("particle math mapping requires spatial geometry - dimension >= 1");
}
StructureMapping[] structureMappings = getSimulationContext().getGeometryContext().getStructureMappings();
for (int i = 0; i < structureMappings.length; i++) {
if (structureMappings[i] instanceof MembraneMapping) {
if (((MembraneMapping) structureMappings[i]).getCalculateVoltage()) {
throw new MappingException("electric potential not yet supported for particle models");
}
}
}
//
// fail if any events
//
BioEvent[] bioEvents = getSimulationContext().getBioEvents();
if (bioEvents != null && bioEvents.length > 0) {
throw new MappingException("events not yet supported for particle-based models");
}
//
// gather only those reactionSteps that are not "excluded"
//
ReactionSpec[] reactionSpecs = getSimulationContext().getReactionContext().getReactionSpecs();
Vector<ReactionStep> rsList = new Vector<ReactionStep>();
for (int i = 0; i < reactionSpecs.length; i++) {
if (reactionSpecs[i].isExcluded() == false) {
if (reactionSpecs[i].isFast()) {
throw new MappingException("fast reactions not supported for particle models");
}
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);
}
}
//
// temporarily place all variables in a hashtable (before binding) and discarding duplicates (check for equality)
//
VariableHash varHash = new VariableHash();
// //
// // verify that all structures are mapped to geometry classes and all geometry classes are mapped to a structure
// //
// Structure structures[] = getSimulationContext().getGeometryContext().getModel().getStructures();
// for (int i = 0; i < structures.length; i++){
// StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(structures[i]);
// if (sm==null || (sm.getGeometryClass() == null)){
// throw new MappingException("model structure '"+structures[i].getName()+"' not mapped to a geometry subdomain");
// }
// if (sm.getUnitSizeParameter()!=null){
// Expression unitSizeExp = sm.getUnitSizeParameter().getExpression();
// if(unitSizeExp != null)
// {
// try {
// double unitSize = unitSizeExp.evaluateConstant();
// if (unitSize != 1.0){
// throw new MappingException("model structure '"+sm.getStructure().getName()+"' unit size = "+unitSize+" != 1.0 ... partial volume or surface mapping not yet supported for particles");
// }
// }catch (ExpressionException e){
// e.printStackTrace(System.out);
// throw new MappingException("couldn't evaluate unit size for model structure '"+sm.getStructure().getName()+"' : "+e.getMessage());
// }
// }
// }
// }
// {
// GeometryClass[] geometryClass = getSimulationContext().getGeometryContext().getGeometry().getGeometryClasses();
// for (int i = 0; i < geometryClass.length; i++){
// Structure[] mappedStructures = getSimulationContext().getGeometryContext().getStructuresFromGeometryClass(geometryClass[i]);
// if (mappedStructures==null || mappedStructures.length==0){
// throw new MappingException("geometryClass '"+geometryClass[i].getName()+"' not mapped from a model structure");
// }
// }
// }
// deals with model parameters
Model model = getSimulationContext().getModel();
ModelUnitSystem modelUnitSystem = model.getUnitSystem();
ModelParameter[] modelParameters = model.getModelParameters();
// populate in globalParameterVariants hashtable
for (int j = 0; j < modelParameters.length; j++) {
Expression modelParamExpr = modelParameters[j].getExpression();
GeometryClass geometryClass = getDefaultGeometryClass(modelParamExpr);
modelParamExpr = getIdentifierSubstitutions(modelParamExpr, modelParameters[j].getUnitDefinition(), geometryClass);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], geometryClass), modelParamExpr, geometryClass));
}
//
// create new MathDescription (based on simContext's previous MathDescription if possible)
//
MathDescription oldMathDesc = getSimulationContext().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(getSimulationContext().getName() + "_generated");
}
//
// volume particle variables
//
Enumeration<SpeciesContextMapping> enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = enum1.nextElement();
if (scm.getVariable() instanceof ParticleVariable) {
if (!(mathDesc.getVariable(scm.getVariable().getName()) instanceof ParticleVariable)) {
varHash.addVariable(scm.getVariable());
}
}
}
varHash.addVariable(new Constant(getMathSymbol(model.getPI_CONSTANT(), null), getIdentifierSubstitutions(model.getPI_CONSTANT().getExpression(), model.getPI_CONSTANT().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(getSimulationContext().getTemperatureKelvin()), model.getTEMPERATURE().getUnitDefinition(), null)));
//
for (int j = 0; j < structureMappings.length; j++) {
if (structureMappings[j] instanceof MembraneMapping) {
MembraneMapping membraneMapping = (MembraneMapping) structureMappings[j];
GeometryClass geometryClass = membraneMapping.getGeometryClass();
//
// don't calculate voltage, still may need it though
//
Parameter initialVoltageParm = membraneMapping.getInitialVoltageParameter();
Variable voltageFunction = newFunctionOrConstant(getMathSymbol(membraneMapping.getMembrane().getMembraneVoltage(), geometryClass), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), geometryClass), geometryClass);
varHash.addVariable(voltageFunction);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(membraneMapping.getMembrane().getMembraneVoltage(), membraneMapping.getGeometryClass()), getIdentifierSubstitutions(membraneMapping.getInitialVoltageParameter().getExpression(), membraneMapping.getInitialVoltageParameter().getUnitDefinition(), membraneMapping.getGeometryClass()), membraneMapping.getGeometryClass()));
}
}
//
for (int j = 0; j < reactionSteps.length; j++) {
ReactionStep rs = reactionSteps[j];
if (getSimulationContext().getReactionContext().getReactionSpec(rs).isExcluded()) {
continue;
}
Kinetics.KineticsParameter[] parameters = rs.getKinetics().getKineticsParameters();
GeometryClass geometryClass = null;
if (rs.getStructure() != null) {
geometryClass = getSimulationContext().getGeometryContext().getStructureMapping(rs.getStructure()).getGeometryClass();
}
if (parameters != null) {
for (int i = 0; i < parameters.length; i++) {
// Reaction rate, currentDensity, LumpedCurrent and null parameters are not going to displayed in the particle math description.
if (((parameters[i].getRole() == Kinetics.ROLE_CurrentDensity) || (parameters[i].getRole() == Kinetics.ROLE_LumpedCurrent) || (parameters[i].getRole() == Kinetics.ROLE_ReactionRate)) || (parameters[i].getExpression() == null)) {
continue;
}
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[i], geometryClass), getIdentifierSubstitutions(parameters[i].getExpression(), parameters[i].getUnitDefinition(), geometryClass), geometryClass));
}
}
}
//
// initial constants (either function or constant)
//
SpeciesContextSpec[] speciesContextSpecs = getSimulationContext().getReactionContext().getSpeciesContextSpecs();
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpecParameter initParm = null;
Expression initExpr = null;
if (getSimulationContext().isUsingConcentration()) {
initParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
initExpr = new Expression(initParm.getExpression());
// if (speciesContextSpecs[i].getSpeciesContext().getStructure() instanceof Feature) {
// initExpr = Expression.div(initExpr, new Expression(model.getKMOLE, getNameScope())).flatten();
// }
} else {
initParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialCount);
initExpr = new Expression(initParm.getExpression());
}
if (initExpr != null) {
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
String[] symbols = initExpr.getSymbols();
// Check if 'initExpr' has other speciesContexts in its expression, need to replace it with 'spContext_init'
for (int j = 0; symbols != null && j < symbols.length; j++) {
// if symbol is a speciesContext, replacing it with a reference to initial condition for that speciesContext.
SpeciesContext spC = null;
SymbolTableEntry ste = initExpr.getSymbolBinding(symbols[j]);
if (ste instanceof SpeciesContextSpecProxyParameter) {
SpeciesContextSpecProxyParameter spspp = (SpeciesContextSpecProxyParameter) ste;
if (spspp.getTarget() instanceof SpeciesContext) {
spC = (SpeciesContext) spspp.getTarget();
SpeciesContextSpec spcspec = getSimulationContext().getReactionContext().getSpeciesContextSpec(spC);
SpeciesContextSpecParameter spCInitParm = spcspec.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
// if initConc param expression is null, try initCount
if (spCInitParm.getExpression() == null) {
spCInitParm = spcspec.getParameterFromRole(SpeciesContextSpec.ROLE_InitialCount);
}
// need to get init condn expression, but can't get it from getMathSymbol() (mapping between bio and math), hence get it as below.
Expression scsInitExpr = new Expression(spCInitParm, getNameScope());
// scsInitExpr.bindExpression(this);
initExpr.substituteInPlace(new Expression(spC.getName()), scsInitExpr);
}
}
}
// now create the appropriate function for the current speciesContextSpec.
varHash.addVariable(newFunctionOrConstant(getMathSymbol(initParm, sm.getGeometryClass()), getIdentifierSubstitutions(initExpr, initParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter diffParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_DiffusionRate);
if (diffParm != null) {
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(diffParm, sm.getGeometryClass()), getIdentifierSubstitutions(diffParm.getExpression(), diffParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter bc_xm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueXm);
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
if (bc_xm != null && (bc_xm.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_xm, sm.getGeometryClass()), getIdentifierSubstitutions(bc_xm.getExpression(), bc_xm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_xp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueXp);
if (bc_xp != null && (bc_xp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_xp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_xp.getExpression(), bc_xp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_ym = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueYm);
if (bc_ym != null && (bc_ym.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_ym, sm.getGeometryClass()), getIdentifierSubstitutions(bc_ym.getExpression(), bc_ym.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_yp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueYp);
if (bc_yp != null && (bc_yp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_yp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_yp.getExpression(), bc_yp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_zm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueZm);
if (bc_zm != null && (bc_zm.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_zm, sm.getGeometryClass()), getIdentifierSubstitutions(bc_zm.getExpression(), bc_zm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_zp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueZp);
if (bc_zp != null && (bc_zp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_zp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_zp.getExpression(), bc_zp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter advection_velX = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityX);
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
GeometryClass geometryClass = sm.getGeometryClass();
if (advection_velX != null && (advection_velX.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velX, geometryClass), getIdentifierSubstitutions(advection_velX.getExpression(), advection_velX.getUnitDefinition(), geometryClass), geometryClass));
}
SpeciesContextSpec.SpeciesContextSpecParameter advection_velY = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityY);
if (advection_velY != null && (advection_velY.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velY, geometryClass), getIdentifierSubstitutions(advection_velY.getExpression(), advection_velY.getUnitDefinition(), geometryClass), geometryClass));
}
SpeciesContextSpec.SpeciesContextSpecParameter advection_velZ = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityZ);
if (advection_velZ != null && (advection_velZ.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velZ, geometryClass), getIdentifierSubstitutions(advection_velZ.getExpression(), advection_velZ.getUnitDefinition(), geometryClass), geometryClass));
}
}
//
// 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());
}
}
//
// conversion factors
//
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.getKMILLIVOLTS(), null), getIdentifierSubstitutions(model.getKMILLIVOLTS().getExpression(), model.getKMILLIVOLTS().getUnitDefinition(), null)));
varHash.addVariable(new Constant(getMathSymbol(model.getK_GHK(), null), getIdentifierSubstitutions(model.getK_GHK().getExpression(), model.getK_GHK().getUnitDefinition(), null)));
//
for (int i = 0; i < structureMappings.length; i++) {
StructureMapping sm = structureMappings[i];
if (getSimulationContext().getGeometry().getDimension() == 0) {
StructureMappingParameter sizeParm = sm.getSizeParameter();
if (sizeParm != null && sizeParm.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(sizeParm, sm.getGeometryClass()), getIdentifierSubstitutions(sizeParm.getExpression(), sizeParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
} else {
if (sm instanceof MembraneMapping) {
MembraneMapping mm = (MembraneMapping) sm;
StructureMappingParameter volFrac = mm.getVolumeFractionParameter();
if (volFrac != null && volFrac.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(volFrac, sm.getGeometryClass()), getIdentifierSubstitutions(volFrac.getExpression(), volFrac.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
StructureMappingParameter surfToVol = mm.getSurfaceToVolumeParameter();
if (surfToVol != null && surfToVol.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(surfToVol, sm.getGeometryClass()), getIdentifierSubstitutions(surfToVol.getExpression(), surfToVol.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
}
} else {
Parameter parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitArea);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitVolume);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitArea);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitVolume);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_Size);
if (parm != null && parm.getExpression() != null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
}
}
}
//
// functions
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() == null && scm.getDependencyExpression() != null) {
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(scm.getSpeciesContext().getStructure());
Variable dependentVariable = newFunctionOrConstant(getMathSymbol(scm.getSpeciesContext(), sm.getGeometryClass()), getIdentifierSubstitutions(scm.getDependencyExpression(), scm.getSpeciesContext().getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass());
dependentVariable.setDomain(new Domain(sm.getGeometryClass()));
varHash.addVariable(dependentVariable);
}
}
//
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
varHash.addVariable(newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass()));
}
}
//
// set Variables to MathDescription all at once with the order resolved by "VariableHash"
//
mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
//
if (getSimulationContext().getGeometryContext().getGeometry() != null) {
try {
mathDesc.setGeometry(getSimulationContext().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");
}
//
// create subdomains (volume and surfaces)
//
GeometryClass[] geometryClasses = getSimulationContext().getGeometryContext().getGeometry().getGeometryClasses();
for (int k = 0; k < geometryClasses.length; k++) {
if (geometryClasses[k] instanceof SubVolume) {
SubVolume subVolume = (SubVolume) geometryClasses[k];
//
// get priority of subDomain
//
// now does not have to match spatial feature, *BUT* needs to be unique
int priority = k;
//
// create subDomain
//
CompartmentSubDomain subDomain = new CompartmentSubDomain(subVolume.getName(), priority);
mathDesc.addSubDomain(subDomain);
//
// assign boundary condition types
//
StructureMapping[] mappedSMs = getSimulationContext().getGeometryContext().getStructureMappings(subVolume);
FeatureMapping mappedFM = null;
for (int i = 0; i < mappedSMs.length; i++) {
if (mappedSMs[i] instanceof FeatureMapping) {
if (mappedFM != null) {
lg.warn("WARNING:::: MathMapping.refreshMathDescription() ... assigning boundary condition types not unique");
}
mappedFM = (FeatureMapping) mappedSMs[i];
}
}
if (mappedFM != null) {
subDomain.setBoundaryConditionXm(mappedFM.getBoundaryConditionTypeXm());
subDomain.setBoundaryConditionXp(mappedFM.getBoundaryConditionTypeXp());
if (getSimulationContext().getGeometry().getDimension() > 1) {
subDomain.setBoundaryConditionYm(mappedFM.getBoundaryConditionTypeYm());
subDomain.setBoundaryConditionYp(mappedFM.getBoundaryConditionTypeYp());
}
if (getSimulationContext().getGeometry().getDimension() > 2) {
subDomain.setBoundaryConditionZm(mappedFM.getBoundaryConditionTypeZm());
subDomain.setBoundaryConditionZp(mappedFM.getBoundaryConditionTypeZp());
}
}
} else if (geometryClasses[k] instanceof SurfaceClass) {
SurfaceClass surfaceClass = (SurfaceClass) geometryClasses[k];
// determine membrane inside and outside subvolume
// this preserves backward compatibility so that membrane subdomain
// inside and outside correspond to structure hierarchy when present
Pair<SubVolume, SubVolume> ret = DiffEquMathMapping.computeBoundaryConditionSource(model, simContext, surfaceClass);
SubVolume innerSubVolume = ret.one;
SubVolume outerSubVolume = ret.two;
//
// create subDomain
//
CompartmentSubDomain outerCompartment = mathDesc.getCompartmentSubDomain(outerSubVolume.getName());
CompartmentSubDomain innerCompartment = mathDesc.getCompartmentSubDomain(innerSubVolume.getName());
MembraneSubDomain memSubDomain = new MembraneSubDomain(innerCompartment, outerCompartment, surfaceClass.getName());
mathDesc.addSubDomain(memSubDomain);
}
}
//
// create Particle Contexts for all Particle Variables
//
Enumeration<SpeciesContextMapping> enumSCM = getSpeciesContextMappings();
Expression unitFactor = getUnitFactor(modelUnitSystem.getStochasticSubstanceUnit().divideBy(modelUnitSystem.getVolumeSubstanceUnit()));
while (enumSCM.hasMoreElements()) {
SpeciesContextMapping scm = enumSCM.nextElement();
SpeciesContext sc = scm.getSpeciesContext();
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure());
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
if (scm.getVariable() instanceof ParticleVariable && scm.getDependencyExpression() == null) {
ParticleVariable particleVariable = (ParticleVariable) scm.getVariable();
//
// initial distribution of particles
//
ArrayList<ParticleInitialCondition> particleInitialConditions = new ArrayList<ParticleInitialCondition>();
ParticleInitialCondition pic = null;
if (getSimulationContext().isUsingConcentration()) {
Expression initialDistribution = scs.getInitialConcentrationParameter().getExpression() == null ? null : new Expression(getMathSymbol(scs.getInitialConcentrationParameter(), sm.getGeometryClass()));
if (particleVariable instanceof VolumeParticleVariable) {
initialDistribution = Expression.mult(initialDistribution, unitFactor);
}
pic = new ParticleInitialConditionConcentration(initialDistribution);
} else {
Expression initialCount = scs.getInitialCountParameter().getExpression() == null ? null : new Expression(getMathSymbol(scs.getInitialCountParameter(), sm.getGeometryClass()));
if (initialCount == null) {
throw new MappingException("initialCount not defined for speciesContext " + scs.getSpeciesContext().getName());
}
Expression locationX = new Expression("u");
Expression locationY = new Expression("u");
Expression locationZ = new Expression("u");
pic = new ParticleInitialConditionCount(initialCount, locationX, locationY, locationZ);
}
particleInitialConditions.add(pic);
//
// diffusion
//
Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), sm.getGeometryClass()));
Expression driftXExp = null;
if (scs.getVelocityXParameter().getExpression() != null) {
driftXExp = new Expression(getMathSymbol(scs.getVelocityXParameter(), sm.getGeometryClass()));
} else {
SpatialQuantity[] velX_quantities = scs.getVelocityQuantities(QuantityComponent.X);
if (velX_quantities.length > 0) {
int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(sm.getGeometryClass()).length;
if (velX_quantities.length == 1 && numRegions == 1) {
driftXExp = new Expression(getMathSymbol(velX_quantities[0], sm.getGeometryClass()));
} else {
throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
}
}
}
Expression driftYExp = null;
if (scs.getVelocityYParameter().getExpression() != null) {
driftYExp = new Expression(getMathSymbol(scs.getVelocityYParameter(), sm.getGeometryClass()));
} else {
SpatialQuantity[] velY_quantities = scs.getVelocityQuantities(QuantityComponent.Y);
if (velY_quantities.length > 0) {
int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(sm.getGeometryClass()).length;
if (velY_quantities.length == 1 && numRegions == 1) {
driftYExp = new Expression(getMathSymbol(velY_quantities[0], sm.getGeometryClass()));
} else {
throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
}
}
}
Expression driftZExp = null;
if (scs.getVelocityZParameter().getExpression() != null) {
driftZExp = new Expression(getMathSymbol(scs.getVelocityZParameter(), sm.getGeometryClass()));
} else {
SpatialQuantity[] velZ_quantities = scs.getVelocityQuantities(QuantityComponent.Z);
if (velZ_quantities.length > 0) {
int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(sm.getGeometryClass()).length;
if (velZ_quantities.length == 1 && numRegions == 1) {
driftZExp = new Expression(getMathSymbol(velZ_quantities[0], sm.getGeometryClass()));
} else {
throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
}
}
}
ParticleProperties particleProperties = new ParticleProperties(particleVariable, diffusion, driftXExp, driftYExp, driftZExp, particleInitialConditions);
GeometryClass myGC = sm.getGeometryClass();
if (myGC == null) {
throw new MappingException("Application '" + getSimulationContext().getName() + "'\nGeometry->StructureMapping->(" + sm.getStructure().getTypeName() + ")'" + sm.getStructure().getName() + "' must be mapped to geometry domain.\n(see 'Problems' tab)");
}
SubDomain subDomain = mathDesc.getSubDomain(myGC.getName());
subDomain.addParticleProperties(particleProperties);
}
}
for (ReactionStep reactionStep : reactionSteps) {
Kinetics kinetics = reactionStep.getKinetics();
StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(reactionStep.getStructure());
GeometryClass reactionStepGeometryClass = sm.getGeometryClass();
SubDomain subdomain = mathDesc.getSubDomain(reactionStepGeometryClass.getName());
KineticsParameter reactionRateParameter = null;
if (kinetics instanceof LumpedKinetics) {
reactionRateParameter = ((LumpedKinetics) kinetics).getLumpedReactionRateParameter();
} else {
reactionRateParameter = ((DistributedKinetics) kinetics).getReactionRateParameter();
}
// macroscopic_irreversible/Microscopic_irreversible for bimolecular membrane reactions. They will NOT go through MassAction solver.
if (kinetics.getKineticsDescription().equals(KineticsDescription.Macroscopic_irreversible) || kinetics.getKineticsDescription().equals(KineticsDescription.Microscopic_irreversible)) {
Expression radiusExp = getIdentifierSubstitutions(reactionStep.getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_Binding_Radius).getExpression(), modelUnitSystem.getBindingRadiusUnit(), reactionStepGeometryClass);
if (radiusExp != null) {
Expression expCopy = new Expression(radiusExp);
try {
MassActionSolver.substituteParameters(expCopy, true).evaluateConstant();
} catch (ExpressionException e) {
throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Problem in binding radius of " + reactionStep.getName() + ": '" + radiusExp.infix() + "', " + e.getMessage()));
}
} else {
throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Binding radius of " + reactionStep.getName() + " is null."));
}
List<ParticleVariable> reactantParticles = new ArrayList<ParticleVariable>();
List<ParticleVariable> productParticles = new ArrayList<ParticleVariable>();
List<Action> forwardActions = new ArrayList<Action>();
for (ReactionParticipant rp : reactionStep.getReactionParticipants()) {
SpeciesContext sc = rp.getSpeciesContext();
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
GeometryClass scGeometryClass = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass();
String varName = getMathSymbol(sc, scGeometryClass);
Variable var = mathDesc.getVariable(varName);
if (var instanceof ParticleVariable) {
ParticleVariable particle = (ParticleVariable) var;
if (rp instanceof Reactant) {
reactantParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (radiusExp != null) {
forwardActions.add(Action.createDestroyAction(particle));
}
}
}
} else if (rp instanceof Product) {
productParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (radiusExp != null) {
forwardActions.add(Action.createCreateAction(particle));
}
}
}
}
} else {
throw new MappingException("particle variable '" + varName + "' not found");
}
}
JumpProcessRateDefinition bindingRadius = new InteractionRadius(radiusExp);
// get jump process name
String jpName = TokenMangler.mangleToSName(reactionStep.getName());
// only for NFSim/Rules for now.
ProcessSymmetryFactor processSymmetryFactor = null;
if (forwardActions.size() > 0) {
ParticleJumpProcess forwardProcess = new ParticleJumpProcess(jpName, reactantParticles, bindingRadius, forwardActions, processSymmetryFactor);
subdomain.addParticleJumpProcess(forwardProcess);
}
} else // other type of reactions
{
/* 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;
// Using the MassActionFunction to write out the math description
MassActionSolver.MassActionFunction maFunc = null;
if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction) || kinetics.getKineticsDescription().equals(KineticsDescription.General) || kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
Expression rateExp = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate).getExpression();
Parameter forwardRateParameter = null;
Parameter reverseRateParameter = null;
if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KForward);
reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KReverse);
} else if (kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
}
maFunc = MassActionSolver.solveMassAction(forwardRateParameter, reverseRateParameter, 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();
}
}
}
if (maFunc != null) {
// if the reaction has forward rate (Mass action,HMMs), or don't have either forward or reverse rate (some other rate laws--like general)
// we process it as forward reaction
List<ParticleVariable> reactantParticles = new ArrayList<ParticleVariable>();
List<ParticleVariable> productParticles = new ArrayList<ParticleVariable>();
List<Action> forwardActions = new ArrayList<Action>();
List<Action> reverseActions = new ArrayList<Action>();
List<ReactionParticipant> reactants = maFunc.getReactants();
List<ReactionParticipant> products = maFunc.getProducts();
for (ReactionParticipant rp : reactants) {
SpeciesContext sc = rp.getSpeciesContext();
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
GeometryClass scGeometryClass = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass();
String varName = getMathSymbol(sc, scGeometryClass);
Variable var = mathDesc.getVariable(varName);
if (var instanceof ParticleVariable) {
ParticleVariable particle = (ParticleVariable) var;
reactantParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (forwardRate != null) {
forwardActions.add(Action.createDestroyAction(particle));
}
if (reverseRate != null) {
reverseActions.add(Action.createCreateAction(particle));
}
}
}
} else {
throw new MappingException("particle variable '" + varName + "' not found");
}
}
for (ReactionParticipant rp : products) {
SpeciesContext sc = rp.getSpeciesContext();
SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
GeometryClass scGeometryClass = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass();
String varName = getMathSymbol(sc, scGeometryClass);
Variable var = mathDesc.getVariable(varName);
if (var instanceof ParticleVariable) {
ParticleVariable particle = (ParticleVariable) var;
productParticles.add(particle);
if (!scs.isConstant() && !scs.isForceContinuous()) {
for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
if (forwardRate != null) {
forwardActions.add(Action.createCreateAction(particle));
}
if (reverseRate != null) {
reverseActions.add(Action.createDestroyAction(particle));
}
}
}
} else {
throw new MappingException("particle variable '" + varName + "' not found");
}
}
//
// There are two unit conversions required:
//
// 1) convert entire reaction rate from vcell reaction units to Smoldyn units (molecules/lengthunit^dim/timeunit)
// (where dim is 2 for membrane reactions and 3 for volume reactions)
//
// for forward rates:
// 2) convert each reactant from Smoldyn units (molecules/lengthunit^dim) to VCell units
// (where dim is 2 for membrane reactants and 3 for volume reactants)
//
// or
//
// for reverse rates:
// 2) convert each product from Smoldyn units (molecules/lengthunit^dim) to VCell units
// (where dim is 2 for membrane products and 3 for volume products)
//
RationalNumber reactionLocationDim = new RationalNumber(reactionStep.getStructure().getDimension());
VCUnitDefinition timeUnit = modelUnitSystem.getTimeUnit();
VCUnitDefinition smoldynReactionSizeUnit = modelUnitSystem.getLengthUnit().raiseTo(reactionLocationDim);
VCUnitDefinition smoldynSubstanceUnit = modelUnitSystem.getStochasticSubstanceUnit();
VCUnitDefinition smoldynReactionRateUnit = smoldynSubstanceUnit.divideBy(smoldynReactionSizeUnit).divideBy(timeUnit);
VCUnitDefinition vcellReactionRateUnit = reactionRateParameter.getUnitDefinition();
VCUnitDefinition reactionUnitFactor = smoldynReactionRateUnit.divideBy(vcellReactionRateUnit);
if (forwardRate != null) {
VCUnitDefinition smoldynReactantsUnit = modelUnitSystem.getInstance_DIMENSIONLESS();
// start with factor to translate entire reaction rate.
VCUnitDefinition forwardUnitFactor = reactionUnitFactor;
//
for (ReactionParticipant reactant : maFunc.getReactants()) {
VCUnitDefinition vcellReactantUnit = reactant.getSpeciesContext().getUnitDefinition();
boolean bForceContinuous = simContext.getReactionContext().getSpeciesContextSpec(reactant.getSpeciesContext()).isForceContinuous();
VCUnitDefinition smoldynReactantUnit = null;
if (bForceContinuous) {
// reactant is continuous (vcell units)
smoldynReactantUnit = reactant.getSpeciesContext().getUnitDefinition();
} else {
// reactant is a particle (smoldyn units)
RationalNumber reactantLocationDim = new RationalNumber(reactant.getStructure().getDimension());
VCUnitDefinition smoldynReactantSize = modelUnitSystem.getLengthUnit().raiseTo(reactantLocationDim);
smoldynReactantUnit = smoldynSubstanceUnit.divideBy(smoldynReactantSize);
}
// keep track of units of all reactants
smoldynReactantsUnit = smoldynReactantsUnit.multiplyBy(smoldynReactantUnit);
RationalNumber reactantStoichiometry = new RationalNumber(reactant.getStoichiometry());
VCUnitDefinition reactantUnitFactor = (vcellReactantUnit.divideBy(smoldynReactantUnit)).raiseTo(reactantStoichiometry);
// accumulate unit factors for all reactants
forwardUnitFactor = forwardUnitFactor.multiplyBy(reactantUnitFactor);
}
forwardRate = Expression.mult(forwardRate, getUnitFactor(forwardUnitFactor));
VCUnitDefinition smoldynExpectedForwardRateUnit = smoldynReactionRateUnit.divideBy(smoldynReactantsUnit);
// get probability
Expression exp = getIdentifierSubstitutions(forwardRate, smoldynExpectedForwardRateUnit, reactionStepGeometryClass).flatten();
JumpProcessRateDefinition partRateDef = new MacroscopicRateConstant(exp);
// create particle jump process
String jpName = TokenMangler.mangleToSName(reactionStep.getName());
// only for NFSim/Rules for now.
ProcessSymmetryFactor processSymmetryFactor = null;
if (forwardActions.size() > 0) {
ParticleJumpProcess forwardProcess = new ParticleJumpProcess(jpName, reactantParticles, partRateDef, forwardActions, processSymmetryFactor);
subdomain.addParticleJumpProcess(forwardProcess);
}
}
// end of forward rate not null
if (reverseRate != null) {
VCUnitDefinition smoldynProductsUnit = modelUnitSystem.getInstance_DIMENSIONLESS();
// start with factor to translate entire reaction rate.
VCUnitDefinition reverseUnitFactor = reactionUnitFactor;
//
for (ReactionParticipant product : maFunc.getProducts()) {
VCUnitDefinition vcellProductUnit = product.getSpeciesContext().getUnitDefinition();
boolean bForceContinuous = simContext.getReactionContext().getSpeciesContextSpec(product.getSpeciesContext()).isForceContinuous();
VCUnitDefinition smoldynProductUnit = null;
if (bForceContinuous) {
smoldynProductUnit = product.getSpeciesContext().getUnitDefinition();
} else {
RationalNumber productLocationDim = new RationalNumber(product.getStructure().getDimension());
VCUnitDefinition smoldynProductSize = modelUnitSystem.getLengthUnit().raiseTo(productLocationDim);
smoldynProductUnit = smoldynSubstanceUnit.divideBy(smoldynProductSize);
}
// keep track of units of all products
smoldynProductsUnit = smoldynProductsUnit.multiplyBy(smoldynProductUnit);
RationalNumber productStoichiometry = new RationalNumber(product.getStoichiometry());
VCUnitDefinition productUnitFactor = (vcellProductUnit.divideBy(smoldynProductUnit)).raiseTo(productStoichiometry);
// accumulate unit factors for all products
reverseUnitFactor = reverseUnitFactor.multiplyBy(productUnitFactor);
}
reverseRate = Expression.mult(reverseRate, getUnitFactor(reverseUnitFactor));
VCUnitDefinition smoldynExpectedReverseRateUnit = smoldynReactionRateUnit.divideBy(smoldynProductsUnit);
// get probability
Expression exp = getIdentifierSubstitutions(reverseRate, smoldynExpectedReverseRateUnit, reactionStepGeometryClass).flatten();
JumpProcessRateDefinition partProbRate = new MacroscopicRateConstant(exp);
// get jump process name
String jpName = TokenMangler.mangleToSName(reactionStep.getName() + "_reverse");
// only for NFSim/Rules for now.
ProcessSymmetryFactor processSymmetryFactor = null;
if (reverseActions.size() > 0) {
ParticleJumpProcess reverseProcess = new ParticleJumpProcess(jpName, productParticles, partProbRate, reverseActions, processSymmetryFactor);
subdomain.addParticleJumpProcess(reverseProcess);
}
}
// end of reverse rate not null
}
// end of maFunc not null
}
// end of reaction step for loop
}
//
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
Variable variable = newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass());
if (mathDesc.getVariable(variable.getName()) == null) {
mathDesc.addVariable(variable);
}
}
}
if (!mathDesc.isValid()) {
lg.warn(mathDesc.getVCML_database());
throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
}
if (lg.isDebugEnabled()) {
System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string begin ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
System.out.println(mathDesc.getVCML());
System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string end ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
}
}
Also used :
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
LumpedKinetics(cbit.vcell.model.LumpedKinetics)
MathDescription(cbit.vcell.math.MathDescription)
ArrayList(java.util.ArrayList)
Product(cbit.vcell.model.Product)
SpeciesContext(cbit.vcell.model.SpeciesContext)
StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter)
Reactant(cbit.vcell.model.Reactant)
ExpressionException(cbit.vcell.parser.ExpressionException)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
SubVolume(cbit.vcell.geometry.SubVolume)
Vector(java.util.Vector)
SpatialQuantity(cbit.vcell.mapping.spatial.SpatialObject.SpatialQuantity)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
JumpProcessRateDefinition(cbit.vcell.math.JumpProcessRateDefinition)
InteractionRadius(cbit.vcell.math.InteractionRadius)
ParticleJumpProcess(cbit.vcell.math.ParticleJumpProcess)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
ParticleInitialCondition(cbit.vcell.math.ParticleProperties.ParticleInitialCondition)
ReactionStep(cbit.vcell.model.ReactionStep)
ParticleProperties(cbit.vcell.math.ParticleProperties)
Kinetics(cbit.vcell.model.Kinetics)
DistributedKinetics(cbit.vcell.model.DistributedKinetics)
LumpedKinetics(cbit.vcell.model.LumpedKinetics)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
Domain(cbit.vcell.math.Variable.Domain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
ReactionParticipant(cbit.vcell.model.ReactionParticipant)
GeometryClass(cbit.vcell.geometry.GeometryClass)
Action(cbit.vcell.math.Action)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable)
ParticleVariable(cbit.vcell.math.ParticleVariable)
Variable(cbit.vcell.math.Variable)
SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter)
SurfaceClass(cbit.vcell.geometry.SurfaceClass)
VariableHash(cbit.vcell.math.VariableHash)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable)
ParticleVariable(cbit.vcell.math.ParticleVariable)
MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant)
Constant(cbit.vcell.math.Constant)
SymbolTableEntry(cbit.vcell.parser.SymbolTableEntry)
MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant)
RationalNumber(ucar.units_vcell.RationalNumber)
ModelUnitSystem(cbit.vcell.model.ModelUnitSystem)
Pair(org.vcell.util.Pair)
ParticleInitialConditionConcentration(cbit.vcell.math.ParticleProperties.ParticleInitialConditionConcentration)
ProcessSymmetryFactor(cbit.vcell.math.ParticleJumpProcess.ProcessSymmetryFactor)
Expression(cbit.vcell.parser.Expression)
VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable)
MathException(cbit.vcell.math.MathException)
Model(cbit.vcell.model.Model)
StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter)
Parameter(cbit.vcell.model.Parameter)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
MassActionSolver(cbit.vcell.model.MassActionSolver)
ParticleInitialConditionCount(cbit.vcell.math.ParticleProperties.ParticleInitialConditionCount)
Example 4 with VariableHash
use of cbit.vcell.math.VariableHash in project vcell by virtualcell.
the class XmlReader method readParameters.
private void readParameters(List<Element> parameterElements, ParameterContext parameterContext, HashMap<String, ParameterRoleEnum> roleHash, ParameterRoleEnum userDefinedRole, HashSet<String> xmlRolesTagsToIgnore, Model model) throws XmlParseException {
String contextName = parameterContext.getNameScope().getName();
try {
//
// prepopulate varHash with reserved symbols
//
VariableHash varHash = new VariableHash();
addResevedSymbols(varHash, model);
//
for (Element xmlParam : parameterElements) {
String parsedParamName = unMangle(xmlParam.getAttributeValue(XMLTags.NameAttrTag));
String parsedRoleString = xmlParam.getAttributeValue(XMLTags.ParamRoleAttrTag);
String parsedExpressionString = xmlParam.getText();
//
if (xmlRolesTagsToIgnore.contains(parsedRoleString)) {
varHash.removeVariable(parsedParamName);
continue;
}
Expression paramExp = null;
if (parsedExpressionString.trim().length() > 0) {
paramExp = unMangleExpression(parsedExpressionString);
}
if (varHash.getVariable(parsedParamName) == null) {
Domain domain = null;
varHash.addVariable(new Function(parsedParamName, paramExp, domain));
} else {
if (model.getReservedSymbolByName(parsedParamName) != null) {
varHash.removeVariable(parsedParamName);
Domain domain = null;
varHash.addVariable(new Function(parsedParamName, paramExp, domain));
}
}
//
// get the parameter for this xml role string
//
ParameterRoleEnum paramRole = roleHash.get(parsedRoleString);
if (paramRole == null) {
throw new XmlParseException("parameter '" + parsedParamName + "' has unexpected role '" + parsedRoleString + "' in '" + contextName + "'");
}
//
if (paramRole != userDefinedRole) {
LocalParameter paramWithSameRole = parameterContext.getLocalParameterFromRole(paramRole);
if (paramWithSameRole == null) {
throw new XmlParseException("can't find parameter with role '" + parsedRoleString + "' in '" + contextName + "'");
}
//
if (!paramWithSameRole.getName().equals(parsedParamName)) {
//
// first rename other parameters with same name
//
LocalParameter paramWithSameNameButDifferentRole = parameterContext.getLocalParameterFromName(parsedParamName);
if (paramWithSameNameButDifferentRole != null) {
//
// find available name
//
int n = 0;
String newName = parsedParamName + "_" + n++;
while (parameterContext.getEntry(newName) != null) {
newName = parsedParamName + "_" + n++;
}
parameterContext.renameLocalParameter(parsedParamName, newName);
}
//
// then rename parameter with correct role
//
parameterContext.renameLocalParameter(paramWithSameRole.getName(), parsedParamName);
}
}
}
//
// create unresolved parameters for all unresolved symbols
//
String unresolvedSymbol = varHash.getFirstUnresolvedSymbol();
while (unresolvedSymbol != null) {
try {
Domain domain = null;
// will turn into an UnresolvedParameter.
varHash.addVariable(new Function(unresolvedSymbol, new Expression(0.0), domain));
} catch (MathException e) {
e.printStackTrace(System.out);
throw new XmlParseException(e.getMessage());
}
parameterContext.addUnresolvedParameter(unresolvedSymbol);
unresolvedSymbol = varHash.getFirstUnresolvedSymbol();
}
//
// in topological order, add parameters to model (getting units also).
// note that all pre-defined parameters already have the correct names
// here we set expressions on pre-defined parameters and add user-defined parameters
//
Variable[] sortedVariables = varHash.getTopologicallyReorderedVariables();
ModelUnitSystem modelUnitSystem = model.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 < parameterElements.size(); j++) {
Element tempParam = (Element) parameterElements.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);
}
LocalParameter tempParam = parameterContext.getLocalParameterFromName(paramFunction.getName());
if (tempParam == null) {
tempParam = parameterContext.addLocalParameter(paramFunction.getName(), new Expression(0.0), userDefinedRole, unit, userDefinedRole.getDescription());
parameterContext.setParameterValue(tempParam, paramFunction.getExpression(), true);
} else {
if (tempParam.getExpression() != null) {
// if the expression is null, it should remain null.
parameterContext.setParameterValue(tempParam, paramFunction.getExpression(), true);
}
tempParam.setUnitDefinition(unit);
}
}
}
} catch (PropertyVetoException | ExpressionException | MathException e) {
e.printStackTrace(System.out);
throw new XmlParseException("Exception while setting parameters for '" + contextName + "': " + e.getMessage(), e);
}
}
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)
VariableHash(cbit.vcell.math.VariableHash)
Element(org.jdom.Element)
ExpressionException(cbit.vcell.parser.ExpressionException)
LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter)
PropertyVetoException(java.beans.PropertyVetoException)
AnnotatedFunction(cbit.vcell.solver.AnnotatedFunction)
Function(cbit.vcell.math.Function)
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
Expression(cbit.vcell.parser.Expression)
MathException(cbit.vcell.math.MathException)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
FilamentSubDomain(cbit.vcell.math.FilamentSubDomain)
PointSubDomain(cbit.vcell.math.PointSubDomain)
Domain(cbit.vcell.math.Variable.Domain)
ParameterRoleEnum(cbit.vcell.mapping.ParameterContext.ParameterRoleEnum)
ModelUnitSystem(cbit.vcell.model.ModelUnitSystem)
Example 5 with VariableHash
use of cbit.vcell.math.VariableHash in project vcell by virtualcell.
the class MathMapping_4_8 method refreshMathDescription.
/**
* This method was created in VisualAge.
*/
private void refreshMathDescription() throws MappingException, MatrixException, MathException, ExpressionException, ModelException {
// All sizes must be set for new ODE models and ratios must be set for old ones.
simContext.checkValidity();
//
// temporarily place all variables in a hashtable (before binding) and discarding duplicates (check for equality)
//
VariableHash varHash = new VariableHash();
StructureMapping[] structureMappings = simContext.getGeometryContext().getStructureMappings();
Model model = simContext.getModel();
StructureTopology structTopology = model.getStructureTopology();
//
// 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 subdomain");
}
if (sm != null && (sm instanceof MembraneMapping) && ((MembraneMapping) sm).getVolumeFractionParameter() != null) {
Expression volFractExp = ((MembraneMapping) sm).getVolumeFractionParameter().getExpression();
if (volFractExp != null) {
try {
double volFract = volFractExp.evaluateConstant();
if (volFract >= 1.0) {
throw new MappingException("model structure '" + structTopology.getInsideFeature(((MembraneMapping) sm).getMembrane()).getName() + "' has volume fraction >= 1.0");
}
} catch (ExpressionException e) {
}
}
}
}
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 subdomain '" + subVolumes[i].getName() + "' not mapped from a model structure");
}
}
// deals with model parameters
Hashtable<VolVariable, EventAssignmentInitParameter> eventVolVarHash = new Hashtable<VolVariable, EventAssignmentInitParameter>();
ModelParameter[] modelParameters = model.getModelParameters();
if (simContext.getGeometry().getDimension() == 0) {
//
// global parameters from model (that presently are constants)
//
BioEvent[] bioEvents = simContext.getBioEvents();
ArrayList<SymbolTableEntry> eventAssignTargets = new ArrayList<SymbolTableEntry>();
if (bioEvents != null && bioEvents.length > 0) {
for (BioEvent be : bioEvents) {
for (EventAssignment ea : be.getEventAssignments()) {
if (!eventAssignTargets.contains(ea.getTarget())) {
eventAssignTargets.add(ea.getTarget());
}
}
}
}
for (int j = 0; j < modelParameters.length; j++) {
Expression modelParamExpr = getIdentifierSubstitutions(modelParameters[j].getExpression(), modelParameters[j].getUnitDefinition(), null);
if (eventAssignTargets.contains(modelParameters[j])) {
EventAssignmentInitParameter eap = null;
try {
eap = addEventAssignmentInitParameter(modelParameters[j].getName(), modelParameters[j].getExpression(), PARAMETER_ROLE_EVENTASSIGN_INITCONDN, modelParameters[j].getUnitDefinition());
} catch (PropertyVetoException e) {
e.printStackTrace(System.out);
throw new MappingException(e.getMessage());
}
// varHash.addVariable(newFunctionOrConstant(getMathSymbol(eap, null), modelParamExpr));
VolVariable volVar = new VolVariable(modelParameters[j].getName(), nullDomain);
varHash.addVariable(volVar);
eventVolVarHash.put(volVar, eap);
} else {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], null), modelParamExpr));
}
}
} else {
//
for (int pass = 0; pass < 2; pass++) {
for (int j = 0; j < modelParameters.length; j++) {
Hashtable<String, Expression> structMappingVariantsHash = new Hashtable<String, Expression>();
for (int k = 0; k < structureMappings.length; k++) {
String paramVariantName = null;
Expression paramVariantExpr = null;
if (modelParameters[j].getExpression().getSymbols() == null) {
paramVariantName = modelParameters[j].getName();
paramVariantExpr = getIdentifierSubstitutions(modelParameters[j].getExpression(), modelParameters[j].getUnitDefinition(), null);
} else {
paramVariantName = modelParameters[j].getName() + "_" + TokenMangler.fixTokenStrict(structureMappings[k].getStructure().getName());
// if the expression has symbols that do not belong in that structureMapping, do not create the variant.
Expression exp1 = modelParameters[j].getExpression();
Expression flattenedModelParamExpr = substituteGlobalParameters(exp1);
String[] symbols = flattenedModelParamExpr.getSymbols();
boolean bValid = true;
Structure sm_struct = structureMappings[k].getStructure();
if (symbols != null) {
for (int ii = 0; ii < symbols.length; ii++) {
SpeciesContext sc = model.getSpeciesContext(symbols[ii]);
if (sc != null) {
// symbol[ii] is a speciesContext, check its structure with structureMapping[k].structure. If they are the same or
// if it is the adjacent membrane(s), allow variant expression to be created. Else, continue.
Structure sp_struct = sc.getStructure();
if (sp_struct.compareEqual(sm_struct)) {
bValid = bValid && true;
} else {
// if the 2 structures are not the same, are they adjacent? then 'bValid' is true, else false.
if ((sm_struct instanceof Feature) && (sp_struct instanceof Membrane)) {
Feature sm_feature = (Feature) sm_struct;
Membrane sp_mem = (Membrane) sp_struct;
if (sp_mem.compareEqual(structTopology.getParentStructure(sm_feature)) || (structTopology.getInsideFeature(sp_mem).compareEqual(sm_feature) || structTopology.getOutsideFeature(sp_mem).compareEqual(sm_feature))) {
bValid = bValid && true;
} else {
bValid = bValid && false;
break;
}
} else if ((sm_struct instanceof Membrane) && (sp_struct instanceof Feature)) {
Feature sp_feature = (Feature) sp_struct;
Membrane sm_mem = (Membrane) sm_struct;
if (sm_mem.compareEqual(structTopology.getParentStructure(sp_feature)) || (structTopology.getInsideFeature(sm_mem).compareEqual(sp_feature) || structTopology.getOutsideFeature(sm_mem).compareEqual(sp_feature))) {
bValid = bValid && true;
} else {
bValid = bValid && false;
break;
}
} else {
bValid = bValid && false;
break;
}
}
}
}
}
if (bValid) {
if (pass == 0) {
paramVariantExpr = new Expression("VCELL_TEMPORARY_EXPRESSION_PLACEHOLDER");
} else {
paramVariantExpr = getIdentifierSubstitutions(modelParameters[j].getExpression(), modelParameters[j].getUnitDefinition(), structureMappings[k]);
}
}
}
if (paramVariantExpr != null) {
structMappingVariantsHash.put(paramVariantName, paramVariantExpr);
}
}
globalParamVariantsHash.put(modelParameters[j], structMappingVariantsHash);
}
}
//
for (int j = 0; j < modelParameters.length; j++) {
if (modelParameters[j].getExpression().getSymbols() == null) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], null), getIdentifierSubstitutions(modelParameters[j].getExpression(), modelParameters[j].getUnitDefinition(), null)));
} else {
Hashtable<String, Expression> smVariantsHash = globalParamVariantsHash.get(modelParameters[j]);
for (int k = 0; k < structureMappings.length; k++) {
String variantName = modelParameters[j].getName() + "_" + TokenMangler.fixTokenStrict(structureMappings[k].getStructure().getName());
Expression variantExpr = smVariantsHash.get(variantName);
if (variantExpr != null) {
varHash.addVariable(newFunctionOrConstant(variantName, variantExpr));
}
}
}
}
}
//
// 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");
}
//
// volume variables
//
Enumeration<SpeciesContextMapping> enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = enum1.nextElement();
if (scm.getVariable() instanceof VolVariable) {
if (!(mathDesc.getVariable(scm.getVariable().getName()) instanceof VolVariable)) {
varHash.addVariable(scm.getVariable());
}
}
}
//
// membrane variables
//
enum1 = getSpeciesContextMappings();
while (enum1.hasMoreElements()) {
SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
if (scm.getVariable() instanceof MemVariable) {
varHash.addVariable(scm.getVariable());
}
}
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)));
//
// only calculate potential if at least one MembraneMapping has CalculateVoltage == true
//
boolean bCalculatePotential = false;
for (int i = 0; i < structureMappings.length; i++) {
if (structureMappings[i] instanceof MembraneMapping) {
if (((MembraneMapping) structureMappings[i]).getCalculateVoltage()) {
bCalculatePotential = true;
}
}
}
// (simContext.getGeometry().getDimension() == 0);
potentialMapping = new PotentialMapping(simContext, this);
potentialMapping.computeMath();
if (bCalculatePotential) {
//
// copy functions for currents and constants for capacitances
//
ElectricalDevice[] devices = potentialMapping.getElectricalDevices();
for (int j = 0; j < devices.length; j++) {
if (devices[j] instanceof MembraneElectricalDevice) {
MembraneElectricalDevice membraneElectricalDevice = (MembraneElectricalDevice) devices[j];
MembraneMapping memMapping = membraneElectricalDevice.getMembraneMapping();
Parameter specificCapacitanceParm = memMapping.getParameterFromRole(MembraneMapping.ROLE_SpecificCapacitance);
varHash.addVariable(new Constant(getMathSymbol(specificCapacitanceParm, memMapping), getIdentifierSubstitutions(specificCapacitanceParm.getExpression(), specificCapacitanceParm.getUnitDefinition(), memMapping)));
ElectricalDevice.ElectricalDeviceParameter transmembraneCurrentParm = membraneElectricalDevice.getParameterFromRole(ElectricalDevice.ROLE_TransmembraneCurrent);
ElectricalDevice.ElectricalDeviceParameter totalCurrentParm = membraneElectricalDevice.getParameterFromRole(ElectricalDevice.ROLE_TotalCurrent);
ElectricalDevice.ElectricalDeviceParameter capacitanceParm = membraneElectricalDevice.getParameterFromRole(ElectricalDevice.ROLE_Capacitance);
if (totalCurrentParm != null && /* totalCurrentDensityParm.getExpression()!=null && */
memMapping.getCalculateVoltage()) {
Expression totalCurrentDensityExp = (totalCurrentParm.getExpression() != null) ? (totalCurrentParm.getExpression()) : (new Expression(0.0));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrentParm, membraneElectricalDevice.getMembraneMapping()), getIdentifierSubstitutions(totalCurrentDensityExp, totalCurrentParm.getUnitDefinition(), membraneElectricalDevice.getMembraneMapping())));
}
if (transmembraneCurrentParm != null && transmembraneCurrentParm.getExpression() != null && memMapping.getCalculateVoltage()) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(transmembraneCurrentParm, membraneElectricalDevice.getMembraneMapping()), getIdentifierSubstitutions(transmembraneCurrentParm.getExpression(), transmembraneCurrentParm.getUnitDefinition(), membraneElectricalDevice.getMembraneMapping())));
}
if (capacitanceParm != null && capacitanceParm.getExpression() != null && memMapping.getCalculateVoltage()) {
StructureMappingParameter sizeParameter = membraneElectricalDevice.getMembraneMapping().getSizeParameter();
if (simContext.getGeometry().getDimension() == 0 && (sizeParameter.getExpression() == null || sizeParameter.getExpression().isZero())) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(capacitanceParm, membraneElectricalDevice.getMembraneMapping()), getIdentifierSubstitutions(Expression.mult(memMapping.getNullSizeParameterValue(), specificCapacitanceParm.getExpression()), capacitanceParm.getUnitDefinition(), membraneElectricalDevice.getMembraneMapping())));
} else {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(capacitanceParm, membraneElectricalDevice.getMembraneMapping()), getIdentifierSubstitutions(capacitanceParm.getExpression(), capacitanceParm.getUnitDefinition(), membraneElectricalDevice.getMembraneMapping())));
}
}
//
if (membraneElectricalDevice.getDependentVoltageExpression() == null) {
// is Voltage Independent?
StructureMapping.StructureMappingParameter initialVoltageParm = memMapping.getInitialVoltageParameter();
varHash.addVariable(newFunctionOrConstant(getMathSymbol(initialVoltageParm, memMapping), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), memMapping)));
} else //
// membrane forced potential
//
{
varHash.addVariable(newFunctionOrConstant(getMathSymbol(memMapping.getMembrane().getMembraneVoltage(), memMapping), getIdentifierSubstitutions(membraneElectricalDevice.getDependentVoltageExpression(), memMapping.getMembrane().getMembraneVoltage().getUnitDefinition(), memMapping)));
}
} else if (devices[j] instanceof CurrentClampElectricalDevice) {
CurrentClampElectricalDevice currentClampDevice = (CurrentClampElectricalDevice) devices[j];
// total current = current source (no capacitance)
Parameter totalCurrentParm = currentClampDevice.getParameterFromRole(CurrentClampElectricalDevice.ROLE_TotalCurrent);
Parameter currentParm = currentClampDevice.getParameterFromRole(CurrentClampElectricalDevice.ROLE_TransmembraneCurrent);
// Parameter dependentVoltage = currentClampDevice.getCurrentClampStimulus().getVoltageParameter();
varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrentParm, null), getIdentifierSubstitutions(totalCurrentParm.getExpression(), totalCurrentParm.getUnitDefinition(), null)));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(currentParm, null), getIdentifierSubstitutions(currentParm.getExpression(), currentParm.getUnitDefinition(), null)));
// varHash.addVariable(newFunctionOrConstant(getMathSymbol(dependentVoltage,null),getIdentifierSubstitutions(currentClampDevice.getDependentVoltageExpression(),dependentVoltage.getUnitDefinition(),null)));
//
// add user-defined parameters
//
ElectricalDevice.ElectricalDeviceParameter[] parameters = currentClampDevice.getParameters();
for (int k = 0; k < parameters.length; k++) {
if (parameters[k].getExpression() != null) {
// guards against voltage parameters that are "variable".
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[k], null), getIdentifierSubstitutions(parameters[k].getExpression(), parameters[k].getUnitDefinition(), null)));
}
}
} else if (devices[j] instanceof VoltageClampElectricalDevice) {
VoltageClampElectricalDevice voltageClampDevice = (VoltageClampElectricalDevice) devices[j];
// total current = current source (no capacitance)
Parameter totalCurrent = voltageClampDevice.getParameterFromRole(VoltageClampElectricalDevice.ROLE_TotalCurrent);
Parameter totalCurrentParm = voltageClampDevice.getParameterFromRole(VoltageClampElectricalDevice.ROLE_TotalCurrent);
Parameter voltageParm = voltageClampDevice.getParameterFromRole(VoltageClampElectricalDevice.ROLE_Voltage);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrent, null), getIdentifierSubstitutions(totalCurrent.getExpression(), totalCurrent.getUnitDefinition(), null)));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrentParm, null), getIdentifierSubstitutions(totalCurrentParm.getExpression(), totalCurrentParm.getUnitDefinition(), null)));
varHash.addVariable(newFunctionOrConstant(getMathSymbol(voltageParm, null), getIdentifierSubstitutions(voltageParm.getExpression(), voltageParm.getUnitDefinition(), null)));
//
// add user-defined parameters
//
ElectricalDevice.ElectricalDeviceParameter[] parameters = voltageClampDevice.getParameters();
for (int k = 0; k < parameters.length; k++) {
if (parameters[k].getRole() == ElectricalDevice.ROLE_UserDefined) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[k], null), getIdentifierSubstitutions(parameters[k].getExpression(), parameters[k].getUnitDefinition(), null)));
}
}
}
}
} else {
//
for (int j = 0; j < structureMappings.length; j++) {
if (structureMappings[j] instanceof MembraneMapping) {
MembraneMapping memMapping = (MembraneMapping) structureMappings[j];
varHash.addVariable(newFunctionOrConstant(getMathSymbol(memMapping.getMembrane().getMembraneVoltage(), memMapping), getIdentifierSubstitutions(memMapping.getInitialVoltageParameter().getExpression(), memMapping.getInitialVoltageParameter().getUnitDefinition(), memMapping)));
}
}
}
//
for (int j = 0; j < structureMappings.length; j++) {
if (structureMappings[j] instanceof MembraneMapping) {
MembraneMapping membraneMapping = (MembraneMapping) structureMappings[j];
Membrane.MembraneVoltage membraneVoltage = membraneMapping.getMembrane().getMembraneVoltage();
ElectricalDevice[] membraneDevices = potentialMapping.getElectricalDevices(membraneMapping.getMembrane());
// ElectricalDevice membraneDevice = null;
for (int i = 0; i < membraneDevices.length; i++) {
if (membraneDevices[i].hasCapacitance() && membraneDevices[i].getDependentVoltageExpression() == null) {
if (membraneMapping.getCalculateVoltage() && bCalculatePotential) {
if (getResolved(membraneMapping)) {
//
if (mathDesc.getVariable(Membrane.MEMBRANE_VOLTAGE_REGION_NAME) == null) {
// varHash.addVariable(new MembraneRegionVariable(MembraneVoltage.MEMBRANE_VOLTAGE_REGION_NAME));
varHash.addVariable(new MembraneRegionVariable(getMathSymbol(membraneVoltage, membraneMapping), nullDomain));
}
} else {
//
// spatially unresolved membrane, and must solve for potential ... make VolVariable for this compartment
//
varHash.addVariable(new VolVariable(getMathSymbol(membraneVoltage, membraneMapping), nullDomain));
}
Parameter initialVoltageParm = membraneMapping.getInitialVoltageParameter();
Variable initVoltageFunction = newFunctionOrConstant(getMathSymbol(initialVoltageParm, membraneMapping), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), membraneMapping));
varHash.addVariable(initVoltageFunction);
} else {
//
// don't calculate voltage, still may need it though
//
Parameter initialVoltageParm = membraneMapping.getInitialVoltageParameter();
Variable voltageFunction = newFunctionOrConstant(getMathSymbol(membraneMapping.getMembrane().getMembraneVoltage(), membraneMapping), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), membraneMapping));
varHash.addVariable(voltageFunction);
}
}
}
}
}
//
for (int j = 0; j < reactionSteps.length; j++) {
ReactionStep rs = reactionSteps[j];
if (simContext.getReactionContext().getReactionSpec(rs).isExcluded()) {
continue;
}
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].getRole() == Kinetics.ROLE_LumpedCurrent)) && (parameters[i].getExpression() == null || parameters[i].getExpression().isZero())) {
continue;
}
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[i], sm), getIdentifierSubstitutions(parameters[i].getExpression(), parameters[i].getUnitDefinition(), sm)));
}
}
}
//
// initial constants (either function or constant)
//
SpeciesContextSpec[] speciesContextSpecs = simContext.getReactionContext().getSpeciesContextSpecs();
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpecParameter initParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
if (initParm != null) {
Expression initExpr = new Expression(initParm.getExpression());
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
String[] symbols = initExpr.getSymbols();
// Check if 'initExpr' has other speciesContexts in its expression, need to replace it with 'spContext_init'
for (int j = 0; symbols != null && j < symbols.length; j++) {
// if symbol is a speciesContext, replacing it with a reference to initial condition for that speciesContext.
SpeciesContext spC = null;
SymbolTableEntry ste = initExpr.getSymbolBinding(symbols[j]);
if (ste instanceof SpeciesContextSpecProxyParameter) {
SpeciesContextSpecProxyParameter spspp = (SpeciesContextSpecProxyParameter) ste;
if (spspp.getTarget() instanceof SpeciesContext) {
spC = (SpeciesContext) spspp.getTarget();
SpeciesContextSpec spcspec = simContext.getReactionContext().getSpeciesContextSpec(spC);
SpeciesContextSpecParameter spCInitParm = spcspec.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
// if initConc param expression is null, try initCount
if (spCInitParm.getExpression() == null) {
spCInitParm = spcspec.getParameterFromRole(SpeciesContextSpec.ROLE_InitialCount);
}
// need to get init condn expression, but can't get it from getMathSymbol() (mapping between bio and math), hence get it as below.
Expression scsInitExpr = new Expression(spCInitParm, getNameScope());
// scsInitExpr.bindExpression(this);
initExpr.substituteInPlace(new Expression(spC.getName()), scsInitExpr);
}
}
}
// now create the appropriate function for the current speciesContextSpec.
varHash.addVariable(newFunctionOrConstant(getMathSymbol(initParm, sm), getIdentifierSubstitutions(initExpr, initParm.getUnitDefinition(), sm)));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextMapping scm = getSpeciesContextMapping(speciesContextSpecs[i].getSpeciesContext());
SpeciesContextSpec.SpeciesContextSpecParameter diffParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_DiffusionRate);
if (diffParm != null && (scm.isPDERequired())) {
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(diffParm, sm), getIdentifierSubstitutions(diffParm.getExpression(), diffParm.getUnitDefinition(), sm)));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter bc_xm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueXm);
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
if (bc_xm != null && (bc_xm.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_xm, sm), getIdentifierSubstitutions(bc_xm.getExpression(), bc_xm.getUnitDefinition(), sm)));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_xp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueXp);
if (bc_xp != null && (bc_xp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_xp, sm), getIdentifierSubstitutions(bc_xp.getExpression(), bc_xp.getUnitDefinition(), sm)));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_ym = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueYm);
if (bc_ym != null && (bc_ym.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_ym, sm), getIdentifierSubstitutions(bc_ym.getExpression(), bc_ym.getUnitDefinition(), sm)));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_yp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueYp);
if (bc_yp != null && (bc_yp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_yp, sm), getIdentifierSubstitutions(bc_yp.getExpression(), bc_yp.getUnitDefinition(), sm)));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_zm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueZm);
if (bc_zm != null && (bc_zm.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_zm, sm), getIdentifierSubstitutions(bc_zm.getExpression(), bc_zm.getUnitDefinition(), sm)));
}
SpeciesContextSpec.SpeciesContextSpecParameter bc_zp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueZp);
if (bc_zp != null && (bc_zp.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_zp, sm), getIdentifierSubstitutions(bc_zp.getExpression(), bc_zp.getUnitDefinition(), sm)));
}
}
//
for (int i = 0; i < speciesContextSpecs.length; i++) {
SpeciesContextSpec.SpeciesContextSpecParameter advection_velX = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityX);
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
if (advection_velX != null && (advection_velX.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velX, sm), getIdentifierSubstitutions(advection_velX.getExpression(), advection_velX.getUnitDefinition(), sm)));
}
SpeciesContextSpec.SpeciesContextSpecParameter advection_velY = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityY);
if (advection_velY != null && (advection_velY.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velY, sm), getIdentifierSubstitutions(advection_velY.getExpression(), advection_velY.getUnitDefinition(), sm)));
}
SpeciesContextSpec.SpeciesContextSpecParameter advection_velZ = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityZ);
if (advection_velZ != null && (advection_velZ.getExpression() != null)) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velZ, sm), getIdentifierSubstitutions(advection_velZ.getExpression(), advection_velZ.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());
}
}
//
// conversion factors
//
varHash.addVariable(new Constant(getMathSymbol(model.getKMOLE(), null), getIdentifierSubstitutions(model.getKMOLE().getExpression(), model.getKMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(model.getN_PMOLE().getName(), getIdentifierSubstitutions(model.getN_PMOLE().getExpression(), model.getN_PMOLE().getUnitDefinition(), null)));
varHash.addVariable(new Constant(model.getKMILLIVOLTS().getName(), getIdentifierSubstitutions(model.getKMILLIVOLTS().getExpression(), model.getKMILLIVOLTS().getUnitDefinition(), null)));
varHash.addVariable(new Constant(model.getK_GHK().getName(), getIdentifierSubstitutions(model.getK_GHK().getExpression(), model.getK_GHK().getUnitDefinition(), null)));
//
// geometric functions
//
ModelUnitSystem modelUnitSystem = simContext.getModel().getUnitSystem();
VCUnitDefinition lengthInverseUnit = modelUnitSystem.getLengthUnit().getInverse();
for (int i = 0; i < structureMappings.length; i++) {
StructureMapping sm = structureMappings[i];
Parameter parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumeFraction);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm)));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_SurfaceToVolumeRatio);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm)));
}
if (sm instanceof MembraneMapping && !getResolved(sm)) {
MembraneMapping mm = (MembraneMapping) sm;
parm = ((MembraneMapping) sm).getVolumeFractionParameter();
if (parm.getExpression() == null) {
throw new MappingException("volume fraction not specified for feature '" + structTopology.getInsideFeature(mm.getMembrane()).getName() + "', please refer to Structure Mapping in Application '" + simContext.getName() + "'");
}
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), modelUnitSystem.getInstance_DIMENSIONLESS(), sm)));
parm = mm.getSurfaceToVolumeParameter();
if (parm.getExpression() == null) {
throw new MappingException("surface to volume ratio not specified for membrane '" + mm.getMembrane().getName() + "', please refer to Structure Mapping in Application '" + simContext.getName() + "'");
}
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), lengthInverseUnit, sm)));
}
StructureMappingParameter sizeParm = sm.getSizeParameter();
if (sizeParm != null) {
if (simContext.getGeometry().getDimension() == 0) {
if (sizeParm.getExpression() != null) {
try {
double value = sizeParm.getExpression().evaluateConstant();
varHash.addVariable(new Constant(getMathSymbol(sizeParm, sm), new Expression(value)));
} catch (ExpressionException e) {
// varHash.addVariable(new Function(getMathSymbol(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.");
}
}
} else {
String compartmentName = null;
VCUnitDefinition sizeUnit = sm.getSizeParameter().getUnitDefinition();
String sizeFunctionName = null;
if (sm instanceof MembraneMapping) {
MembraneMapping mm = (MembraneMapping) sm;
if (getResolved(mm)) {
FeatureMapping fm_inside = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(structTopology.getInsideFeature(mm.getMembrane()));
FeatureMapping fm_outside = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(structTopology.getOutsideFeature(mm.getMembrane()));
compartmentName = getSubVolume(fm_inside).getName() + "_" + getSubVolume(fm_outside).getName();
sizeFunctionName = MathFunctionDefinitions.Function_regionArea_current.getFunctionName();
} else {
FeatureMapping fm_inside = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(structTopology.getInsideFeature(mm.getMembrane()));
FeatureMapping fm_outside = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(structTopology.getOutsideFeature(mm.getMembrane()));
if (getSubVolume(fm_inside) == getSubVolume(fm_outside)) {
compartmentName = getSubVolume(fm_inside).getName();
sizeFunctionName = MathFunctionDefinitions.Function_regionVolume_current.getFunctionName();
} else {
throw new RuntimeException("unexpected structure mapping for membrane '" + mm.getMembrane().getName() + "'");
}
}
} else if (sm instanceof FeatureMapping) {
FeatureMapping fm = (FeatureMapping) sm;
compartmentName = getSubVolume(fm).getName();
sizeFunctionName = MathFunctionDefinitions.Function_regionVolume_current.getFunctionName();
} else {
throw new RuntimeException("structure mapping " + sm.getClass().getName() + " not yet supported");
}
Expression totalVolumeCorrection = sm.getStructureSizeCorrection(simContext, this);
Expression sizeFunctionExpression = Expression.function(sizeFunctionName, new Expression[] { new Expression("'" + compartmentName + "'") });
sizeFunctionExpression.bindExpression(mathDesc);
varHash.addVariable(newFunctionOrConstant(getMathSymbol(sizeParm, sm), getIdentifierSubstitutions(Expression.mult(totalVolumeCorrection, sizeFunctionExpression), sizeUnit, sm)));
parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitArea);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm)));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitVolume);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm)));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitArea);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm)));
}
parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitVolume);
if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm)));
}
}
}
}
//
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
varHash.addVariable(newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], null), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), null)));
}
//
// functions
//
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());
varHash.addVariable(newFunctionOrConstant(getMathSymbol(scm.getSpeciesContext(), sm), getIdentifierSubstitutions(scm.getDependencyExpression(), scm.getSpeciesContext().getUnitDefinition(), sm)));
}
}
//
// set Variables to MathDescription all at once with the order resolved by "VariableHash"
//
mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
//
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");
}
//
// volume subdomains
//
subVolumes = simContext.getGeometryContext().getGeometry().getGeometrySpec().getSubVolumes();
VCUnitDefinition timeUnit = modelUnitSystem.getTimeUnit();
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;
}
//
// create subDomain
//
CompartmentSubDomain subDomain = new CompartmentSubDomain(subVolume.getName(), priority);
mathDesc.addSubDomain(subDomain);
//
if (spatialFeature != null) {
FeatureMapping fm = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(spatialFeature);
subDomain.setBoundaryConditionXm(fm.getBoundaryConditionTypeXm());
subDomain.setBoundaryConditionXp(fm.getBoundaryConditionTypeXp());
if (simContext.getGeometry().getDimension() > 1) {
subDomain.setBoundaryConditionYm(fm.getBoundaryConditionTypeYm());
subDomain.setBoundaryConditionYp(fm.getBoundaryConditionTypeYp());
}
if (simContext.getGeometry().getDimension() > 2) {
subDomain.setBoundaryConditionZm(fm.getBoundaryConditionTypeZm());
subDomain.setBoundaryConditionZp(fm.getBoundaryConditionTypeZp());
}
}
//
// create equations
//
VolumeStructureAnalyzer structureAnalyzer = getVolumeStructureAnalyzer(subVolume);
Enumeration<SpeciesContextMapping> enumSCM = getSpeciesContextMappings();
while (enumSCM.hasMoreElements()) {
SpeciesContextMapping scm = enumSCM.nextElement();
//
if (scm.getVariable() instanceof VolVariable && scm.getDependencyExpression() == null) {
SpeciesContext sc = scm.getSpeciesContext();
StructureMapping sm = simContext.getGeometryContext().getStructureMapping(sc.getStructure());
SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(sc);
VolVariable variable = (VolVariable) scm.getVariable();
Equation equation = null;
if ((scm.isPDERequired()) && sm instanceof FeatureMapping) {
//
if (getSubVolume((FeatureMapping) sm) == subVolume) {
//
// species context belongs to this subDomain
//
Expression initial = new Expression(getMathSymbol(scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration), sm));
Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()));
Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), sm));
equation = new PdeEquation(variable, initial, rate, diffusion);
((PdeEquation) equation).setBoundaryXm((scs.getBoundaryXmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXmParameter(), sm)));
((PdeEquation) equation).setBoundaryXp((scs.getBoundaryXpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXpParameter(), sm)));
((PdeEquation) equation).setBoundaryYm((scs.getBoundaryYmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYmParameter(), sm)));
((PdeEquation) equation).setBoundaryYp((scs.getBoundaryYpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYpParameter(), sm)));
((PdeEquation) equation).setBoundaryZm((scs.getBoundaryZmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZmParameter(), sm)));
((PdeEquation) equation).setBoundaryZp((scs.getBoundaryZpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZpParameter(), sm)));
((PdeEquation) equation).setVelocityX((scs.getVelocityXParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getVelocityXParameter(), sm)));
((PdeEquation) equation).setVelocityY((scs.getVelocityYParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getVelocityYParameter(), sm)));
((PdeEquation) equation).setVelocityZ((scs.getVelocityZParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getVelocityZParameter(), sm)));
subDomain.replaceEquation(equation);
} else {
Expression initial = new Expression(0.0);
Expression rate = new Expression(0.0);
Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), sm));
equation = new PdeEquation(variable, initial, rate, diffusion);
if (subDomain.getEquation(variable) == null) {
subDomain.addEquation(equation);
}
}
} else {
//
// ODE
//
SubVolume mappedSubVolume = null;
if (sm instanceof FeatureMapping) {
mappedSubVolume = getSubVolume((FeatureMapping) sm);
} else if (sm instanceof MembraneMapping) {
// membrane is mapped to that of the inside feature
FeatureMapping featureMapping = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(structTopology.getInsideFeature((Membrane) sm.getStructure()));
mappedSubVolume = getSubVolume(featureMapping);
}
if (mappedSubVolume == subVolume) {
//
// species context belongs to this subDomain
//
Expression initial = new Expression(getMathSymbol(scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration), null));
Expression rate = (scm.getRate() == null) ? new Expression(0.0) : getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()));
equation = new OdeEquation(variable, initial, rate);
subDomain.replaceEquation(equation);
} else {
Expression initial = new Expression(0.0);
Expression rate = new Expression(0.0);
equation = new OdeEquation(variable, initial, rate);
if (subDomain.getEquation(variable) == null) {
subDomain.addEquation(equation);
}
}
}
}
}
//
// create fast system (if neccessary)
//
SpeciesContextMapping[] fastSpeciesContextMappings = structureAnalyzer.getFastSpeciesContextMappings();
VCUnitDefinition subDomainUnit = modelUnitSystem.getVolumeConcentrationUnit();
if (fastSpeciesContextMappings != null) {
FastSystem fastSystem = new FastSystem(mathDesc);
for (int i = 0; i < fastSpeciesContextMappings.length; i++) {
SpeciesContextMapping scm = fastSpeciesContextMappings[i];
if (scm.getFastInvariant() == null) {
//
// independant-fast variable, create a fastRate object
//
Expression rate = getIdentifierSubstitutions(scm.getFastRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(getResolvedFeature(subVolume)));
FastRate fastRate = new FastRate(rate);
fastSystem.addFastRate(fastRate);
} else {
//
// dependant-fast variable, create a fastInvariant object
//
Expression rate = getIdentifierSubstitutions(scm.getFastInvariant(), subDomainUnit, simContext.getGeometryContext().getStructureMapping(getResolvedFeature(subVolume)));
FastInvariant fastInvariant = new FastInvariant(rate);
fastSystem.addFastInvariant(fastInvariant);
}
}
subDomain.setFastSystem(fastSystem);
// constructor calls the 'refresh' method which constructs depemdency matrix, dependent/independent vars and pseudoconstants, etc.
FastSystemAnalyzer fs_analyzer = new FastSystemAnalyzer(fastSystem, mathDesc);
}
//
// create ode's for voltages to be calculated on unresolved membranes mapped to this subVolume
//
Structure[] localStructures = getStructures(subVolume);
for (int sIndex = 0; sIndex < localStructures.length; sIndex++) {
if (localStructures[sIndex] instanceof Membrane) {
Membrane membrane = (Membrane) localStructures[sIndex];
MembraneMapping membraneMapping = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(membrane);
if (!getResolved(membraneMapping) && membraneMapping.getCalculateVoltage()) {
MembraneElectricalDevice capacitiveDevice = potentialMapping.getCapacitiveDevice(membrane);
if (capacitiveDevice.getDependentVoltageExpression() == null) {
VolVariable vVar = (VolVariable) mathDesc.getVariable(getMathSymbol(capacitiveDevice.getVoltageSymbol(), membraneMapping));
Expression initExp = new Expression(getMathSymbol(capacitiveDevice.getMembraneMapping().getInitialVoltageParameter(), membraneMapping));
subDomain.addEquation(new OdeEquation(vVar, initExp, getIdentifierSubstitutions(potentialMapping.getOdeRHS(capacitiveDevice, this), membrane.getMembraneVoltage().getUnitDefinition().divideBy(timeUnit), membraneMapping)));
} else {
//
//
//
}
}
}
}
}
//
for (int k = 0; k < subVolumes.length; k++) {
SubVolume subVolume = (SubVolume) subVolumes[k];
//
// if there is a spatially resolved membrane surrounding this subVolume, then create a membraneSubDomain
//
structures = getStructures(subVolume);
Membrane membrane = null;
if (structures != null) {
for (int j = 0; j < structures.length; j++) {
if (structures[j] instanceof Membrane && getResolved(simContext.getGeometryContext().getStructureMapping(structures[j]))) {
membrane = (Membrane) structures[j];
}
}
}
if (membrane == null) {
continue;
}
SubVolume outerSubVolume = getSubVolume(((FeatureMapping) simContext.getGeometryContext().getStructureMapping(structTopology.getOutsideFeature(membrane))));
SubVolume innerSubVolume = getSubVolume(((FeatureMapping) simContext.getGeometryContext().getStructureMapping(structTopology.getInsideFeature(membrane))));
if (innerSubVolume != subVolume) {
throw new MappingException("membrane " + membrane.getName() + " improperly mapped to inner subVolume " + innerSubVolume.getName());
}
//
// get priority of subDomain
//
// Feature spatialFeature = simContext.getGeometryContext().getResolvedFeature(subVolume);
// int priority = spatialFeature.getPriority();
//
// create subDomain
//
CompartmentSubDomain outerCompartment = mathDesc.getCompartmentSubDomain(outerSubVolume.getName());
CompartmentSubDomain innerCompartment = mathDesc.getCompartmentSubDomain(innerSubVolume.getName());
SurfaceClass surfaceClass = simContext.getGeometry().getGeometrySurfaceDescription().getSurfaceClass(innerSubVolume, outerSubVolume);
MembraneSubDomain memSubDomain = new MembraneSubDomain(innerCompartment, outerCompartment, surfaceClass.getName());
mathDesc.addSubDomain(memSubDomain);
//
// create equations for membrane-bound molecular species
//
MembraneStructureAnalyzer membraneStructureAnalyzer = getMembraneStructureAnalyzer(membrane);
Enumeration<SpeciesContextMapping> enumSCM = getSpeciesContextMappings();
while (enumSCM.hasMoreElements()) {
SpeciesContextMapping scm = enumSCM.nextElement();
SpeciesContext sc = scm.getSpeciesContext();
SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(sc);
//
if ((scm.getVariable() instanceof MemVariable) && scm.getDependencyExpression() == null) {
//
// independant variable, create an equation object
//
Equation equation = null;
MemVariable variable = (MemVariable) scm.getVariable();
MembraneMapping mm = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(sc.getStructure());
if (scm.isPDERequired()) {
//
if (mm.getMembrane() == membrane) {
//
// species context belongs to this subDomain
//
Expression initial = new Expression(getMathSymbol(scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration), mm));
Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()));
Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), mm));
equation = new PdeEquation(variable, initial, rate, diffusion);
((PdeEquation) equation).setBoundaryXm((scs.getBoundaryXmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXmParameter(), mm)));
((PdeEquation) equation).setBoundaryXp((scs.getBoundaryXpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXpParameter(), mm)));
((PdeEquation) equation).setBoundaryYm((scs.getBoundaryYmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYmParameter(), mm)));
((PdeEquation) equation).setBoundaryYp((scs.getBoundaryYpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYpParameter(), mm)));
((PdeEquation) equation).setBoundaryZm((scs.getBoundaryZmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZmParameter(), mm)));
((PdeEquation) equation).setBoundaryZp((scs.getBoundaryZpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZpParameter(), mm)));
memSubDomain.replaceEquation(equation);
} else {
Expression initial = new Expression(0.0);
Expression rate = new Expression(0.0);
Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), mm));
equation = new PdeEquation(variable, initial, rate, diffusion);
if (memSubDomain.getEquation(variable) == null) {
memSubDomain.addEquation(equation);
}
}
} else {
//
if (mm.getMembrane() == membrane) {
//
// species context belongs to this subDomain
//
Expression initial = new Expression(getMathSymbol(scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration), null));
Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()));
equation = new OdeEquation(variable, initial, rate);
memSubDomain.replaceEquation(equation);
} else {
Expression initial = new Expression(0.0);
Expression rate = new Expression(0.0);
equation = new OdeEquation(variable, initial, rate);
if (memSubDomain.getEquation(variable) == null) {
memSubDomain.addEquation(equation);
}
}
}
}
}
//
// create dummy jump conditions for all volume variables that diffuse and/or advect
//
Enumeration<SpeciesContextMapping> enum_scm = getSpeciesContextMappings();
while (enum_scm.hasMoreElements()) {
SpeciesContextMapping scm = enum_scm.nextElement();
if (scm.isPDERequired()) {
// Species species = scm.getSpeciesContext().getSpecies();
Variable var = scm.getVariable();
if (var instanceof VolVariable && (scm.isPDERequired())) {
JumpCondition jc = memSubDomain.getJumpCondition((VolVariable) var);
if (jc == null) {
// System.out.println("MathMapping.refreshMathDescription(), adding jump condition for diffusing variable "+var.getName()+" on membrane "+membraneStructureAnalyzer.getMembrane().getName());
jc = new JumpCondition((VolVariable) var);
memSubDomain.addJumpCondition(jc);
}
}
}
}
//
// create jump conditions for any volume variables that bind to membrane or have explicitly defined fluxes
//
ResolvedFlux[] resolvedFluxes = membraneStructureAnalyzer.getResolvedFluxes();
if (resolvedFluxes != null) {
for (int i = 0; i < resolvedFluxes.length; i++) {
Species species = resolvedFluxes[i].getSpecies();
SpeciesContext sc = simContext.getReactionContext().getModel().getSpeciesContext(species, structTopology.getInsideFeature(membraneStructureAnalyzer.getMembrane()));
if (sc == null) {
sc = simContext.getReactionContext().getModel().getSpeciesContext(species, structTopology.getOutsideFeature(membraneStructureAnalyzer.getMembrane()));
}
SpeciesContextMapping scm = getSpeciesContextMapping(sc);
// if (scm.getVariable() instanceof VolVariable && scm.isDiffusing()){
if (scm.getVariable() instanceof VolVariable && ((MembraneStructureAnalyzer.bNoFluxIfFixed || (scm.isPDERequired())))) {
if (MembraneStructureAnalyzer.bNoFluxIfFixed && !scm.isPDERequired()) {
MembraneStructureAnalyzer.bNoFluxIfFixedExercised = true;
}
JumpCondition jc = memSubDomain.getJumpCondition((VolVariable) scm.getVariable());
if (jc == null) {
jc = new JumpCondition((VolVariable) scm.getVariable());
memSubDomain.addJumpCondition(jc);
}
Expression inFlux = getIdentifierSubstitutions(resolvedFluxes[i].inFluxExpression, resolvedFluxes[i].getUnitDefinition(), simContext.getGeometryContext().getStructureMapping(membraneStructureAnalyzer.getMembrane()));
jc.setInFlux(inFlux);
Expression outFlux = getIdentifierSubstitutions(resolvedFluxes[i].outFluxExpression, resolvedFluxes[i].getUnitDefinition(), simContext.getGeometryContext().getStructureMapping(membraneStructureAnalyzer.getMembrane()));
jc.setOutFlux(outFlux);
} else {
throw new MappingException("APPLICATION " + simContext.getName() + " : " + scm.getSpeciesContext().getName() + " has spatially resolved flux at membrane " + membrane.getName() + ", but doesn't diffuse in compartment " + scm.getSpeciesContext().getStructure().getName());
}
}
}
//
// create fast system (if neccessary)
//
SpeciesContextMapping[] fastSpeciesContextMappings = membraneStructureAnalyzer.getFastSpeciesContextMappings();
if (fastSpeciesContextMappings != null) {
FastSystem fastSystem = new FastSystem(mathDesc);
for (int i = 0; i < fastSpeciesContextMappings.length; i++) {
SpeciesContextMapping scm = fastSpeciesContextMappings[i];
if (scm.getFastInvariant() == null) {
//
// independant-fast variable, create a fastRate object
//
VCUnitDefinition rateUnit = scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit);
MembraneMapping membraneMapping = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(membraneStructureAnalyzer.getMembrane());
FastRate fastRate = new FastRate(getIdentifierSubstitutions(scm.getFastRate(), rateUnit, membraneMapping));
fastSystem.addFastRate(fastRate);
} else {
//
// dependant-fast variable, create a fastInvariant object
//
VCUnitDefinition invariantUnit = scm.getSpeciesContext().getUnitDefinition();
MembraneMapping membraneMapping = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(membraneStructureAnalyzer.getMembrane());
FastInvariant fastInvariant = new FastInvariant(getIdentifierSubstitutions(scm.getFastInvariant(), invariantUnit, membraneMapping));
fastSystem.addFastInvariant(fastInvariant);
}
}
memSubDomain.setFastSystem(fastSystem);
// constructor calls the 'refresh' method which constructs depemdency matrix, dependent/independent vars and pseudoconstants, etc.
FastSystemAnalyzer fs_analyzer = new FastSystemAnalyzer(fastSystem, mathDesc);
}
//
// create Membrane-region equations for potential of this resolved membrane
//
MembraneMapping membraneMapping = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(membrane);
if (membraneMapping.getCalculateVoltage()) {
ElectricalDevice[] membraneDevices = potentialMapping.getElectricalDevices(membrane);
int numCapacitiveDevices = 0;
MembraneElectricalDevice capacitiveDevice = null;
for (int i = 0; i < membraneDevices.length; i++) {
if (membraneDevices[i] instanceof MembraneElectricalDevice) {
numCapacitiveDevices++;
capacitiveDevice = (MembraneElectricalDevice) membraneDevices[i];
}
}
if (numCapacitiveDevices != 1) {
throw new MappingException("expecting 1 capacitive electrical device on graph edge for membrane " + membrane.getName() + ", found '" + numCapacitiveDevices + "'");
}
if (mathDesc.getVariable(getMathSymbol(capacitiveDevice.getVoltageSymbol(), membraneMapping)) instanceof MembraneRegionVariable) {
MembraneRegionVariable vVar = (MembraneRegionVariable) mathDesc.getVariable(getMathSymbol(capacitiveDevice.getVoltageSymbol(), membraneMapping));
Parameter initialVoltageParm = capacitiveDevice.getMembraneMapping().getInitialVoltageParameter();
Expression initExp = getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), capacitiveDevice.getMembraneMapping());
MembraneRegionEquation vEquation = new MembraneRegionEquation(vVar, initExp);
vEquation.setMembraneRateExpression(getIdentifierSubstitutions(potentialMapping.getOdeRHS(capacitiveDevice, this), membrane.getMembraneVoltage().getUnitDefinition().divideBy(timeUnit), capacitiveDevice.getMembraneMapping()));
memSubDomain.addEquation(vEquation);
}
}
}
// create equations for event assign targets that are model params/strutureSize, etc.
Set<VolVariable> hashKeySet = eventVolVarHash.keySet();
Iterator<VolVariable> volVarsIter = hashKeySet.iterator();
// working under teh assumption that we are dealing with non-spatial math, hence only one compartment domain!
SubDomain subDomain = mathDesc.getSubDomains().nextElement();
while (volVarsIter.hasNext()) {
VolVariable volVar = volVarsIter.next();
EventAssignmentInitParameter eap = eventVolVarHash.get(volVar);
Expression rateExpr = new Expression(0.0);
Equation equation = new OdeEquation(volVar, new Expression(getMathSymbol(eap, null)), rateExpr);
subDomain.addEquation(equation);
}
// events - add events to math desc and odes for event assignments that have parameters as target variables
BioEvent[] bioevents = simContext.getBioEvents();
if (bioevents != null && bioevents.length > 0) {
for (BioEvent be : bioevents) {
// transform the bioEvent trigger/delay to math Event
Expression mathTriggerExpr = getIdentifierSubstitutions(be.generateTriggerExpression(), modelUnitSystem.getInstance_DIMENSIONLESS(), null);
Delay mathDelay = null;
if (be.getParameter(BioEventParameterType.TriggerDelay) != null) {
boolean bUseValsFromTriggerTime = be.getUseValuesFromTriggerTime();
Expression mathDelayExpr = getIdentifierSubstitutions(be.getParameter(BioEventParameterType.TriggerDelay).getExpression(), timeUnit, null);
mathDelay = new Delay(bUseValsFromTriggerTime, mathDelayExpr);
}
// now deal with (bio)event Assignment translation to math EventAssignment
ArrayList<EventAssignment> eventAssignments = be.getEventAssignments();
ArrayList<Event.EventAssignment> mathEventAssignmentsList = new ArrayList<Event.EventAssignment>();
for (EventAssignment ea : eventAssignments) {
SymbolTableEntry ste = simContext.getEntry(ea.getTarget().getName());
VCUnitDefinition eventAssignVarUnit = ste.getUnitDefinition();
Variable variable = varHash.getVariable(ste.getName());
Event.EventAssignment mathEA = new Event.EventAssignment(variable, getIdentifierSubstitutions(ea.getAssignmentExpression(), eventAssignVarUnit, null));
mathEventAssignmentsList.add(mathEA);
}
// use the translated trigger, delay and event assignments to create (math) event
Event mathEvent = new Event(be.getName(), mathTriggerExpr, mathDelay, mathEventAssignmentsList);
mathDesc.addEvent(mathEvent);
}
}
if (!mathDesc.isValid()) {
throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
}
// System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string begin ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
// System.out.println(mathDesc.getVCML());
// System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string end ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
}
Also used :
MembraneMapping(cbit.vcell.mapping.MembraneMapping)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
ArrayList(java.util.ArrayList)
SpeciesContext(cbit.vcell.model.SpeciesContext)
StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter)
Feature(cbit.vcell.model.Feature)
MemVariable(cbit.vcell.math.MemVariable)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
SubVolume(cbit.vcell.geometry.SubVolume)
Vector(java.util.Vector)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
StructureTopology(cbit.vcell.model.Model.StructureTopology)
ReactionSpec(cbit.vcell.mapping.ReactionSpec)
FastInvariant(cbit.vcell.math.FastInvariant)
PropertyVetoException(java.beans.PropertyVetoException)
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
FastSystem(cbit.vcell.math.FastSystem)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
ReactionStep(cbit.vcell.model.ReactionStep)
MembraneRegionEquation(cbit.vcell.math.MembraneRegionEquation)
SurfaceClass(cbit.vcell.geometry.SurfaceClass)
VariableHash(cbit.vcell.math.VariableHash)
StructureMapping(cbit.vcell.mapping.StructureMapping)
FeatureMapping(cbit.vcell.mapping.FeatureMapping)
Structure(cbit.vcell.model.Structure)
ModelUnitSystem(cbit.vcell.model.ModelUnitSystem)
Hashtable(java.util.Hashtable)
StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter)
Expression(cbit.vcell.parser.Expression)
Model(cbit.vcell.model.Model)
ProxyParameter(cbit.vcell.model.ProxyParameter)
StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter)
Parameter(cbit.vcell.model.Parameter)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter)
LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
BioEvent(cbit.vcell.mapping.BioEvent)
Event(cbit.vcell.math.Event)
BioEvent(cbit.vcell.mapping.BioEvent)
MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable)
MathDescription(cbit.vcell.math.MathDescription)
SpeciesContextMapping(cbit.vcell.mapping.SpeciesContextMapping)
SpeciesContextSpec(cbit.vcell.mapping.SpeciesContextSpec)
ExpressionException(cbit.vcell.parser.ExpressionException)
Delay(cbit.vcell.math.Event.Delay)
MappingException(cbit.vcell.mapping.MappingException)
PropertyVetoException(java.beans.PropertyVetoException)
PdeEquation(cbit.vcell.math.PdeEquation)
CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain)
SubDomain(cbit.vcell.math.SubDomain)
MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)
Species(cbit.vcell.model.Species)
SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)
EventAssignment(cbit.vcell.mapping.BioEvent.EventAssignment)
VolVariable(cbit.vcell.math.VolVariable)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
OdeEquation(cbit.vcell.math.OdeEquation)
JumpCondition(cbit.vcell.math.JumpCondition)
MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable)
VolVariable(cbit.vcell.math.VolVariable)
MemVariable(cbit.vcell.math.MemVariable)
Variable(cbit.vcell.math.Variable)
SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter)
Constant(cbit.vcell.math.Constant)
SymbolTableEntry(cbit.vcell.parser.SymbolTableEntry)
Membrane(cbit.vcell.model.Membrane)
EventAssignment(cbit.vcell.mapping.BioEvent.EventAssignment)
OdeEquation(cbit.vcell.math.OdeEquation)
MembraneRegionEquation(cbit.vcell.math.MembraneRegionEquation)
PdeEquation(cbit.vcell.math.PdeEquation)
Equation(cbit.vcell.math.Equation)
FastRate(cbit.vcell.math.FastRate)
Aggregations
VariableHash (cbit.vcell.math.VariableHash)16
Expression (cbit.vcell.parser.Expression)15
Variable (cbit.vcell.math.Variable)14
CompartmentSubDomain (cbit.vcell.math.CompartmentSubDomain)12
Constant (cbit.vcell.math.Constant)12
Function (cbit.vcell.math.Function)12
VolVariable (cbit.vcell.math.VolVariable)11
MathDescription (cbit.vcell.math.MathDescription)10
ExpressionException (cbit.vcell.parser.ExpressionException)10
Domain (cbit.vcell.math.Variable.Domain)9
PropertyVetoException (java.beans.PropertyVetoException)9
Vector (java.util.Vector)9
MathException (cbit.vcell.math.MathException)8
MemVariable (cbit.vcell.math.MemVariable)7
MembraneRegionVariable (cbit.vcell.math.MembraneRegionVariable)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
SubVolume (cbit.vcell.geometry.SubVolume)6
