use of cbit.vcell.parser.Expression in project vcell by virtualcell.
the class VariableHashTest method main.
/**
* Starts the application.
* @param args an array of command-line arguments
*/
public static void main(java.lang.String[] args) {
try {
VariableHash hash = new VariableHash();
hash.addVariable(new Constant("D", new Expression("4")));
hash.addVariable(new Constant("E", new Expression("D")));
hash.addVariable(new Constant("B", new Expression("C")));
hash.addVariable(new Constant("F", new Expression("A+B+C+E")));
hash.addVariable(new Constant("C", new Expression("D+5")));
hash.addVariable(new Constant("A", new Expression("B+C")));
Domain domain = null;
hash.addVariable(new VolVariable("V1", domain));
hash.addVariable(new VolVariable("V2", domain));
hash.addVariable(new Function("B1", new Expression("C1+V3"), domain));
hash.addVariable(new Function("F1", new Expression("A+B+C+E/V1+V2-C1+B1"), domain));
hash.addVariable(new Function("C1", new Expression("D+5+B1"), domain));
hash.addVariable(new VolVariable("V3", domain));
// //
// // unsorted list
// //
// Variable unsortedVars[] = hash.getTopologicallyReorderedVariables();
// for (int i = 0; i < unsortedVars.length; i++){
// System.out.println(unsortedVars[i]);
// }
//
// Topologically sorted list
//
Variable[] sortedVars = hash.getTopologicallyReorderedVariables();
for (int i = 0; i < sortedVars.length; i++) {
System.out.println(sortedVars[i]);
}
//
// Alphabetically sorted list
//
sortedVars = hash.getAlphabeticallyOrderedVariables();
for (int i = 0; i < sortedVars.length; i++) {
System.out.println(sortedVars[i]);
}
} catch (Throwable e) {
e.printStackTrace(System.out);
}
}
use of cbit.vcell.parser.Expression in project vcell by virtualcell.
the class ElectricalCircuitGraph method getCircuitGraph.
/**
* Insert the method's description here.
* Creation date: (2/19/2002 11:24:04 AM)
* @return cbit.vcell.mapping.potential.Graph
* @param simContext cbit.vcell.mapping.SimulationContext
*/
public static Graph getCircuitGraph(SimulationContext simContext, AbstractMathMapping mathMapping) throws ExpressionException {
Graph graph = new Graph();
Model model = simContext.getModel();
//
// add nodes to the graph (one for each Feature)
//
Structure[] structures = model.getStructures();
for (int i = 0; i < structures.length; i++) {
if (structures[i] instanceof Feature) {
graph.addNode(new Node(structures[i].getName(), structures[i]));
}
}
//
// add edges for all current clamp electrodes (always have dependent voltages)
//
ElectricalStimulus[] stimuli = simContext.getElectricalStimuli();
Electrode groundElectrode = simContext.getGroundElectrode();
for (int i = 0; i < stimuli.length; i++) {
ElectricalStimulus stimulus = stimuli[i];
//
// get electrodes
//
Electrode probeElectrode = stimulus.getElectrode();
if (probeElectrode == null) {
throw new RuntimeException("null electrode for electrical stimulus");
}
if (groundElectrode == null) {
throw new RuntimeException("null ground electrode for electrical stimulus");
}
// if (!membraneMapping.getResolved()){
Node groundNode = graph.getNode(groundElectrode.getFeature().getName());
Node probeNode = graph.getNode(probeElectrode.getFeature().getName());
if (stimulus instanceof CurrentDensityClampStimulus) {
CurrentDensityClampStimulus ccStimulus = (CurrentDensityClampStimulus) stimulus;
ElectricalDevice device = new CurrentClampElectricalDevice(ccStimulus, mathMapping);
Edge edge = new Edge(probeNode, groundNode, device);
graph.addEdge(edge);
} else if (stimulus instanceof TotalCurrentClampStimulus) {
TotalCurrentClampStimulus ccStimulus = (TotalCurrentClampStimulus) stimulus;
ElectricalDevice device = new CurrentClampElectricalDevice(ccStimulus, mathMapping);
Edge edge = new Edge(probeNode, groundNode, device);
graph.addEdge(edge);
}
// }
}
//
// add edges for all membranes
//
ElectricalTopology electricalTopology = simContext.getModel().getElectricalTopology();
for (int i = 0; i < structures.length; i++) {
if (structures[i] instanceof Membrane) {
Membrane membrane = (Membrane) structures[i];
MembraneMapping membraneMapping = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(membrane);
Feature positiveFeature = electricalTopology.getPositiveFeature(membrane);
Feature negativeFeature = electricalTopology.getNegativeFeature(membrane);
if (positiveFeature != null && negativeFeature != null) {
Node insideNode = graph.getNode(positiveFeature.getName());
Node outsideNode = graph.getNode(negativeFeature.getName());
//
// getTotalMembraneCurrent() already converts to "outwardCurrent" so that same convention as voltage
//
Expression currentSource = getTotalMembraneCurrent(simContext, membrane, mathMapping);
MembraneElectricalDevice device = new MembraneElectricalDevice(membraneMapping, mathMapping);
device.getParameterFromRole(ElectricalDevice.ROLE_TransmembraneCurrent).setExpression(currentSource);
Edge edge = new Edge(insideNode, outsideNode, device);
graph.addEdge(edge);
}
}
}
//
for (int i = 0; i < stimuli.length; i++) {
ElectricalStimulus stimulus = stimuli[i];
//
// get electrodes
//
Electrode probeElectrode = stimulus.getElectrode();
if (probeElectrode == null) {
throw new RuntimeException("null electrode for electrical stimulus");
}
if (groundElectrode == null) {
throw new RuntimeException("null ground electrode for electrical stimulus");
}
// if (!membraneMapping.getResolved()){
Node groundNode = graph.getNode(groundElectrode.getFeature().getName());
Node probeNode = graph.getNode(probeElectrode.getFeature().getName());
if (stimulus instanceof VoltageClampStimulus) {
VoltageClampStimulus vcStimulus = (VoltageClampStimulus) stimulus;
ElectricalDevice device = new VoltageClampElectricalDevice(vcStimulus, mathMapping);
Edge edge = new Edge(probeNode, groundNode, device);
graph.addEdge(edge);
}
// }
}
// System.out.println(graph);
return graph;
}
use of cbit.vcell.parser.Expression 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());
// }
// }
// }
}
use of cbit.vcell.parser.Expression in project vcell by virtualcell.
the class FastSystemAnalyzer method checkLinearity.
/**
* @exception java.lang.Exception The exception description.
*/
private void checkLinearity() throws MathException, ExpressionException {
Enumeration<Variable> enum1 = fastVarList.elements();
while (enum1.hasMoreElements()) {
Variable var = enum1.nextElement();
//
// d invariant
// for each variable, make sure ------------- = constant;
// d Var
//
Enumeration<FastInvariant> enum_fi = fastSystem.getFastInvariants();
while (enum_fi.hasMoreElements()) {
FastInvariant fi = enum_fi.nextElement();
Expression exp = fi.getFunction().differentiate(var.getName());
exp = MathUtilities.substituteFunctions(exp, this).flatten();
if (!exp.isNumeric()) {
// If expression is in terms of 'x','y','z' - then its ok - relax the constant requirement.
String[] symbols = exp.getSymbols();
for (int i = 0; i < symbols.length; i++) {
if (!symbols[i].equals(ReservedVariable.X.getName()) && !symbols[i].equals(ReservedVariable.Y.getName()) && !symbols[i].equals(ReservedVariable.Z.getName()) && !symbols[i].equals(ReservedVariable.TIME.getName())) {
throw new MathException("FastInvariant " + fi.getFunction().toString() + " isn't linear, d/d(" + var.getName() + ") = " + exp.toString());
}
}
}
}
}
}
use of cbit.vcell.parser.Expression in project vcell by virtualcell.
the class MathMapping_4_8 method getResidualVolumeFraction.
public Expression getResidualVolumeFraction(FeatureMapping featureMapping) throws ExpressionException {
Expression exp = new Expression(1.0);
Structure[] structures = simContext.getGeometryContext().getModel().getStructures();
StructureTopology structTopology = simContext.getModel().getStructureTopology();
for (int i = 0; i < structures.length; i++) {
//
// for each membrane that is distributed within this feature, subtract that volume fraction
// ????? beware, 1 - v1 - v2 ... can result in a negative number if we're not careful.
//
Structure struct = structures[i];
if (struct instanceof Membrane) {
if ((structTopology.getOutsideFeature((Membrane) struct)) == featureMapping.getFeature()) {
MembraneMapping mm = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(struct);
if (getResolved(mm) == false) {
exp = Expression.add(exp, Expression.negate(new Expression(mm.getVolumeFractionParameter(), simContext.getNameScope())));
}
}
}
}
return exp;
}
Aggregations