Example 36 with VCUnitDefinition
use of cbit.vcell.units.VCUnitDefinition in project vcell by virtualcell.
the class PotentialMapping method determineLumpedEquations.
/**
* Insert the method's description here.
* Creation date: (2/20/2002 9:02:42 AM)
* @return cbit.vcell.mathmodel.MathModel
* @param circuitGraph cbit.vcell.mapping.potential.Graph
*/
private void determineLumpedEquations(Graph graph, double temperatureKelvin) throws ExpressionException, MatrixException, MappingException, MathException {
//
// traverses graph and calculates RHS expressions for all capacitive devices (dV/dt)
//
// calculate dependent voltages as functions of independent voltages.
//
//
Graph[] spanningTrees = graph.getSpanningForest();
if (!bSilent) {
System.out.println("spanning tree(s):");
for (int i = 0; i < spanningTrees.length; i++) {
System.out.println(i + ") " + spanningTrees[i]);
}
}
Path[] fundamentalCycles = graph.getFundamentalCycles();
if (!bSilent) {
System.out.println("fundamental cycles:");
for (int i = 0; i < fundamentalCycles.length; i++) {
System.out.println(" " + fundamentalCycles[i]);
}
}
//
// print out basic device information
//
fieldEdges = graph.getEdges();
//
if (!bSilent) {
System.out.println("\n\n applying KVL to <<ALL>> fundamental cycles\n");
}
Path[] graphCycles = graph.getFundamentalCycles();
RationalExpMatrix voltageMatrix = new RationalExpMatrix(graphCycles.length, graph.getNumEdges());
for (int i = 0; i < graphCycles.length; i++) {
Edge[] cycleEdges = graphCycles[i].getEdges();
Node[] nodesTraversed = graphCycles[i].getNodesTraversed();
Expression exp = new Expression(0.0);
//
for (int j = 0; j < cycleEdges.length; j++) {
int edgeIndex = graph.getIndex(cycleEdges[j]);
Expression voltage = new Expression(((ElectricalDevice) cycleEdges[j].getData()).getVoltageSymbol(), fieldMathMapping.getNameScope());
if (cycleEdges[j].getNode1().equals(nodesTraversed[j])) {
// going same direction
exp = Expression.add(exp, voltage);
voltageMatrix.set_elem(i, edgeIndex, 1);
} else {
// going opposite direction
exp = Expression.add(exp, Expression.negate(voltage));
voltageMatrix.set_elem(i, edgeIndex, -1);
}
}
if (!bSilent) {
System.out.println(exp.flatten().infix() + " = 0.0");
}
}
if (!bSilent) {
voltageMatrix.show();
}
if (voltageMatrix.getNumRows() > 0) {
RationalExpMatrix vPermutationMatrix = new RationalExpMatrix(voltageMatrix.getNumRows(), voltageMatrix.getNumRows());
voltageMatrix.gaussianElimination(vPermutationMatrix);
if (!bSilent) {
System.out.println("reduced matrix");
voltageMatrix.show();
}
} else {
voltageMatrix = null;
}
//
// declare dependent voltages as functions of independent voltages
//
// 1) Always use Voltage-Clamps as independent voltages
// 2) Try to use Capacitive devices as independent voltages
//
// solve for current-clamp voltages and redundant/constrained capacitive voltages as function of (1) and (2).
//
Expression[] dependentVoltageExpressions = new Expression[fieldEdges.length];
//
// make sure assumptions hold regarding edge ordering, otherwise wrong dependent voltages will be selected
//
verifyEdgeOrdering();
for (int i = 0; voltageMatrix != null && i < voltageMatrix.getNumRows(); i++) {
//
// find first '1.0' element, this column is the next 'dependent' voltage
//
int column = -1;
for (int j = i; j < voltageMatrix.getNumCols(); j++) {
RationalExp elem = voltageMatrix.get(i, j);
if (elem.isConstant() && elem.getConstant().doubleValue() == 1.0) {
column = j;
break;
}
}
if (column != -1) {
//
// get electrical device of dependent voltage
//
ElectricalDevice device = (ElectricalDevice) fieldEdges[column].getData();
//
// form dependency expression
//
StringBuffer buffer = new StringBuffer();
for (int j = column + 1; j < graph.getNumEdges(); j++) {
if (!voltageMatrix.get_elem(i, j).isZero()) {
ElectricalDevice colDevice = (ElectricalDevice) fieldEdges[j].getData();
Expression voltage = new Expression(colDevice.getVoltageSymbol(), fieldMathMapping.getNameScope());
buffer.append(" + " + voltageMatrix.get(i, j).minus().infixString() + '*' + voltage.infix());
}
}
dependentVoltageExpressions[column] = (new Expression(buffer.toString())).flatten();
dependentVoltageExpressions[column].bindExpression(device);
device.setDependentVoltageExpression(dependentVoltageExpressions[column]);
}
}
if (!bSilent) {
for (int i = 0; i < dependentVoltageExpressions.length; i++) {
System.out.println("dependentVoltageExpressions[" + i + "] = " + dependentVoltageExpressions[i]);
}
}
//
if (!bSilent) {
System.out.println("\n\nSOLVE FOR TOTAL CURRENTS IN TERMS OF APPLIED CURRENTS");
System.out.println("\n\n 1) applying KVL to all fundamental \"capacitive\" and Voltage-clamp cycles (after current-clamp edges are removed)\n");
}
Graph capacitorGraph = new Graph();
for (int i = 0; i < fieldEdges.length; i++) {
ElectricalDevice device = (ElectricalDevice) fieldEdges[i].getData();
if (device.hasCapacitance() || device.isVoltageSource()) {
capacitorGraph.addEdge(fieldEdges[i]);
}
}
Path[] capacitorGraphCycles = capacitorGraph.getFundamentalCycles();
RationalExpMatrix capacitorGraphVoltageMatrix = new RationalExpMatrix(capacitorGraphCycles.length, graph.getNumEdges());
for (int i = 0; i < capacitorGraphCycles.length; i++) {
Edge[] cycleEdges = capacitorGraphCycles[i].getEdges();
Node[] nodesTraversed = capacitorGraphCycles[i].getNodesTraversed();
Expression exp = new Expression(0.0);
//
for (int j = 0; j < cycleEdges.length; j++) {
int edgeIndex = graph.getIndex(cycleEdges[j]);
Expression voltage = new Expression(((ElectricalDevice) cycleEdges[j].getData()).getVoltageSymbol(), fieldMathMapping.getNameScope());
if (cycleEdges[j].getNode1().equals(nodesTraversed[j])) {
// going same direction
exp = Expression.add(exp, voltage);
capacitorGraphVoltageMatrix.set_elem(i, edgeIndex, 1);
} else {
// going opposite direction
exp = Expression.add(exp, Expression.negate(voltage));
capacitorGraphVoltageMatrix.set_elem(i, edgeIndex, -1);
}
}
if (!bSilent) {
System.out.println(exp.flatten().infix() + " = 0.0");
}
}
if (!bSilent) {
capacitorGraphVoltageMatrix.show();
}
//
if (!bSilent) {
System.out.println("\n\n 2) applying KCL to all nodes (except one) .. n-1 nodes of full graph --- CONSERVATION OF TOTAL CURRENT\n");
}
Node[] nodes = graph.getNodes();
RationalExpMatrix kclMatrix = new RationalExpMatrix(graph.getNumNodes() - 1, graph.getNumEdges());
for (int i = 0; i < nodes.length - 1; i++) {
if (graph.getDegree(nodes[i]) > 0) {
Edge[] adjacentEdges = graph.getAdjacentEdges(nodes[i]);
Expression exp = new Expression(0.0);
//
for (int j = 0; j < adjacentEdges.length; j++) {
int edgeIndex = graph.getIndex(adjacentEdges[j]);
Expression totalCurrent = new Expression(((ElectricalDevice) adjacentEdges[j].getData()).getTotalCurrentSymbol(), fieldMathMapping.getNameScope());
if (adjacentEdges[j].getNode1().equals(nodes[i])) {
exp = Expression.add(exp, Expression.negate(totalCurrent));
kclMatrix.set_elem(i, edgeIndex, -1);
} else {
exp = Expression.add(exp, totalCurrent);
kclMatrix.set_elem(i, edgeIndex, 1);
}
}
if (!bSilent) {
System.out.println(exp.flatten().infix() + " = 0.0");
}
} else {
throw new MappingException("compartment '" + nodes[i].getName() + "' is electrically isolated, cannot generate ciruit equations for application '" + fieldSimContext.getName() + "'. \n\nPlease specify electrical properties for all membranes (see Structures tab in Physiology).");
}
}
if (!bSilent) {
kclMatrix.show();
}
//
if (!bSilent) {
System.out.println("\n\n 3) form total 'current' matrix\n");
}
int numNonCapacitiveEdges = 0;
for (int i = 0; i < fieldEdges.length; i++) {
ElectricalDevice device = (ElectricalDevice) fieldEdges[i].getData();
if (!device.hasCapacitance()) {
numNonCapacitiveEdges++;
}
}
int cmat_rows = kclMatrix.getNumRows() + capacitorGraphVoltageMatrix.getNumRows() + numNonCapacitiveEdges;
RationalExpMatrix currentMatrix = new RationalExpMatrix(cmat_rows, 3 * graph.getNumEdges());
//
// order edges for current elimination (unknown voltage-clamp currents as well as all total currents).
//
int[] cIndex = new int[fieldEdges.length];
int ci = 0;
for (int i = 0; i < fieldEdges.length; i++) {
if (((ElectricalDevice) fieldEdges[i].getData()).isVoltageSource()) {
cIndex[ci++] = i;
}
}
for (int i = 0; i < fieldEdges.length; i++) {
if (!((ElectricalDevice) fieldEdges[i].getData()).isVoltageSource()) {
cIndex[ci++] = i;
}
}
if (ci != fieldEdges.length) {
throw new RuntimeException("error computing current indexes");
}
if (!bSilent) {
System.out.println("reordered devices for current matrix elimination");
for (int i = 0; i < cIndex.length; i++) {
Expression voltage = new Expression(((ElectricalDevice) fieldEdges[cIndex[i]].getData()).getVoltageSymbol(), fieldMathMapping.getNameScope());
System.out.println(i + ") = device[" + cIndex[i] + "] = " + voltage.infix());
}
}
int row = 0;
//
for (int i = 0; i < kclMatrix.getNumRows(); i++) {
for (int j = 0; j < kclMatrix.getNumCols(); j++) {
// entry for i's
currentMatrix.set_elem(row, j, kclMatrix.get(i, cIndex[j]));
}
row++;
}
//
for (int i = 0; i < fieldEdges.length; i++) {
ElectricalDevice device = (ElectricalDevice) fieldEdges[cIndex[i]].getData();
if (!device.hasCapacitance()) {
currentMatrix.set_elem(row, i, 1);
currentMatrix.set_elem(row, i + graph.getNumEdges(), -1);
row++;
}
}
//
// add current terms of (i - F)*1000/C form using KVL relationships from "capacitor graph"
//
ModelUnitSystem modelUnitSystem = fieldMathMapping.getSimulationContext().getModel().getUnitSystem();
VCUnitDefinition potentialUnit = modelUnitSystem.getVoltageUnit();
VCUnitDefinition timeUnit = modelUnitSystem.getTimeUnit();
for (int i = 0; i < capacitorGraphVoltageMatrix.getNumRows(); i++) {
for (int j = 0; j < capacitorGraphVoltageMatrix.getNumCols(); j++) {
ElectricalDevice device = (ElectricalDevice) fieldEdges[cIndex[j]].getData();
RationalExp coefficient = capacitorGraphVoltageMatrix.get(i, cIndex[j]);
if (device.hasCapacitance()) {
//
// replace dVi/dt with 1000/Ci * Ii + 1000/Ci * Fi
//
SymbolTableEntry capacitanceParameter = ((MembraneElectricalDevice) device).getCapacitanceParameter();
Expression capacitance = new Expression(capacitanceParameter, fieldMathMapping.getNameScope());
String Cname = capacitance.infix();
VCUnitDefinition unitFactor = potentialUnit.divideBy(timeUnit).multiplyBy(capacitanceParameter.getUnitDefinition()).divideBy(device.getTotalCurrentSymbol().getUnitDefinition());
RationalExp unitFactorExp = fieldMathMapping.getUnitFactorAsRationalExp(unitFactor);
// entry for i's
currentMatrix.set_elem(row, j, coefficient.mult(unitFactorExp.div(new RationalExp(Cname))));
// entry for F's
currentMatrix.set_elem(row, j + graph.getNumEdges(), coefficient.minus().mult(unitFactorExp).div(new RationalExp(Cname)));
} else if (device.isVoltageSource()) {
//
// directly insert "symbolic" dVi/dt into the new matrix
//
currentMatrix.set_elem(row, j + 2 * graph.getNumEdges(), coefficient);
}
}
row++;
}
if (!bSilent) {
currentMatrix.show();
}
if (currentMatrix.getNumRows() > 0) {
RationalExpMatrix cPermutationMatrix = new RationalExpMatrix(currentMatrix.getNumRows(), currentMatrix.getNumRows());
currentMatrix.gaussianElimination(cPermutationMatrix);
if (!bSilent) {
System.out.println("reduced matrix");
currentMatrix.show();
}
}
//
if (!bSilent) {
System.out.println("\n\n total currents for each device\n");
}
Expression[] totalCurrents = new Expression[fieldEdges.length];
for (int i = 0; i < fieldEdges.length; i++) {
ElectricalDevice device = (ElectricalDevice) fieldEdges[cIndex[i]].getData();
StringBuffer buffer = new StringBuffer("0.0");
//
for (int j = 0; j < graph.getNumEdges(); j++) {
RationalExp coefficient = currentMatrix.get(i, j + graph.getNumEdges());
if (!coefficient.isZero()) {
ElectricalDevice colDevice = (ElectricalDevice) fieldEdges[cIndex[j]].getData();
Expression source = new Expression(colDevice.getSourceSymbol(), fieldMathMapping.getNameScope());
buffer.append(" + " + coefficient.minus() + "*" + source.infix());
}
}
//
for (int j = 0; j < graph.getNumEdges(); j++) {
RationalExp coefficient = currentMatrix.get(i, j + 2 * graph.getNumEdges());
if (!coefficient.isZero()) {
VoltageClampElectricalDevice colDevice = (VoltageClampElectricalDevice) fieldEdges[cIndex[j]].getData();
Expression timeDeriv = colDevice.getVoltageClampStimulus().getVoltageParameter().getExpression().differentiate("t");
timeDeriv = timeDeriv.flatten();
timeDeriv.bindExpression(colDevice);
timeDeriv.renameBoundSymbols(colDevice.getNameScope().getParent());
buffer.append(" + " + coefficient.minus() + "*" + timeDeriv.infix());
}
}
totalCurrents[cIndex[i]] = (new Expression(buffer.toString())).flatten();
totalCurrents[cIndex[i]].bindExpression(device.getNameScope().getParent().getScopedSymbolTable());
device.getParameterFromRole(ElectricalDevice.ROLE_TotalCurrent).setExpression(totalCurrents[cIndex[i]]);
}
if (!bSilent) {
for (int i = 0; i < totalCurrents.length; i++) {
System.out.println("totalCurrents[" + i + "] = " + totalCurrents[cIndex[i]].toString());
}
}
//
if (!bSilent) {
System.out.println("\n\n capacitive currents for each device\n");
}
Expression[] capacitiveCurrents = new Expression[fieldEdges.length];
for (int i = 0; i < fieldEdges.length; i++) {
ElectricalDevice device = (ElectricalDevice) fieldEdges[i].getData();
if (device instanceof MembraneElectricalDevice) {
MembraneElectricalDevice membraneElectricalDevice = (MembraneElectricalDevice) device;
Expression source = new Expression(membraneElectricalDevice.getSourceSymbol(), fieldMathMapping.getNameScope());
capacitiveCurrents[i] = Expression.add(totalCurrents[i], Expression.negate(source));
capacitiveCurrents[i].bindExpression(membraneElectricalDevice);
// membraneElectricalDevice.setCapacitiveCurrentExpression(capacitiveCurrents[i]);
} else {
// device.setCapacitiveCurrentExpression(new Expression(0.0));
}
}
if (!bSilent) {
for (int i = 0; i < capacitiveCurrents.length; i++) {
System.out.println("capacitiveCurrents[" + i + "] = " + ((capacitiveCurrents[cIndex[i]] == null) ? "0.0" : capacitiveCurrents[cIndex[i]].infix()));
}
}
//
// display equations for independent voltages.
//
// if (!bSilent){
// for (int i = 0; i < graph.getNumEdges(); i++){
// ElectricalDevice device = (ElectricalDevice)graph.getEdges()[i].getData();
// //
// // membrane ode
// //
// if (device.hasCapacitance() && dependentVoltageExpressions[i]==null){
// Expression initExp = new Expression(0.0);
// System.out.println(device.getInitialVoltageFunction().getVCML());
// System.out.println((new cbit.vcell.math.OdeEquation(new cbit.vcell.math.VolVariable(device.getVName()),new Expression(device.getInitialVoltageFunction().getName()),new Expression(fieldCapacitiveCurrent[i].flatten().toString()+"*"+cbit.vcell.model.ReservedSymbol.KMILLIVOLTS.getName()+"/"+device.getCapName()))).getVCML());
// }
// //
// // membrane forced potential
// //
// if (device.hasCapacitance() && dependentVoltageExpressions[i]!=null){
// System.out.println((new Function(device.getVName(),dependentVoltageExpressions[i])).getVCML());
// System.out.println((new Function(device.getSourceName(),device.getCurrentSourceExpression()).getVCML()));
// }
// //
// // current clamp
// //
// if (!device.hasCapacitance() && !device.isVoltageSource()){
// System.out.println((new Function(device.getSourceName(),device.getCurrentSourceExpression()).getVCML()));
// System.out.println((new Function(device.getVName(),dependentVoltageExpressions[i])).getVCML());
// }
// //
// // voltage clamp
// //
// if (!device.hasCapacitance() && device.isVoltageSource()){
// System.out.println((new Function(device.getIName(),totalCurrents[i])).getVCML());
// System.out.println((new Function(device.getVName(),device.getInitialVoltageFunction().getExpression())).getVCML());
// }
// }
// }
}
Also used :
Path(cbit.util.graph.Path)
Node(cbit.util.graph.Node)
RationalExp(cbit.vcell.matrix.RationalExp)
MappingException(cbit.vcell.mapping.MappingException)
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
SymbolTableEntry(cbit.vcell.parser.SymbolTableEntry)
Graph(cbit.util.graph.Graph)
RationalExpMatrix(cbit.vcell.matrix.RationalExpMatrix)
Expression(cbit.vcell.parser.Expression)
Edge(cbit.util.graph.Edge)
ModelUnitSystem(cbit.vcell.model.ModelUnitSystem)
Example 37 with VCUnitDefinition
use of cbit.vcell.units.VCUnitDefinition in project vcell by virtualcell.
the class MathMapping_4_8 method refreshKFluxParameters.
/**
* Insert the method's description here.
* Creation date: (6/13/2004 8:10:34 AM)
*/
private void refreshKFluxParameters() throws ExpressionException {
MathMapping_4_8.MathMappingParameter[] newMathMappingParameters = (MathMappingParameter[]) fieldMathMappingParameters.clone();
//
for (int i = 0; i < newMathMappingParameters.length; i++) {
if (newMathMappingParameters[i].getRole() == PARAMETER_ROLE_KFLUX) {
newMathMappingParameters = (MathMappingParameter[]) BeanUtils.removeElement(newMathMappingParameters, newMathMappingParameters[i]);
}
}
//
// Add new KFlux Parameters
//
StructureMapping[] structureMappings = simContext.getGeometryContext().getStructureMappings();
StructureTopology structTopology = simContext.getModel().getStructureTopology();
VCUnitDefinition lengthInverseUnit = simContext.getModel().getX().getUnitDefinition().getInverse();
for (int i = 0; i < structureMappings.length; i++) {
if (structureMappings[i] instanceof MembraneMapping && !getResolved(structureMappings[i])) {
MembraneMapping membraneMapping = (MembraneMapping) structureMappings[i];
//
// add "inside" flux correction
//
Expression insideCorrectionExp = getInsideFluxCorrectionExpression(membraneMapping);
insideCorrectionExp.bindExpression(this);
Membrane memMappingMembrane = membraneMapping.getMembrane();
Feature insideFeature = structTopology.getInsideFeature(memMappingMembrane);
String membraneNameWithScope = membraneMapping.getNameScope().getName();
String insideName = "KFlux_" + membraneNameWithScope + "_" + insideFeature.getNameScope().getName();
KFluxParameter insideKFluxParameter = new KFluxParameter(insideName, insideCorrectionExp, lengthInverseUnit, membraneMapping, insideFeature);
newMathMappingParameters = (MathMappingParameter[]) BeanUtils.addElement(newMathMappingParameters, insideKFluxParameter);
//
// add "outside" flux correction
//
Expression outsideCorrectionExp = getOutsideFluxCorrectionExpression(simContext, membraneMapping);
outsideCorrectionExp.bindExpression(this);
Feature outsideFeature = structTopology.getOutsideFeature(memMappingMembrane);
String outsideName = "KFlux_" + membraneNameWithScope + "_" + outsideFeature.getNameScope().getName();
KFluxParameter outsideKFluxParameter = new KFluxParameter(outsideName, outsideCorrectionExp, lengthInverseUnit, membraneMapping, outsideFeature);
newMathMappingParameters = (MathMappingParameter[]) BeanUtils.addElement(newMathMappingParameters, outsideKFluxParameter);
}
}
try {
setMathMapppingParameters(newMathMappingParameters);
} catch (java.beans.PropertyVetoException e) {
e.printStackTrace(System.out);
throw new RuntimeException(e.getMessage());
}
}
Also used :
MembraneMapping(cbit.vcell.mapping.MembraneMapping)
StructureTopology(cbit.vcell.model.Model.StructureTopology)
StructureMapping(cbit.vcell.mapping.StructureMapping)
Feature(cbit.vcell.model.Feature)
PropertyVetoException(java.beans.PropertyVetoException)
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
Expression(cbit.vcell.parser.Expression)
Membrane(cbit.vcell.model.Membrane)
Example 38 with VCUnitDefinition
use of cbit.vcell.units.VCUnitDefinition 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)
Example 39 with VCUnitDefinition
use of cbit.vcell.units.VCUnitDefinition in project vcell by virtualcell.
the class AbstractMathMapping method getFluxCorrectionParameter.
/**
* This method was created in VisualAge.
* @return Expression
*/
public KFluxParameter getFluxCorrectionParameter(StructureMapping sourceStructureMapping, StructureMapping targetStructureMapping) throws MappingException, ExpressionException {
for (int i = 0; i < fieldMathMappingParameters.length; i++) {
if (fieldMathMappingParameters[i] instanceof KFluxParameter) {
KFluxParameter kfluxParameter = (KFluxParameter) fieldMathMappingParameters[i];
if (kfluxParameter.getSourceStructureMapping() == sourceStructureMapping && kfluxParameter.getTargetStructureMapping() == targetStructureMapping) {
return kfluxParameter;
}
}
}
//
// not found, add new parameter
//
String sourceName = sourceStructureMapping.getStructure().getNameScope().getName();
String targetName = targetStructureMapping.getStructure().getNameScope().getName();
Parameter sourceSizeParameter = null;
Parameter targetSizeParameter = null;
if (sourceStructureMapping.getGeometryClass() instanceof CompartmentSubVolume) {
sourceSizeParameter = sourceStructureMapping.getSizeParameter();
} else {
sourceSizeParameter = sourceStructureMapping.getUnitSizeParameter();
}
if (targetStructureMapping.getGeometryClass() instanceof CompartmentSubVolume) {
targetSizeParameter = targetStructureMapping.getSizeParameter();
if (targetSizeParameter == null || targetSizeParameter.getExpression() == null) {
throw new MappingException("structure mapping sizes not set for application " + simContext.getName());
}
} else {
targetSizeParameter = targetStructureMapping.getUnitSizeParameter();
}
Expression fluxCorrectionExp = Expression.div(new Expression(sourceSizeParameter, simContext.getNameScope()), new Expression(targetSizeParameter, simContext.getNameScope()));
VCUnitDefinition sourceSizeUnit = sourceSizeParameter.getUnitDefinition();
VCUnitDefinition targetSizeUnit = targetSizeParameter.getUnitDefinition();
VCUnitDefinition unit = sourceSizeUnit.divideBy(targetSizeUnit);
fluxCorrectionExp.bindExpression(this);
String parameterName = PARAMETER_K_FLUX_PREFIX + sourceName + "_" + targetName;
KFluxParameter kFluxParameter = new KFluxParameter(parameterName, fluxCorrectionExp, unit, sourceStructureMapping, targetStructureMapping);
MathMappingParameter[] newMathMappingParameters = (MathMappingParameter[]) BeanUtils.addElement(this.fieldMathMappingParameters, kFluxParameter);
try {
setMathMapppingParameters(newMathMappingParameters);
} catch (java.beans.PropertyVetoException e) {
e.printStackTrace(System.out);
throw new RuntimeException(e.getMessage());
}
return kFluxParameter;
}
Also used :
PropertyVetoException(java.beans.PropertyVetoException)
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
CompartmentSubVolume(cbit.vcell.geometry.CompartmentSubVolume)
Expression(cbit.vcell.parser.Expression)
ProxyParameter(cbit.vcell.model.ProxyParameter)
Parameter(cbit.vcell.model.Parameter)
KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter)
LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter)
ModelParameter(cbit.vcell.model.Model.ModelParameter)
Example 40 with VCUnitDefinition
use of cbit.vcell.units.VCUnitDefinition in project vcell by virtualcell.
the class AbstractStochMathMapping method getExpressionConcToExpectedCount.
/**
* getExpressionConcToAmt : converts the concentration expression ('concExpr') to an expression of the number of particles.
* If argument 'speciesContext' is on a membrane, particlesExpr = concExpr * size_of_Mem. If 'speciesContext' is in
* feature, particlesExpr = (concExpr * size_of_Feature)/KMOLE.
* @param concExpr
* @param speciesContext
* @return
* @throws MappingException
* @throws ExpressionException
*/
protected Expression getExpressionConcToExpectedCount(Expression concExpr, Structure structure) throws MappingException, ExpressionException {
Expression exp = Expression.mult(concExpr, new Expression(structure.getStructureSize(), getNameScope()));
ModelUnitSystem unitSystem = getSimulationContext().getModel().getUnitSystem();
VCUnitDefinition substanceUnit = unitSystem.getSubstanceUnit(structure);
Expression unitFactor = getUnitFactor(unitSystem.getStochasticSubstanceUnit().divideBy(substanceUnit));
Expression particlesExpr = Expression.mult(exp, unitFactor);
return particlesExpr;
}
Also used :
VCUnitDefinition(cbit.vcell.units.VCUnitDefinition)
Expression(cbit.vcell.parser.Expression)
ModelUnitSystem(cbit.vcell.model.ModelUnitSystem)
Aggregations
VCUnitDefinition (cbit.vcell.units.VCUnitDefinition)113
Expression (cbit.vcell.parser.Expression)73
ModelUnitSystem (cbit.vcell.model.ModelUnitSystem)36
ExpressionException (cbit.vcell.parser.ExpressionException)26
PropertyVetoException (java.beans.PropertyVetoException)24
LocalParameter (cbit.vcell.mapping.ParameterContext.LocalParameter)21
Element (org.jdom.Element)20
KineticsParameter (cbit.vcell.model.Kinetics.KineticsParameter)19
ArrayList (java.util.ArrayList)17
ModelParameter (cbit.vcell.model.Model.ModelParameter)16
Model (cbit.vcell.model.Model)14
Parameter (cbit.vcell.model.Parameter)14
SpeciesContext (cbit.vcell.model.SpeciesContext)13
SymbolTableEntry (cbit.vcell.parser.SymbolTableEntry)13
Membrane (cbit.vcell.model.Membrane)12
Structure (cbit.vcell.model.Structure)12
StructureMapping (cbit.vcell.mapping.StructureMapping)11
MathException (cbit.vcell.math.MathException)10
ReactionStep (cbit.vcell.model.ReactionStep)10
Feature (cbit.vcell.model.Feature)9
