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

Example 1 with RationalExp

use of cbit.vcell.matrix.RationalExp 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) {
            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 2 with RationalExp

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

the class RationalExpUtils method getRationalExp.

private static RationalExp getRationalExp(SimpleNode node, boolean bAllowFunctions) throws ExpressionException {
    if (node == null) {
        return null;
    }
    if (node instanceof ASTAndNode || node instanceof ASTOrNode || node instanceof ASTNotNode || node instanceof ASTRelationalNode || node instanceof DerivativeNode || node instanceof ASTLaplacianNode) {
        throw new ExpressionException("sub-expression " + node.infixString(SimpleNode.LANGUAGE_DEFAULT) + " cannot be translated to a Rational Expression");
    } else if (node instanceof ASTFuncNode) {
        if (((ASTFuncNode) node).getFunction() == FunctionType.POW) {
            try {
                double constantExponent = node.jjtGetChild(1).evaluateConstant();
                if (constantExponent != Math.floor(constantExponent)) {
                    throw new ExpressionException("sub-expression " + node.infixString(SimpleNode.LANGUAGE_DEFAULT) + " cannot be translated to a Rational Expression, exponent is not an integer");
                }
                int intExponent = (int) constantExponent;
                if (intExponent == 0) {
                    return RationalExp.ONE;
                }
                RationalExp base = getRationalExp((SimpleNode) node.jjtGetChild(0), bAllowFunctions);
                if (intExponent < 0) {
                    base = base.inverse();
                    intExponent = -intExponent;
                }
                RationalExp baseUnit = new RationalExp(base);
                for (int i = 1; i < intExponent; i++) {
                    base = base.mult(baseUnit);
                }
                return base;
            } catch (ExpressionException e) {
                throw new ExpressionException("sub-expression " + node.infixString(SimpleNode.LANGUAGE_DEFAULT) + " cannot be translated to a Rational Expression, exponent is not constant");
            }
        } else {
            if (bAllowFunctions) {
                return new RationalExp(((ASTFuncNode) node).infixString(SimpleNode.LANGUAGE_DEFAULT));
            } else {
                throw new ExpressionException("sub-expression " + node.infixString(SimpleNode.LANGUAGE_DEFAULT) + " cannot be translated to a Rational Expression");
            }
        }
    } else if (node instanceof ASTPowerNode) {
        try {
            double constantExponent = node.jjtGetChild(1).evaluateConstant();
            if (constantExponent != Math.floor(constantExponent)) {
                throw new ExpressionException("sub-expression " + node.infixString(SimpleNode.LANGUAGE_DEFAULT) + " cannot be translated to a Rational Expression, exponent is not an integer");
            }
            int intExponent = (int) constantExponent;
            if (intExponent == 0) {
                return RationalExp.ONE;
            }
            RationalExp base = getRationalExp((SimpleNode) node.jjtGetChild(0), bAllowFunctions);
            if (intExponent < 0) {
                base = base.inverse();
                intExponent = -intExponent;
            }
            RationalExp baseUnit = new RationalExp(base);
            for (int i = 1; i < intExponent; i++) {
                base = base.mult(baseUnit);
            }
            return base;
        } catch (ExpressionException e) {
            throw new ExpressionException("sub-expression " + node.infixString(SimpleNode.LANGUAGE_DEFAULT) + " cannot be translated to a Rational Expression, exponent is not constant");
        }
    } else if (node instanceof ASTAddNode) {
        RationalExp exp = RationalExp.ZERO;
        for (int i = 0; i < node.jjtGetNumChildren(); i++) {
            exp = exp.add(getRationalExp((SimpleNode) node.jjtGetChild(i), bAllowFunctions));
        }
        return exp;
    } else if (node instanceof ASTMinusTermNode) {
        return getRationalExp((SimpleNode) node.jjtGetChild(0), bAllowFunctions).minus();
    } else if (node instanceof ASTMultNode) {
        if (node.jjtGetNumChildren() == 1) {
            return getRationalExp((SimpleNode) node.jjtGetChild(0), bAllowFunctions);
        }
        RationalExp exp = RationalExp.ONE;
        for (int i = 0; i < node.jjtGetNumChildren(); i++) {
            exp = exp.mult(getRationalExp((SimpleNode) node.jjtGetChild(i), bAllowFunctions));
        }
        return exp;
    } else if (node instanceof ASTInvertTermNode) {
        SimpleNode child = (SimpleNode) node.jjtGetChild(0);
        return getRationalExp(child, bAllowFunctions).inverse();
    } else if (node instanceof ASTFloatNode) {
        // return TBD instead of dimensionless.
        RationalNumber r = RationalNumber.getApproximateFraction(((ASTFloatNode) node).value.doubleValue());
        return new RationalExp(r);
    // } else if (node instanceof ASTFuncNode) {
    // String functionName = ((ASTFuncNode)node).getName();
    // if (functionName.equalsIgnoreCase("pow")) {
    // SimpleNode child0 =  (SimpleNode)node.jjtGetChild(0);
    // SimpleNode child1 =  (SimpleNode)node.jjtGetChild(1);
    // VCUnitDefinition unit0 = getRationalExp(child0);
    // VCUnitDefinition unit1 = getRationalExp(child1);
    // if (!unit1.compareEqual(VCUnitDefinition.UNIT_DIMENSIONLESS) && !unit1.compareEqual(VCUnitDefinition.UNIT_TBD)){
    // throw new VCUnitException("exponent of '"+node.infixString(SimpleNode.LANGUAGE_DEFAULT,SimpleNode.NAMESCOPE_DEFAULT)+"' has units of "+unit0);
    // }
    // if (unit0.compareEqual(VCUnitDefinition.UNIT_DIMENSIONLESS) || unit0.isTBD()){
    // return unit0;
    // }
    // try {
    // double d = ((SimpleNode)node.jjtGetChild(1)).evaluateConstant();
    // RationalNumber rn = RationalNumber.getApproximateFraction(d);
    // return unit0.raiseTo(rn);
    // }catch(ExpressionException e){
    // return VCUnitDefinition.UNIT_TBD;  // ????? don't know the unit now
    // }
    // }
    // } else if (node instanceof ASTPowerNode) {
    // SimpleNode child0 =  (SimpleNode)node.jjtGetChild(0);
    // SimpleNode child1 =  (SimpleNode)node.jjtGetChild(1);
    // VCUnitDefinition unit0 = getRationalExp(child0);
    // VCUnitDefinition unit1 = getRationalExp(child1);
    // if (!unit1.compareEqual(VCUnitDefinition.UNIT_DIMENSIONLESS) && !unit1.compareEqual(VCUnitDefinition.UNIT_TBD)){
    // throw new VCUnitException("exponent of '"+node.infixString(SimpleNode.LANGUAGE_DEFAULT,SimpleNode.NAMESCOPE_DEFAULT)+"' has units of "+unit0);
    // }
    // if (unit0.compareEqual(VCUnitDefinition.UNIT_DIMENSIONLESS)){
    // return VCUnitDefinition.UNIT_DIMENSIONLESS;
    // }
    // boolean bConstantExponent = false;
    // double exponentValue = 1;
    // try {
    // exponentValue = child1.evaluateConstant();
    // bConstantExponent = true;
    // }catch(ExpressionException e){
    // bConstantExponent = false;
    // }
    // if (bConstantExponent){ //
    // if (unit0.isTBD()){
    // return VCUnitDefinition.UNIT_TBD;
    // }else{
    // RationalNumber rn = RationalNumber.getApproximateFraction(exponentValue);
    // return unit0.raiseTo(rn);
    // }
    // }else{
    // return VCUnitDefinition.UNIT_TBD;
    // }
    } else if (node instanceof ASTIdNode) {
        return new RationalExp(((ASTIdNode) node).name);
    } else {
        throw new ExpressionException("node type " + node.getClass().toString() + " not supported yet");
    }
}

Example 3 with RationalExp

use of cbit.vcell.matrix.RationalExp 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 4 with RationalExp

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

the class AbstractMathMapping method getUnitFactorAsRationalExp.

/**
 * This method was created in VisualAge.
 * @return Expression
 */
public RationalExp getUnitFactorAsRationalExp(VCUnitDefinition unitFactor) {
    if (unitFactor.isEquivalent(getSimulationContext().getModel().getUnitSystem().getInstance_DIMENSIONLESS())) {
        return RationalExp.ONE;
    }
    for (MathMappingParameter p : fieldMathMappingParameters) {
        if (p instanceof UnitFactorParameter && p.getUnitDefinition().isEquivalent(unitFactor)) {
            return new RationalExp(p.getName());
        }
    }
    RationalNumber factor = unitFactor.getDimensionlessScale();
    String name = PARAMETER_K_UNITFACTOR_PREFIX + TokenMangler.fixTokenStrict(unitFactor.getSymbol().replace("-", "_neg_"));
    UnitFactorParameter unitFactorParameter = new UnitFactorParameter(name, new Expression(factor), unitFactor);
    MathMappingParameter[] newMathMappingParameters = (MathMappingParameter[]) BeanUtils.addElement(this.fieldMathMappingParameters, unitFactorParameter);
    try {
        setMathMapppingParameters(newMathMappingParameters);
    } catch (java.beans.PropertyVetoException e) {
        e.printStackTrace(System.out);
        throw new RuntimeException(e.getMessage());
    }
    return new RationalExp(factor);
}

Example 5 with RationalExp

use of cbit.vcell.matrix.RationalExp 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;
        try {
            // matrix.show();
            solution = matrix.solveLinearExpressions();
            // solution[0] is forward rate.
            solution[0] = solution[0].simplify();
            // solution[1] is reverse rate.
            solution[1] = solution[1].simplify();
        } 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."));
    }
    // 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;
}