Search in sources :

Example 1 with MatrixException

use of cbit.vcell.matrix.MatrixException 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;
}
Also used : ArrayList(java.util.ArrayList) RationalExp(cbit.vcell.matrix.RationalExp) Matchable(org.vcell.util.Matchable) ExpressionException(cbit.vcell.parser.ExpressionException) MatrixException(cbit.vcell.matrix.MatrixException) Expression(cbit.vcell.parser.Expression) RationalExpMatrix(cbit.vcell.matrix.RationalExpMatrix) CommentStringTokenizer(org.vcell.util.CommentStringTokenizer) PrintWriter(java.io.PrintWriter) HashSet(java.util.HashSet)

Example 2 with MatrixException

use of cbit.vcell.matrix.MatrixException 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("");
}
Also used : MatrixException(cbit.vcell.matrix.MatrixException) ReservedVariable(cbit.vcell.math.ReservedVariable) Variable(cbit.vcell.math.Variable) VolVariable(cbit.vcell.math.VolVariable) MemVariable(cbit.vcell.math.MemVariable) RationalExpMatrix(cbit.vcell.matrix.RationalExpMatrix) Expression(cbit.vcell.parser.Expression) MathException(cbit.vcell.math.MathException) RationalExp(cbit.vcell.matrix.RationalExp) FastInvariant(cbit.vcell.math.FastInvariant)

Example 3 with MatrixException

use of cbit.vcell.matrix.MatrixException 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();
}
Also used : Variable(cbit.vcell.math.Variable) MathDescription(cbit.vcell.math.MathDescription) MatrixException(cbit.vcell.matrix.MatrixException) RationalExpMatrix(cbit.vcell.matrix.RationalExpMatrix) HashSet(java.util.HashSet) Discontinuity(cbit.vcell.parser.Discontinuity) SimulationSymbolTable(cbit.vcell.solver.SimulationSymbolTable) Equation(cbit.vcell.math.Equation) FastRate(cbit.vcell.math.FastRate) RationalExp(cbit.vcell.matrix.RationalExp) FastInvariant(cbit.vcell.math.FastInvariant) FastSystemAnalyzer(cbit.vcell.mapping.FastSystemAnalyzer) Simulation(cbit.vcell.solver.Simulation) FastSystem(cbit.vcell.math.FastSystem) Expression(cbit.vcell.parser.Expression) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) MathException(cbit.vcell.math.MathException)

Example 4 with MatrixException

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

the class ModelOptimizationMapping method getRemappedReferenceData.

/**
 * Gets the constraintData property (cbit.vcell.opt.ConstraintData) value.
 * @return The constraintData property value.
 * @see #setConstraintData
 */
private ReferenceData getRemappedReferenceData(MathMapping mathMapping) throws MappingException {
    if (modelOptimizationSpec.getReferenceData() == null) {
        return null;
    }
    // 
    // make sure time is mapped
    // 
    ReferenceData refData = modelOptimizationSpec.getReferenceData();
    ReferenceDataMappingSpec[] refDataMappingSpecs = modelOptimizationSpec.getReferenceDataMappingSpecs();
    RowColumnResultSet rowColResultSet = new RowColumnResultSet();
    Vector<SymbolTableEntry> modelObjectList = new Vector<SymbolTableEntry>();
    Vector<double[]> dataList = new Vector<double[]>();
    // 
    // find bound columns, (time is always mapped to the first column)
    // 
    int mappedColumnCount = 0;
    for (int i = 0; i < refDataMappingSpecs.length; i++) {
        SymbolTableEntry modelObject = refDataMappingSpecs[i].getModelObject();
        if (modelObject != null) {
            int mappedColumnIndex = mappedColumnCount;
            if (modelObject instanceof Model.ReservedSymbol && ((ReservedSymbol) modelObject).isTime()) {
                mappedColumnIndex = 0;
            }
            String origRefDataColumnName = refDataMappingSpecs[i].getReferenceDataColumnName();
            int origRefDataColumnIndex = refData.findColumn(origRefDataColumnName);
            if (origRefDataColumnIndex < 0) {
                throw new RuntimeException("reference data column named '" + origRefDataColumnName + "' not found");
            }
            double[] columnData = refData.getDataByColumn(origRefDataColumnIndex);
            if (modelObjectList.contains(modelObject)) {
                throw new RuntimeException("multiple reference data columns mapped to same model object '" + modelObject.getName() + "'");
            }
            modelObjectList.insertElementAt(modelObject, mappedColumnIndex);
            dataList.insertElementAt(columnData, mappedColumnIndex);
            mappedColumnCount++;
        }
    }
    // 
    if (modelObjectList.size() == 0) {
        throw new RuntimeException("reference data was not associated with model");
    }
    if (modelObjectList.size() == 1) {
        throw new RuntimeException("reference data was not associated with model, must map time and at least one other column");
    }
    boolean bFoundTimeVar = false;
    for (SymbolTableEntry ste : modelObjectList) {
        if (ste instanceof Model.ReservedSymbol && ((ReservedSymbol) ste).isTime()) {
            bFoundTimeVar = true;
            break;
        }
    }
    if (!bFoundTimeVar) {
        throw new RuntimeException("must map time column of reference data to model");
    }
    // 
    for (int i = 0; i < modelObjectList.size(); i++) {
        SymbolTableEntry modelObject = (SymbolTableEntry) modelObjectList.elementAt(i);
        try {
            // Find by name because MathSybolMapping has different 'objects' than refDataMapping 'objects'
            Variable variable = mathMapping.getMathSymbolMapping().findVariableByName(modelObject.getName());
            if (variable != null) {
                String symbol = variable.getName();
                rowColResultSet.addDataColumn(new ODESolverResultSetColumnDescription(symbol));
            } else if (modelObject instanceof Model.ReservedSymbol && ((Model.ReservedSymbol) modelObject).isTime()) {
                Model.ReservedSymbol time = (Model.ReservedSymbol) modelObject;
                String symbol = time.getName();
                rowColResultSet.addDataColumn(new ODESolverResultSetColumnDescription(symbol));
            }
        } catch (MathException | MatrixException | ExpressionException | ModelException e) {
            e.printStackTrace();
            throw new MappingException(e.getMessage(), e);
        }
    }
    // 
    // populate data columns (time and rest)
    // 
    double[] weights = new double[rowColResultSet.getColumnDescriptionsCount()];
    weights[0] = 1.0;
    int numRows = ((double[]) dataList.elementAt(0)).length;
    int numColumns = modelObjectList.size();
    for (int j = 0; j < numRows; j++) {
        double[] row = new double[numColumns];
        for (int i = 0; i < numColumns; i++) {
            row[i] = ((double[]) dataList.elementAt(i))[j];
            if (i > 0) {
                weights[i] += row[i] * row[i];
            }
        }
        rowColResultSet.addRow(row);
    }
    for (int i = 0; i < numColumns; i++) {
        if (weights[i] == 0) {
            weights[i] = 1;
        } else {
            weights[i] = 1 / weights[i];
        }
    }
    SimpleReferenceData remappedRefData = new SimpleReferenceData(rowColResultSet, weights);
    return remappedRefData;
}
Also used : ParameterVariable(cbit.vcell.math.ParameterVariable) Variable(cbit.vcell.math.Variable) ReservedSymbol(cbit.vcell.model.Model.ReservedSymbol) ReservedSymbol(cbit.vcell.model.Model.ReservedSymbol) ExpressionException(cbit.vcell.parser.ExpressionException) MappingException(cbit.vcell.mapping.MappingException) MatrixException(cbit.vcell.matrix.MatrixException) SymbolTableEntry(cbit.vcell.parser.SymbolTableEntry) ODESolverResultSetColumnDescription(cbit.vcell.math.ODESolverResultSetColumnDescription) Vector(java.util.Vector) RowColumnResultSet(cbit.vcell.math.RowColumnResultSet) ModelException(cbit.vcell.model.ModelException) SimpleReferenceData(cbit.vcell.opt.SimpleReferenceData) SimpleReferenceData(cbit.vcell.opt.SimpleReferenceData) ReferenceData(cbit.vcell.opt.ReferenceData) MathException(cbit.vcell.math.MathException) Model(cbit.vcell.model.Model)

Example 5 with MatrixException

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

the class MathVerifier method testMathGeneration.

public MathGenerationResults testMathGeneration(KeyValue simContextKey) throws SQLException, ObjectNotFoundException, DataAccessException, XmlParseException, MappingException, MathException, MatrixException, ExpressionException, ModelException, PropertyVetoException {
    User adminUser = new User(PropertyLoader.ADMINISTRATOR_ACCOUNT, new org.vcell.util.document.KeyValue(PropertyLoader.ADMINISTRATOR_ID));
    if (lg.isTraceEnabled())
        lg.trace("Testing SimContext with key '" + simContextKey + "'");
    // get biomodel refs
    java.sql.Connection con = null;
    java.sql.Statement stmt = null;
    con = conFactory.getConnection(new Object());
    cbit.vcell.modeldb.BioModelSimContextLinkTable bmscTable = cbit.vcell.modeldb.BioModelSimContextLinkTable.table;
    cbit.vcell.modeldb.BioModelTable bmTable = cbit.vcell.modeldb.BioModelTable.table;
    cbit.vcell.modeldb.UserTable userTable = cbit.vcell.modeldb.UserTable.table;
    String sql = "SELECT " + bmscTable.bioModelRef.getQualifiedColName() + "," + bmTable.ownerRef.getQualifiedColName() + "," + userTable.userid.getQualifiedColName() + " FROM " + bmscTable.getTableName() + "," + bmTable.getTableName() + "," + userTable.getTableName() + " WHERE " + bmscTable.simContextRef.getQualifiedColName() + " = " + simContextKey + " AND " + bmTable.id.getQualifiedColName() + " = " + bmscTable.bioModelRef.getQualifiedColName() + " AND " + bmTable.ownerRef.getQualifiedColName() + " = " + userTable.id.getQualifiedColName();
    ArrayList<KeyValue> bioModelKeys = new ArrayList<KeyValue>();
    stmt = con.createStatement();
    User owner = null;
    try {
        ResultSet rset = stmt.executeQuery(sql);
        while (rset.next()) {
            KeyValue key = new KeyValue(rset.getBigDecimal(bmscTable.bioModelRef.getUnqualifiedColName()));
            bioModelKeys.add(key);
            KeyValue ownerRef = new KeyValue(rset.getBigDecimal(bmTable.ownerRef.getUnqualifiedColName()));
            String userid = rset.getString(userTable.userid.getUnqualifiedColName());
            owner = new User(userid, ownerRef);
        }
    } finally {
        if (stmt != null) {
            stmt.close();
        }
        con.close();
    }
    // use the first biomodel...
    if (bioModelKeys.size() == 0) {
        throw new RuntimeException("zombie simContext ... no biomodels");
    }
    BioModelInfo bioModelInfo = dbServerImpl.getBioModelInfo(owner, bioModelKeys.get(0));
    // 
    // read in the BioModel from the database
    // 
    BigString bioModelXML = dbServerImpl.getBioModelXML(owner, bioModelInfo.getVersion().getVersionKey());
    BioModel bioModelFromDB = XmlHelper.XMLToBioModel(new XMLSource(bioModelXML.toString()));
    BioModel bioModelNewMath = XmlHelper.XMLToBioModel(new XMLSource(bioModelXML.toString()));
    bioModelFromDB.refreshDependencies();
    bioModelNewMath.refreshDependencies();
    // 
    // get all Simulations for this model
    // 
    Simulation[] modelSimsFromDB = bioModelFromDB.getSimulations();
    // 
    // ---> only for the SimContext we started with...
    // recompute mathDescription, and verify it is equivalent
    // then check each associated simulation to ensure math overrides are applied in an equivalent manner also.
    // 
    SimulationContext[] simContextsFromDB = bioModelFromDB.getSimulationContexts();
    SimulationContext[] simContextsNewMath = bioModelNewMath.getSimulationContexts();
    SimulationContext simContextFromDB = null;
    SimulationContext simContextNewMath = null;
    for (int k = 0; k < simContextsFromDB.length; k++) {
        // find it...
        if (simContextsFromDB[k].getKey().equals(simContextKey)) {
            simContextFromDB = simContextsFromDB[k];
            simContextNewMath = simContextsNewMath[k];
            break;
        }
    }
    if (simContextFromDB == null) {
        throw new RuntimeException("BioModel referred to by this SimContext does not contain this SimContext");
    } else {
        MathDescription origMathDesc = simContextFromDB.getMathDescription();
        // 
        try {
            if (simContextNewMath.getGeometry().getDimension() > 0 && simContextNewMath.getGeometry().getGeometrySurfaceDescription().getGeometricRegions() == null) {
                simContextNewMath.getGeometry().getGeometrySurfaceDescription().updateAll();
            }
        } catch (Exception e) {
            e.printStackTrace(System.out);
        }
        // 
        // updated mathDescription loaded into copy of bioModel, then test for equivalence.
        // 
        cbit.vcell.mapping.MathMapping mathMapping = simContextNewMath.createNewMathMapping();
        MathDescription mathDesc_latest = mathMapping.getMathDescription();
        MathMapping_4_8 mathMapping_4_8 = new MathMapping_4_8(simContextNewMath);
        MathDescription mathDesc_4_8 = mathMapping_4_8.getMathDescription();
        String issueString = null;
        org.vcell.util.Issue[] issues = mathMapping.getIssues();
        if (issues != null && issues.length > 0) {
            StringBuffer buffer = new StringBuffer("Issues(" + issues.length + "):\n");
            for (int l = 0; l < issues.length; l++) {
                buffer.append(" <<" + issues[l].toString() + ">>\n");
            }
            issueString = buffer.toString();
        }
        simContextNewMath.setMathDescription(mathDesc_latest);
        MathCompareResults mathCompareResults_latest = MathDescription.testEquivalency(SimulationSymbolTable.createMathSymbolTableFactory(), origMathDesc, mathDesc_latest);
        MathCompareResults mathCompareResults_4_8 = null;
        try {
            mathCompareResults_4_8 = MathDescription.testEquivalency(SimulationSymbolTable.createMathSymbolTableFactory(), origMathDesc, mathDesc_4_8);
        } catch (Exception e) {
            e.printStackTrace(System.out);
            mathCompareResults_4_8 = new MathCompareResults(Decision.MathDifferent_FAILURE_UNKNOWN, e.getMessage());
        }
        return new MathGenerationResults(bioModelFromDB, simContextFromDB, origMathDesc, mathDesc_latest, mathCompareResults_latest, mathDesc_4_8, mathCompareResults_4_8);
    }
}
Also used : User(org.vcell.util.document.User) Connection(java.sql.Connection) KeyValue(org.vcell.util.document.KeyValue) MathDescription(cbit.vcell.math.MathDescription) ArrayList(java.util.ArrayList) BigString(org.vcell.util.BigString) BigString(org.vcell.util.BigString) MathCompareResults(cbit.vcell.math.MathCompareResults) ResultSet(java.sql.ResultSet) BioModelInfo(org.vcell.util.document.BioModelInfo) MathMapping_4_8(cbit.vcell.mapping.vcell_4_8.MathMapping_4_8) SimulationContext(cbit.vcell.mapping.SimulationContext) PropertyVetoException(java.beans.PropertyVetoException) MatrixException(cbit.vcell.matrix.MatrixException) ModelException(cbit.vcell.model.ModelException) ObjectNotFoundException(org.vcell.util.ObjectNotFoundException) SQLException(java.sql.SQLException) XmlParseException(cbit.vcell.xml.XmlParseException) DataAccessException(org.vcell.util.DataAccessException) ExpressionException(cbit.vcell.parser.ExpressionException) MappingException(cbit.vcell.mapping.MappingException) MathException(cbit.vcell.math.MathException) KeyValue(org.vcell.util.document.KeyValue) Simulation(cbit.vcell.solver.Simulation) BioModel(cbit.vcell.biomodel.BioModel) XMLSource(cbit.vcell.xml.XMLSource)

Aggregations

MatrixException (cbit.vcell.matrix.MatrixException)7 MathException (cbit.vcell.math.MathException)6 Variable (cbit.vcell.math.Variable)5 Expression (cbit.vcell.parser.Expression)5 RationalExp (cbit.vcell.matrix.RationalExp)4 RationalExpMatrix (cbit.vcell.matrix.RationalExpMatrix)4 ExpressionException (cbit.vcell.parser.ExpressionException)4 MappingException (cbit.vcell.mapping.MappingException)3 FastInvariant (cbit.vcell.math.FastInvariant)3 MathDescription (cbit.vcell.math.MathDescription)3 ModelException (cbit.vcell.model.ModelException)3 SimulationContext (cbit.vcell.mapping.SimulationContext)2 MemVariable (cbit.vcell.math.MemVariable)2 ParameterVariable (cbit.vcell.math.ParameterVariable)2 ReservedVariable (cbit.vcell.math.ReservedVariable)2 VolVariable (cbit.vcell.math.VolVariable)2 ReferenceData (cbit.vcell.opt.ReferenceData)2 SimpleReferenceData (cbit.vcell.opt.SimpleReferenceData)2 Simulation (cbit.vcell.solver.Simulation)2 ArrayList (java.util.ArrayList)2