Search in sources :

Example 81 with Expression

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);
    }
}
Also used : Function(cbit.vcell.math.Function) Variable(cbit.vcell.math.Variable) VolVariable(cbit.vcell.math.VolVariable) Expression(cbit.vcell.parser.Expression) VolVariable(cbit.vcell.math.VolVariable) VariableHash(cbit.vcell.math.VariableHash) Constant(cbit.vcell.math.Constant) Domain(cbit.vcell.math.Variable.Domain)

Example 82 with Expression

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;
}
Also used : MembraneMapping(cbit.vcell.mapping.MembraneMapping) Electrode(cbit.vcell.mapping.Electrode) Node(cbit.util.graph.Node) ElectricalTopology(cbit.vcell.model.Model.ElectricalTopology) CurrentDensityClampStimulus(cbit.vcell.mapping.CurrentDensityClampStimulus) Feature(cbit.vcell.model.Feature) TotalCurrentClampStimulus(cbit.vcell.mapping.TotalCurrentClampStimulus) ElectricalStimulus(cbit.vcell.mapping.ElectricalStimulus) Graph(cbit.util.graph.Graph) Expression(cbit.vcell.parser.Expression) VoltageClampStimulus(cbit.vcell.mapping.VoltageClampStimulus) Model(cbit.vcell.model.Model) Membrane(cbit.vcell.model.Membrane) Structure(cbit.vcell.model.Structure) Edge(cbit.util.graph.Edge)

Example 83 with Expression

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());
// }
// }
// }
}
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 84 with Expression

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());
                    }
                }
            }
        }
    }
}
Also used : ReservedVariable(cbit.vcell.math.ReservedVariable) Variable(cbit.vcell.math.Variable) VolVariable(cbit.vcell.math.VolVariable) MemVariable(cbit.vcell.math.MemVariable) Expression(cbit.vcell.parser.Expression) MathException(cbit.vcell.math.MathException) FastInvariant(cbit.vcell.math.FastInvariant)

Example 85 with Expression

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;
}
Also used : MembraneMapping(cbit.vcell.mapping.MembraneMapping) StructureTopology(cbit.vcell.model.Model.StructureTopology) Expression(cbit.vcell.parser.Expression) Membrane(cbit.vcell.model.Membrane) Structure(cbit.vcell.model.Structure)

Aggregations

Expression (cbit.vcell.parser.Expression)549 ExpressionException (cbit.vcell.parser.ExpressionException)163 VCUnitDefinition (cbit.vcell.units.VCUnitDefinition)76 PropertyVetoException (java.beans.PropertyVetoException)73 Variable (cbit.vcell.math.Variable)69 ArrayList (java.util.ArrayList)58 Vector (java.util.Vector)56 MathException (cbit.vcell.math.MathException)55 VolVariable (cbit.vcell.math.VolVariable)53 SymbolTableEntry (cbit.vcell.parser.SymbolTableEntry)51 SpeciesContext (cbit.vcell.model.SpeciesContext)50 Element (org.jdom.Element)47 KineticsParameter (cbit.vcell.model.Kinetics.KineticsParameter)45 ExpressionBindingException (cbit.vcell.parser.ExpressionBindingException)45 Model (cbit.vcell.model.Model)43 Function (cbit.vcell.math.Function)42 Constant (cbit.vcell.math.Constant)41 ModelParameter (cbit.vcell.model.Model.ModelParameter)41 ModelUnitSystem (cbit.vcell.model.ModelUnitSystem)41 LocalParameter (cbit.vcell.mapping.ParameterContext.LocalParameter)38