Example 1 with Path
use of cbit.util.graph.Path 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 2 with Path
use of cbit.util.graph.Path in project vcell by virtualcell.
the class StructureSorter method sortStructures.
public static Structure[] sortStructures(Model model) {
Graph graph = new Graph();
// add all structures to the graph
for (Structure structure : model.getStructures()) {
graph.addNode(new Node(structure.getName(), structure));
}
// add an edge
for (Structure structure : model.getStructures()) {
Structure parentStructure = null;
if (structure.getSbmlParentStructure() != null) {
parentStructure = structure.getSbmlParentStructure();
} else if (model.getStructureTopology().getParentStructure(structure) != null) {
parentStructure = model.getStructureTopology().getParentStructure(structure);
}
if (parentStructure != null) {
int index1 = graph.getIndex(graph.getNode(structure.getName()));
int index2 = graph.getIndex(graph.getNode(parentStructure.getName()));
StructureRelationshipEdge edge = getStructureRelationshipEdge(graph, index1, index2);
edge.bParentChild = true;
}
}
for (ReactionStep rs : model.getReactionSteps()) {
int index1 = graph.getIndex(graph.getNode(rs.getStructure().getName()));
for (ReactionParticipant rp : rs.getReactionParticipants()) {
if (rp.getStructure() != rs.getStructure()) {
int index2 = graph.getIndex(graph.getNode(rp.getStructure().getName()));
StructureRelationshipEdge edge = getStructureRelationshipEdge(graph, index2, index1);
edge.numConnections++;
}
}
}
// remove edges such that all nodes have degree <= 2
for (Node node : graph.getNodes()) {
Edge[] adjacentEdges = graph.getAdjacentEdges(node);
while (adjacentEdges.length > 2) {
StructureRelationshipEdge lowestCostEdge = getLowestCostEdge(Arrays.asList(adjacentEdges));
graph.remove(lowestCostEdge);
adjacentEdges = graph.getAdjacentEdges(node);
}
}
// for each dependency loop, remove edge with lowest cost
Path[] fundamentalCycles = graph.getFundamentalCycles();
while (fundamentalCycles.length > 0) {
StructureRelationshipEdge lowestCostEdge = getLowestCostEdge(Arrays.asList(fundamentalCycles[0].getEdges()));
graph.remove(lowestCostEdge);
fundamentalCycles = graph.getFundamentalCycles();
}
// all graphs in the forest are linear now ... find path from first to last and line up all trees in list.
Tree[] spanningForest = graph.getSpanningForest();
ArrayList<Structure> structures = new ArrayList<Structure>();
for (Tree tree : spanningForest) {
if (tree.getNodes().length > 1) {
Node start = null;
Node end = null;
for (Node node : tree.getNodes()) {
if (tree.getDegree(node) == 1) {
if (start == null) {
start = node;
} else if (end == null) {
end = node;
} else {
break;
}
}
}
Path path = tree.getTreePath(start, end);
for (Node node : path.getNodesTraversed()) {
structures.add((Structure) node.getData());
}
} else {
structures.add((Structure) tree.getNodes()[0].getData());
}
}
return structures.toArray(new Structure[structures.size()]);
}Aggregations
Edge (cbit.util.graph.Edge)2
Graph (cbit.util.graph.Graph)2
Node (cbit.util.graph.Node)2
Path (cbit.util.graph.Path)2
Tree (cbit.util.graph.Tree)1
MappingException (cbit.vcell.mapping.MappingException)1
RationalExp (cbit.vcell.matrix.RationalExp)1
RationalExpMatrix (cbit.vcell.matrix.RationalExpMatrix)1
ModelUnitSystem (cbit.vcell.model.ModelUnitSystem)1
Expression (cbit.vcell.parser.Expression)1
SymbolTableEntry (cbit.vcell.parser.SymbolTableEntry)1
VCUnitDefinition (cbit.vcell.units.VCUnitDefinition)1
ArrayList (java.util.ArrayList)1
