Examples with RationalExpMatrix cbit.vcell.matrix.RationalExpMatrix used on opensource projects

Example 1 with RationalExpMatrix

use of cbit.vcell.matrix.RationalExpMatrix in project vcell by virtualcell.

the class FastSystemAnalyzer method refreshInvarianceMatrix.

/**
 */
private void refreshInvarianceMatrix() throws MathException, ExpressionException {
    // 
    // the invariance's are expressed in matrix form
    // 
    // |a a a a a -1  0  0|   |x1|   |0|
    // |a a a a a  0 -1  0| * |x2| = |0|
    // |a a a a a  0  0 -1|   |x3|   |0|
    // |x4|
    // |x5|
    // |c1|
    // |c2|
    // |c3|
    // 
    Variable[] vars = new Variable[fastVarList.size()];
    fastVarList.copyInto(vars);
    int numVars = fastVarList.size();
    int rows = fastSystem.getNumFastInvariants();
    int cols = numVars + fastSystem.getNumFastInvariants();
    RationalExpMatrix matrix = new RationalExpMatrix(rows, cols);
    Enumeration<FastInvariant> fastInvariantsEnum = fastSystem.getFastInvariants();
    for (int i = 0; i < rows && fastInvariantsEnum.hasMoreElements(); i++) {
        FastInvariant fi = (FastInvariant) fastInvariantsEnum.nextElement();
        Expression function = fi.getFunction();
        for (int j = 0; j < numVars; j++) {
            Variable var = (Variable) fastVarList.elementAt(j);
            Expression exp = function.differentiate(var.getName());
            exp.bindExpression(null);
            exp = exp.flatten();
            RationalExp coeffRationalExp = RationalExpUtils.getRationalExp(exp);
            matrix.set_elem(i, j, coeffRationalExp);
        }
        matrix.set_elem(i, numVars + i, -1);
    }
    // Print
    System.out.println("origMatrix");
    matrix.show();
    // 
    // gaussian elimination on the matrix give the following representation
    // note that some column pivoting (variable re-ordering) is sometimes required to
    // determine N-r dependent vars
    // 
    // |10i0iccc|
    // |01i0iccc|  where (c)'s are the coefficients for constants of invariances
    // |00i1iccc|        (i)'s are the coefficients for dependent vars in terms of independent vars
    // 
    // Print
    System.out.println("reducedMatrix");
    if (rows > 0) {
        try {
            matrix.gaussianElimination(new RationalExpMatrix(rows, rows));
        } catch (MatrixException e) {
            e.printStackTrace(System.out);
            throw new MathException(e.getMessage());
        }
    }
    matrix.show();
    for (int i = 0; i < vars.length; i++) {
        System.out.print(vars[i].getName() + "  ");
    }
    System.out.println("");
    // 
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < rows; j++) {
            RationalExp rexp = matrix.get(i, j).simplify();
            matrix.set_elem(i, j, rexp);
        }
    }
    for (int i = 0; i < rows; i++) {
        // 
        if (!matrix.get(i, i).isConstant() || matrix.get(i, i).getConstant().doubleValue() != 1) {
            for (int j = i + 1; j < numVars; j++) {
                if (matrix.get(i, j).isConstant() && matrix.get(i, j).getConstant().doubleValue() == 1.0) {
                    for (int ii = 0; ii < rows; ii++) {
                        RationalExp temp = matrix.get(ii, i);
                        matrix.set_elem(ii, i, matrix.get(ii, j));
                        matrix.set_elem(ii, j, temp);
                    }
                    Variable tempVar = vars[i];
                    vars[i] = vars[j];
                    vars[j] = tempVar;
                    break;
                }
            }
        }
    }
    // Print
    for (int i = 0; i < vars.length; i++) {
        System.out.print(vars[i].getName() + "  ");
    }
    System.out.println("");
    matrix.show();
    // 
    // separate into dependent and indepent variables, and chop off identity matrix (left N-r columns)
    // 
    // T       |iiccc|                   T
    // [x1 x2 x4] = -1 * |iiccc| * [x3 x5 c1 c2 c3]
    // |iiccc|
    // 
    // 
    int numInvariants = fastSystem.getNumFastInvariants();
    dependentVarList.removeAllElements();
    for (int i = 0; i < numInvariants; i++) {
        dependentVarList.addElement(vars[i]);
    }
    independentVarList.removeAllElements();
    for (int i = numInvariants; i < vars.length; i++) {
        independentVarList.addElement(vars[i]);
    }
    int new_cols = independentVarList.size() + numInvariants;
    dependencyMatrix = new RationalExpMatrix(rows, new_cols);
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < new_cols; j++) {
            RationalExp rexp = matrix.get(i, j + dependentVarList.size()).simplify().minus();
            dependencyMatrix.set_elem(i, j, rexp);
        }
    }
    // Print
    System.out.println("\n\nDEPENDENCY MATRIX");
    dependencyMatrix.show();
    System.out.print("dependent vars: ");
    for (int i = 0; i < dependentVarList.size(); i++) {
        System.out.print(((Variable) dependentVarList.elementAt(i)).getName() + "  ");
    }
    System.out.println("");
    System.out.print("independent vars: ");
    for (int i = 0; i < independentVarList.size(); i++) {
        System.out.print(((Variable) independentVarList.elementAt(i)).getName() + "  ");
    }
    System.out.println("");
}

Example 2 with RationalExpMatrix

use of cbit.vcell.matrix.RationalExpMatrix 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());
// }
// }
// }
}

Example 3 with RationalExpMatrix

use of cbit.vcell.matrix.RationalExpMatrix in project vcell by virtualcell.

the class StructureSizeSolver method updateAbsoluteStructureSizes.

public static void updateAbsoluteStructureSizes(SimulationContext simContext, Structure struct, double structSize, VCUnitDefinition structSizeUnit) throws Exception {
    StructureMapping[] structMappings = simContext.getGeometryContext().getStructureMappings();
    try {
        StructureTopology structTopology = simContext.getModel().getStructureTopology();
        SolverParameterCollection solverParameters = new SolverParameterCollection();
        ArrayList<String> unknownVars = new ArrayList<String>();
        for (StructureMapping sm : structMappings) {
            System.out.println("Structure " + sm.getStructure().getDisplayName() + " size " + sm.getStructure().getStructureSize().getExpression());
            if (sm.getStructure() instanceof Membrane) {
                MembraneMapping mm = (MembraneMapping) sm;
                StructureMappingParameter svRatioParam = mm.getSurfaceToVolumeParameter();
                StructureMappingParameter volFractParam = mm.getVolumeFractionParameter();
                StructureMappingParameter sizeParam = mm.getSizeParameter();
                solverParameters.add(new SolverParameter(svRatioParam, TokenMangler.mangleToSName("sv_" + mm.getMembrane().getName()), svRatioParam.getExpression().evaluateConstant()));
                solverParameters.add(new SolverParameter(volFractParam, TokenMangler.mangleToSName("vf_" + mm.getMembrane().getName()), volFractParam.getExpression().evaluateConstant()));
            }
            StructureMappingParameter sizeParam = sm.getSizeParameter();
            Double priorKnownValue = null;
            String varName = TokenMangler.mangleToSName("size_" + sm.getStructure().getName());
            if (sizeParam.getExpression() != null) {
                priorKnownValue = sizeParam.getExpression().evaluateConstant();
            } else if (sm.getStructure() == struct) {
                priorKnownValue = structSize;
            } else {
                unknownVars.add(varName);
            }
            solverParameters.add(new SolverParameter(sizeParam, varName, priorKnownValue));
        }
        ArrayList<Expression> expressions = new ArrayList<Expression>();
        for (int i = 0; i < structMappings.length; i++) {
            if (structMappings[i] instanceof MembraneMapping) {
                MembraneMapping membraneMapping = (MembraneMapping) structMappings[i];
                Feature insideFeature = structTopology.getInsideFeature(membraneMapping.getMembrane());
                Feature outsideFeature = structTopology.getOutsideFeature(membraneMapping.getMembrane());
                StructureMappingParameter sizeParameter = membraneMapping.getSizeParameter();
                StructureMappingParameter volFractParameter = membraneMapping.getVolumeFractionParameter();
                StructureMappingParameter surfToVolParameter = membraneMapping.getSurfaceToVolumeParameter();
                // 
                // EC eclosing cyt, which contains er and golgi
                // "(cyt_size+ er_size + golgi_size) * cyt_svRatio - PM_size"  ... implicit equation, exp == 0 is implied
                // 
                Expression sumOfInsideVolumeExp = new Expression(0.0);
                for (int j = 0; j < structMappings.length; j++) {
                    if (structMappings[j] instanceof FeatureMapping && structTopology.enclosedBy(structMappings[j].getStructure(), insideFeature)) {
                        FeatureMapping childFeatureMappingOfInside = ((FeatureMapping) structMappings[j]);
                        sumOfInsideVolumeExp = Expression.add(sumOfInsideVolumeExp, new Expression(solverParameters.getName(childFeatureMappingOfInside.getSizeParameter())));
                    }
                }
                Expression tempExpr = Expression.mult(sumOfInsideVolumeExp, new Expression(solverParameters.getName(surfToVolParameter)));
                tempExpr = Expression.add(tempExpr, new Expression("-" + solverParameters.getName(sizeParameter)));
                expressions.add(tempExpr);
                // 
                // EC eclosing cyt, which contains er and golgi
                // (EC_size + cyt_size + er_size + golgi_size) * cyt_vfRatio - (cyt_size + er_size + golgi_size)  ... implicit equation, exp == 0 is implied
                // 
                Expression sumOfParentVolumeExp = new Expression(0.0);
                for (int j = 0; j < structMappings.length; j++) {
                    if (structMappings[j] instanceof FeatureMapping && structTopology.enclosedBy(structMappings[j].getStructure(), outsideFeature)) {
                        FeatureMapping childFeatureMappingOfParent = ((FeatureMapping) structMappings[j]);
                        sumOfParentVolumeExp = Expression.add(sumOfParentVolumeExp, new Expression(TokenMangler.mangleToSName(solverParameters.getName(childFeatureMappingOfParent.getSizeParameter()))));
                    }
                }
                Expression exp = Expression.mult(sumOfParentVolumeExp, new Expression(solverParameters.getName(volFractParameter)));
                exp = Expression.add(exp, Expression.negate(sumOfInsideVolumeExp));
                expressions.add(exp);
            }
        }
        if (expressions.size() != unknownVars.size()) {
            throw new RuntimeException("number of unknowns is " + unknownVars.size() + ", number of equations is " + expressions.size());
        }
        if (unknownVars.size() == 0 && expressions.size() == 0) {
            StructureMappingParameter sizeParam = simContext.getGeometryContext().getStructureMapping(struct).getSizeParameter();
            sizeParam.setExpression(new Expression(structSize));
            return;
        }
        RationalExp[][] rowColData = new RationalExp[unknownVars.size()][unknownVars.size() + 1];
        for (int row = 0; row < unknownVars.size(); row++) {
            // 
            // verify that there is no "constant" term (without an unknown)
            // 
            // System.out.println("equation("+row+"): "+expressions.get(row).infix());
            Expression constantTerm = new Expression(expressions.get(row));
            for (String var : unknownVars) {
                constantTerm.substituteInPlace(new Expression(var), new Expression(0.0));
            }
            constantTerm = constantTerm.flatten();
            // 
            for (int col = 0; col < unknownVars.size(); col++) {
                Expression equation = new Expression(expressions.get(row));
                String colVariable = unknownVars.get(col);
                Expression deriv = equation.differentiate(colVariable).flatten();
                String[] symbols = deriv.getSymbols();
                if (symbols != null) {
                    for (String symbol : symbols) {
                        if (unknownVars.contains(symbol)) {
                            throw new RuntimeException("equation is not linear in the unknowns");
                        }
                    }
                }
                rowColData[row][col] = RationalExpUtils.getRationalExp(deriv);
            }
            rowColData[row][unknownVars.size()] = RationalExpUtils.getRationalExp(constantTerm).minus();
        }
        RationalExpMatrix rationalExpMatrix = new RationalExpMatrix(rowColData);
        // rationalExpMatrix.show();
        RationalExp[] solutions = rationalExpMatrix.solveLinearExpressions();
        double[] solutionValues = new double[solutions.length];
        for (int i = 0; i < unknownVars.size(); i++) {
            Expression vcSolution = new Expression(solutions[i].infixString());
            String[] symbols = vcSolution.getSymbols();
            if (symbols != null) {
                for (String symbol : symbols) {
                    SolverParameter p = solverParameters.get(symbol);
                    if (p.knownValue == null) {
                        throw new RuntimeException("solution for var " + unknownVars.get(i) + " is a function of unknown var " + p.name);
                    }
                    vcSolution.substituteInPlace(new Expression(symbol), new Expression(p.knownValue.doubleValue()));
                }
            }
            double value = vcSolution.flatten().evaluateConstant();
            // System.out.println(unknownVars.get(i)+" = "+value+" = "+solutions[i].infixString());
            solutionValues[i] = value;
        }
        for (int i = 0; i < unknownVars.size(); i++) {
            SolverParameter p = solverParameters.get(unknownVars.get(i));
            p.parameter.setExpression(new Expression(solutionValues[i]));
        }
        // 
        // set the one known value (if not set already by the gui).
        // 
        StructureMappingParameter sizeParam = simContext.getGeometryContext().getStructureMapping(struct).getSizeParameter();
        sizeParam.setExpression(new Expression(structSize));
        System.out.println("done");
    } catch (ExpressionException e) {
        e.printStackTrace(System.out);
        throw new Exception(e.getMessage());
    }
}

Example 4 with RationalExpMatrix

use of cbit.vcell.matrix.RationalExpMatrix in project vcell by virtualcell.

the class MassActionSolver method solveMassAction.

public static MassActionFunction solveMassAction(Parameter optionalForwardRateParameter, Parameter optionalReverseRateParameter, Expression orgExp, ReactionStep rs) throws ExpressionException, ModelException, DivideByZeroException {
    MassActionFunction maFunc = new MassActionFunction();
    // get reactants, products, overlaps, non-overlap reactants and non-overlap products
    ArrayList<ReactionParticipant> reactants = new ArrayList<ReactionParticipant>();
    ArrayList<ReactionParticipant> products = new ArrayList<ReactionParticipant>();
    ReactionParticipant[] rp = rs.getReactionParticipants();
    // should use this one to compare functional equavalent since this duplicatedExp has all params substituted.
    Expression duplicatedExp = substituteParameters(orgExp, false);
    // separate the reactants and products, fluxes, catalysts
    String rxnName = rs.getName();
    // reaction with membrane current can not be transformed to mass action
    if (rs.getPhysicsOptions() == ReactionStep.PHYSICS_MOLECULAR_AND_ELECTRICAL || rs.getPhysicsOptions() == ReactionStep.PHYSICS_ELECTRICAL_ONLY) {
        throw new ModelException("Kinetics of reaction " + rxnName + " has membrane current. It can not be automatically transformed to Mass Action kinetics.");
    }
    for (int i = 0; i < rp.length; i++) {
        if (rp[i] instanceof Reactant) {
            reactants.add(rp[i]);
        } else if (rp[i] instanceof Product) {
            products.add(rp[i]);
        } else if (rp[i] instanceof Catalyst) {
            String catalystName = rp[i].getSpeciesContext().getName();
            // only if duplictedExp is not a non-linear function of catalystName.
            if (duplicatedExp.hasSymbol(catalystName)) {
                // Only when catalyst appears in ReactionRate, we add catalyst to both reactant and product
                // the stoichiometry should be set to 1.
                ReactionParticipant catalystRP = new ReactionParticipant(null, rs, rp[i].getSpeciesContext(), 1) {

                    public boolean compareEqual(Matchable obj) {
                        ReactionParticipant rp = (ReactionParticipant) obj;
                        if (rp == null) {
                            return false;
                        }
                        if (!Compare.isEqual(getSpecies(), rp.getSpecies())) {
                            return false;
                        }
                        if (!Compare.isEqual(getStructure(), rp.getStructure())) {
                            return false;
                        }
                        if (getStoichiometry() != rp.getStoichiometry()) {
                            return false;
                        }
                        return true;
                    }

                    @Override
                    public void writeTokens(PrintWriter pw) {
                    }

                    @Override
                    public void fromTokens(CommentStringTokenizer tokens, Model model) throws Exception {
                    }
                };
                products.add(catalystRP);
                reactants.add(catalystRP);
            }
        }
    }
    /**
     * The code below is going to solve reaction with kinetics that are NOT Massaction. Or Massaction with catalysts involved.
     */
    // 
    // 2x2 rational matrix
    // 
    // lets define
    // J() is the substituted total rate expression
    // P() as the theoretical "product" term with Kf = 1
    // R() as the theoretical "reactant" term with Kr = 1
    // 
    // Kf * R([1 1 1])  - Kr * P([1 1 1])  = J([1 1 1])
    // Kf * R([2 3 4])  - Kr * P([2 3 4])  = J([2 3 4])
    // 
    // in matrix form,
    // 
    // |                                        |
    // | R([1 1 1])   -P([1 1 1])    J([1 1 1]) |
    // | R([2 3 4])   -P([2 3 4])    J([2 3 4]) |
    // |                                        |
    // 
    // 
    // example: Kf*A*B^2*C - Kr*C*A
    // J() = forwardRate*43/Kabc*A*B^2*C - (myC*5-2)*C*A
    // P() = C*A
    // R() = A*B^2*C
    // 
    // |                                              |
    // | R([1 1 1])  -P([1 1 1])    J([1 1 1])        |
    // | R([2 3 4])  -P([2 3 4])    J([2 3 4])        |
    // |                                              |
    // 
    // |                                                                              |
    // | R([1 1 1])*R([2 3 4])    -P([1 1 1])*R([2 3 4])     J([1 1 1])*R([2 3 4])    |
    // | R([2 3 4])*R([1 1 1])    -P([2 3 4])*R([1 1 1])     J([2 3 4])*R([1 1 1])    |
    // |                                                                              |
    // 
    // 
    // |                                                                                                                    |
    // | R([1 1 1])*R([2 3 4])  -P([1 1 1])*R([2 3 4])                         J([1 1 1])*R([2 3 4])                        |
    // | 0                      -P([2 3 4])*R([1 1 1])+P([1 1 1])*R([2 3 4])   J([2 3 4])*R([1 1 1])-J([1 1 1])*R([2 3 4])  |
    // |                                                                                                                    |
    // 
    // 
    // 
    // 
    Expression forwardExp = null;
    Expression reverseExp = null;
    Expression R_exp = new Expression(1);
    if (reactants.size() > 0) {
        R_exp = new Expression(1.0);
        for (ReactionParticipant reactant : reactants) {
            R_exp = Expression.mult(R_exp, Expression.power(new Expression(reactant.getName()), new Expression(reactant.getStoichiometry())));
        }
    }
    Expression P_exp = new Expression(1);
    if (products.size() > 0) {
        P_exp = new Expression(1.0);
        for (ReactionParticipant product : products) {
            P_exp = Expression.mult(P_exp, Expression.power(new Expression(product.getName()), new Expression(product.getStoichiometry())));
        }
    }
    HashSet<String> reactionParticipantNames = new HashSet<String>();
    for (ReactionParticipant reactionParticipant : rs.getReactionParticipants()) {
        reactionParticipantNames.add(reactionParticipant.getName());
    }
    Expression R_1 = new Expression(R_exp);
    Expression R_2 = new Expression(R_exp);
    Expression P_1 = new Expression(P_exp);
    Expression P_2 = new Expression(P_exp);
    Expression J_1 = new Expression(duplicatedExp);
    Expression J_2 = new Expression(duplicatedExp);
    int index = 0;
    for (String rpName : reactionParticipantNames) {
        R_1.substituteInPlace(new Expression(rpName), new Expression(1.0));
        R_2.substituteInPlace(new Expression(rpName), new Expression(primeIntNumbers[index]));
        P_1.substituteInPlace(new Expression(rpName), new Expression(1.0));
        P_2.substituteInPlace(new Expression(rpName), new Expression(primeIntNumbers[index]));
        J_1.substituteInPlace(new Expression(rpName), new Expression(1.0));
        J_2.substituteInPlace(new Expression(rpName), new Expression(primeIntNumbers[index]));
        index++;
    }
    R_1 = R_1.flatten();
    R_2 = R_2.flatten();
    P_1 = P_1.flatten();
    P_2 = P_2.flatten();
    J_1 = J_1.flatten();
    J_2 = J_2.flatten();
    // e.g. A+B <-> A+B, A+2B <-> A+2B
    if (ExpressionUtils.functionallyEquivalent(R_exp, P_exp, false, 1e-8, 1e-8)) {
        throw new ModelException(VCellErrorMessages.getMassActionSolverMessage(rs.getName(), "Identical reactants and products not supported for stochastic models."));
    }
    if (R_exp.compareEqual(new Expression(1)) && P_exp.compareEqual(new Expression(1))) {
        // no reactants, no products ... nothing to do
        forwardExp = null;
        reverseExp = null;
    } else {
        // both reactants and products
        RationalExpMatrix matrix = new RationalExpMatrix(2, 3);
        matrix.set_elem(0, 0, RationalExpUtils.getRationalExp(R_1, true));
        matrix.set_elem(0, 1, RationalExpUtils.getRationalExp(Expression.negate(P_1), true));
        matrix.set_elem(0, 2, RationalExpUtils.getRationalExp(J_1, true));
        matrix.set_elem(1, 0, RationalExpUtils.getRationalExp(R_2, true));
        matrix.set_elem(1, 1, RationalExpUtils.getRationalExp(Expression.negate(P_2), true));
        matrix.set_elem(1, 2, RationalExpUtils.getRationalExp(J_2, true));
        RationalExp[] solution = null;
        RationalExp[] originalSolution = null;
        try {
            // matrix.show();
            solution = matrix.solveLinearExpressions();
            originalSolution = matrix.solveLinearExpressions();
            // solution[0] is forward rate.
            solution[0] = solution[0].simplify();
            // solution[1] is reverse rate.
            solution[1] = solution[1].simplify();
        } catch (ArithmeticException e) {
            if (e.getMessage() == null || !e.getMessage().startsWith(jscl.math.Expression.FailedToSimplify)) {
                // if failed to simplify, continue with what we have, otherwise rethrow
                throw (e);
            }
        } catch (MatrixException e) {
            e.printStackTrace();
            throw new ModelException(VCellErrorMessages.getMassActionSolverMessage(rs.getName(), "MatrixException: " + e.getMessage()));
        }
        forwardExp = new Expression(solution[0].infixString());
        reverseExp = new Expression(solution[1].infixString());
        // for massAction, if there is no reactant in the forward rate, the forwardExp should be set to null to avoid writing jump process for the forward reaction.
        if (R_exp.compareEqual(new Expression(1)) && rs.getKinetics().getKineticsDescription().equals(KineticsDescription.MassAction)) {
            forwardExp = null;
        }
        // for massAction, if there is no product in the reverse rate, the reverseExp should be set to null to avoid writing jump process for the reverse reaction.
        if (P_exp.compareEqual(new Expression(1)) && rs.getKinetics().getKineticsDescription().equals(KineticsDescription.MassAction)) {
            reverseExp = null;
        }
    }
    if (forwardExp != null) {
        forwardExp.bindExpression(rs);
    }
    if (reverseExp != null) {
        reverseExp.bindExpression(rs);
    }
    // Reconstruct the rate based on the extracted forward rate and reverse rate. If the reconstructed rate is not equivalent to the original rate,
    // it means the original rate is not in the form of Kf*r1^n1*r2^n2-Kr*p1^m1*p2^m2.
    Expression constructedExp = reconstructedRate(forwardExp, reverseExp, reactants, products, rs.getNameScope());
    Expression orgExp_withoutCatalyst = removeCatalystFromExp(duplicatedExp, rs);
    Expression constructedExp_withoutCatalyst = removeCatalystFromExp(constructedExp, rs);
    if (!ExpressionUtils.functionallyEquivalent(orgExp_withoutCatalyst, constructedExp_withoutCatalyst, false, 1e-8, 1e-8)) {
        throw new ModelException(VCellErrorMessages.getMassActionSolverMessage(rs.getName(), "Mathmatical form incompatible with mass action."));
    } else {
        // check if forward rate constant and reverse rate constant both can be evaluated to constants(numbers) after substituting all parameters.
        if (forwardExp != null) {
            Expression forwardExpCopy = new Expression(forwardExp);
            try {
                substituteParameters(forwardExpCopy, true).evaluateConstant();
            } catch (ExpressionException e) {
                throw new ModelException(VCellErrorMessages.getMassActionSolverMessage(rs.getName(), "Problem in forward rate '" + forwardExp.infix() + "', " + e.getMessage()));
            }
            // 
            if (optionalForwardRateParameter != null) {
                Expression forwardRateParameterCopy = new Expression(optionalForwardRateParameter, optionalForwardRateParameter.getNameScope());
                try {
                    substituteParameters(forwardRateParameterCopy, true).evaluateConstant();
                    if (forwardExpCopy.compareEqual(forwardRateParameterCopy)) {
                        forwardExp = new Expression(optionalForwardRateParameter, optionalForwardRateParameter.getNameScope());
                    }
                } catch (ExpressionException e) {
                    // not expecting a problem because forwardExpCopy didn't have a problem, but in any case it is ok to swallow this exception
                    e.printStackTrace(System.out);
                }
            }
        }
        if (reverseExp != null) {
            Expression reverseExpCopy = new Expression(reverseExp);
            try {
                substituteParameters(reverseExpCopy, true).evaluateConstant();
            } catch (ExpressionException e) {
                throw new ModelException(VCellErrorMessages.getMassActionSolverMessage(rs.getName(), "Problem in reverse rate '" + reverseExp.infix() + "', " + e.getMessage()));
            }
            // 
            if (optionalReverseRateParameter != null) {
                Expression reverseRateParameterCopy = new Expression(optionalReverseRateParameter, optionalReverseRateParameter.getNameScope());
                try {
                    substituteParameters(reverseRateParameterCopy, true).evaluateConstant();
                    if (reverseExpCopy.compareEqual(reverseRateParameterCopy)) {
                        reverseExp = new Expression(optionalReverseRateParameter, optionalReverseRateParameter.getNameScope());
                    }
                } catch (ExpressionException e) {
                    // not expecting a problem because reverseExpCopy didn't have a problem, but in any case it is ok to swallow this exception
                    e.printStackTrace(System.out);
                }
            }
        }
        maFunc.setForwardRate(forwardExp);
        maFunc.setReverseRate(reverseExp);
        maFunc.setReactants(reactants);
        maFunc.setProducts(products);
        return maFunc;
    }
}

Example 5 with RationalExpMatrix

use of cbit.vcell.matrix.RationalExpMatrix in project vcell by virtualcell.

the class IDAFileWriter method writeEquations.

/**
 * Insert the method's description here.
 * Creation date: (3/8/00 10:31:52 PM)
 */
protected String writeEquations(HashMap<Discontinuity, String> discontinuityNameMap) throws MathException, ExpressionException {
    Simulation simulation = simTask.getSimulation();
    StringBuffer sb = new StringBuffer();
    MathDescription mathDescription = simulation.getMathDescription();
    if (mathDescription.hasFastSystems()) {
        // 
        // define vector of original variables
        // 
        SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
        CompartmentSubDomain subDomain = (CompartmentSubDomain) mathDescription.getSubDomains().nextElement();
        FastSystem fastSystem = subDomain.getFastSystem();
        FastSystemAnalyzer fs_Analyzer = new FastSystemAnalyzer(fastSystem, simSymbolTable);
        int numIndependent = fs_Analyzer.getNumIndependentVariables();
        int systemDim = mathDescription.getStateVariableNames().size();
        int numDependent = systemDim - numIndependent;
        // 
        // get all variables from fast system (dependent and independent)
        // 
        HashSet<String> fastSystemVarHash = new HashSet<String>();
        Enumeration<Variable> dependentfastSystemVarEnum = fs_Analyzer.getDependentVariables();
        while (dependentfastSystemVarEnum.hasMoreElements()) {
            fastSystemVarHash.add(dependentfastSystemVarEnum.nextElement().getName());
        }
        Enumeration<Variable> independentfastSystemVarEnum = fs_Analyzer.getIndependentVariables();
        while (independentfastSystemVarEnum.hasMoreElements()) {
            fastSystemVarHash.add(independentfastSystemVarEnum.nextElement().getName());
        }
        // 
        // get all equations including for variables that are not in the fastSystem (ode equations for "slow system")
        // 
        RationalExpMatrix origInitVector = new RationalExpMatrix(systemDim, 1);
        RationalExpMatrix origSlowRateVector = new RationalExpMatrix(systemDim, 1);
        RationalExpMatrix origVarColumnVector = new RationalExpMatrix(systemDim, 1);
        Enumeration<Equation> enumEquations = subDomain.getEquations();
        int varIndex = 0;
        while (enumEquations.hasMoreElements()) {
            Equation equation = enumEquations.nextElement();
            origVarColumnVector.set_elem(varIndex, 0, new RationalExp(equation.getVariable().getName()));
            Expression rateExpr = equation.getRateExpression();
            rateExpr.bindExpression(varsSymbolTable);
            rateExpr = MathUtilities.substituteFunctions(rateExpr, varsSymbolTable);
            origSlowRateVector.set_elem(varIndex, 0, new RationalExp("(" + rateExpr.flatten().infix() + ")"));
            Expression initExpr = new Expression(equation.getInitialExpression());
            initExpr.substituteInPlace(new Expression("t"), new Expression(0.0));
            initExpr = MathUtilities.substituteFunctions(initExpr, varsSymbolTable).flatten();
            origInitVector.set_elem(varIndex, 0, new RationalExp("(" + initExpr.flatten().infix() + ")"));
            varIndex++;
        }
        // 
        // make symbolic matrix for fast invariants (from FastSystem's fastInvariants as well as a new fast invariant for each variable not included in the fast system.
        // 
        RationalExpMatrix fastInvarianceMatrix = new RationalExpMatrix(numDependent, systemDim);
        int row = 0;
        for (int i = 0; i < origVarColumnVector.getNumRows(); i++) {
            // 
            if (!fastSystemVarHash.contains(origVarColumnVector.get(i, 0).infixString())) {
                fastInvarianceMatrix.set_elem(row, i, RationalExp.ONE);
                row++;
            }
        }
        Enumeration<FastInvariant> enumFastInvariants = fastSystem.getFastInvariants();
        while (enumFastInvariants.hasMoreElements()) {
            FastInvariant fastInvariant = enumFastInvariants.nextElement();
            Expression fastInvariantExpression = fastInvariant.getFunction();
            for (int col = 0; col < systemDim; col++) {
                Expression coeff = fastInvariantExpression.differentiate(origVarColumnVector.get(col, 0).infixString()).flatten();
                coeff = simSymbolTable.substituteFunctions(coeff);
                fastInvarianceMatrix.set_elem(row, col, RationalExpUtils.getRationalExp(coeff));
            }
            row++;
        }
        for (int i = 0; i < systemDim; i++) {
            sb.append("VAR " + origVarColumnVector.get(i, 0).infixString() + " INIT " + origInitVector.get(i, 0).infixString() + ";\n");
        }
        RationalExpMatrix fullMatrix = null;
        RationalExpMatrix inverseFullMatrix = null;
        RationalExpMatrix newSlowRateVector = null;
        try {
            RationalExpMatrix fastMat = ((RationalExpMatrix) fastInvarianceMatrix.transpose().findNullSpace());
            fullMatrix = new RationalExpMatrix(systemDim, systemDim);
            row = 0;
            for (int i = 0; i < fastInvarianceMatrix.getNumRows(); i++) {
                for (int col = 0; col < systemDim; col++) {
                    fullMatrix.set_elem(row, col, fastInvarianceMatrix.get(i, col));
                }
                row++;
            }
            for (int i = 0; i < fastMat.getNumRows(); i++) {
                for (int col = 0; col < systemDim; col++) {
                    fullMatrix.set_elem(row, col, fastMat.get(i, col));
                }
                row++;
            }
            inverseFullMatrix = new RationalExpMatrix(systemDim, systemDim);
            inverseFullMatrix.identity();
            RationalExpMatrix copyOfFullMatrix = new RationalExpMatrix(fullMatrix);
            copyOfFullMatrix.gaussianElimination(inverseFullMatrix);
            newSlowRateVector = new RationalExpMatrix(numDependent, 1);
            newSlowRateVector.matmul(fastInvarianceMatrix, origSlowRateVector);
        } catch (MatrixException ex) {
            ex.printStackTrace();
            throw new MathException(ex.getMessage());
        }
        sb.append("TRANSFORM\n");
        for (row = 0; row < systemDim; row++) {
            for (int col = 0; col < systemDim; col++) {
                sb.append(fullMatrix.get(row, col).getConstant().doubleValue() + " ");
            }
            sb.append("\n");
        }
        sb.append("INVERSETRANSFORM\n");
        for (row = 0; row < systemDim; row++) {
            for (int col = 0; col < systemDim; col++) {
                sb.append(inverseFullMatrix.get(row, col).getConstant().doubleValue() + " ");
            }
            sb.append("\n");
        }
        int numDifferential = numDependent;
        int numAlgebraic = numIndependent;
        sb.append("RHS DIFFERENTIAL " + numDifferential + " ALGEBRAIC " + numAlgebraic + "\n");
        int equationIndex = 0;
        while (equationIndex < numDependent) {
            // print row of mass matrix followed by slow rate corresponding to fast invariant
            Expression slowRateExp = new Expression(newSlowRateVector.get(equationIndex, 0).infixString()).flatten();
            slowRateExp.bindExpression(simSymbolTable);
            slowRateExp = MathUtilities.substituteFunctions(slowRateExp, varsSymbolTable).flatten();
            Vector<Discontinuity> v = slowRateExp.getDiscontinuities();
            for (Discontinuity od : v) {
                od = getSubsitutedAndFlattened(od, varsSymbolTable);
                String dname = discontinuityNameMap.get(od);
                if (dname == null) {
                    dname = ROOT_VARIABLE_PREFIX + discontinuityNameMap.size();
                    discontinuityNameMap.put(od, dname);
                }
                slowRateExp.substituteInPlace(od.getDiscontinuityExp(), new Expression("(" + dname + "==1)"));
            }
            sb.append(slowRateExp.infix() + ";\n");
            equationIndex++;
        }
        Enumeration<FastRate> enumFastRates = fastSystem.getFastRates();
        while (enumFastRates.hasMoreElements()) {
            // print the fastRate for this row
            Expression fastRateExp = new Expression(enumFastRates.nextElement().getFunction());
            fastRateExp = MathUtilities.substituteFunctions(fastRateExp, varsSymbolTable).flatten();
            Vector<Discontinuity> v = fastRateExp.getDiscontinuities();
            for (Discontinuity od : v) {
                od = getSubsitutedAndFlattened(od, varsSymbolTable);
                String dname = discontinuityNameMap.get(od);
                if (dname == null) {
                    dname = ROOT_VARIABLE_PREFIX + discontinuityNameMap.size();
                    discontinuityNameMap.put(od, dname);
                }
                fastRateExp.substituteInPlace(od.getDiscontinuityExp(), new Expression("(" + dname + "==1)"));
            }
            sb.append(fastRateExp.flatten().infix() + ";\n");
            equationIndex++;
        }
    } else {
        for (int i = 0; i < getStateVariableCount(); i++) {
            StateVariable stateVar = getStateVariable(i);
            Expression initExpr = new Expression(stateVar.getInitialRateExpression());
            initExpr = MathUtilities.substituteFunctions(initExpr, varsSymbolTable);
            initExpr.substituteInPlace(new Expression("t"), new Expression(0.0));
            sb.append("VAR " + stateVar.getVariable().getName() + " INIT " + initExpr.flatten().infix() + ";\n");
        }
        sb.append("TRANSFORM\n");
        for (int row = 0; row < getStateVariableCount(); row++) {
            for (int col = 0; col < getStateVariableCount(); col++) {
                sb.append((row == col ? 1 : 0) + " ");
            }
            sb.append("\n");
        }
        sb.append("INVERSETRANSFORM\n");
        for (int row = 0; row < getStateVariableCount(); row++) {
            for (int col = 0; col < getStateVariableCount(); col++) {
                sb.append((row == col ? 1 : 0) + " ");
            }
            sb.append("\n");
        }
        sb.append("RHS DIFFERENTIAL " + getStateVariableCount() + " ALGEBRAIC 0\n");
        for (int i = 0; i < getStateVariableCount(); i++) {
            StateVariable stateVar = getStateVariable(i);
            Expression rateExpr = new Expression(stateVar.getRateExpression());
            rateExpr = MathUtilities.substituteFunctions(rateExpr, varsSymbolTable).flatten();
            Vector<Discontinuity> v = rateExpr.getDiscontinuities();
            for (Discontinuity od : v) {
                od = getSubsitutedAndFlattened(od, varsSymbolTable);
                String dname = discontinuityNameMap.get(od);
                if (dname == null) {
                    dname = ROOT_VARIABLE_PREFIX + discontinuityNameMap.size();
                    discontinuityNameMap.put(od, dname);
                }
                rateExpr.substituteInPlace(od.getDiscontinuityExp(), new Expression("(" + dname + "==1)"));
            }
            sb.append(rateExpr.infix() + ";\n");
        }
    }
    return sb.toString();
}