Examples with SimulationSymbolTable cbit.vcell.solver.SimulationSymbolTable used on opensource projects

Example 11 with SimulationSymbolTable

use of cbit.vcell.solver.SimulationSymbolTable in project vcell by virtualcell.

the class OdeFileWriter method createStateVariables.

private void createStateVariables() throws Exception {
    Simulation simulation = simTask.getSimulation();
    MathDescription mathDescription = simulation.getMathDescription();
    SolverTaskDescription solverTaskDescription = simulation.getSolverTaskDescription();
    // get Ode's from MathDescription and create ODEStateVariables
    Enumeration<Equation> enum1 = mathDescription.getSubDomains().nextElement().getEquations();
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    while (enum1.hasMoreElements()) {
        Equation equation = enum1.nextElement();
        if (equation instanceof OdeEquation) {
            fieldStateVariables.addElement(new ODEStateVariable((OdeEquation) equation, simSymbolTable));
        } else {
            throw new MathException("encountered non-ode equation, unsupported");
        }
    }
    // Get sensitivity variables
    Variable[] variables = simSymbolTable.getVariables();
    Vector<SensVariable> sensVariables = new Vector<SensVariable>();
    Constant sensitivityParameter = solverTaskDescription.getSensitivityParameter();
    if (sensitivityParameter != null) {
        Constant origSensParam = sensitivityParameter;
        Constant overriddenSensParam = (Constant) simSymbolTable.getVariable(origSensParam.getName());
        for (int i = 0; i < variables.length; i++) {
            if (variables[i] instanceof VolVariable) {
                VolVariable volVariable = (VolVariable) variables[i];
                SensVariable sv = new SensVariable(volVariable, overriddenSensParam);
                sensVariables.addElement(sv);
            }
        }
    }
    if (rateSensitivity == null) {
        rateSensitivity = new RateSensitivity(mathDescription, mathDescription.getSubDomains().nextElement());
    }
    if (jacobian == null) {
        jacobian = new Jacobian(mathDescription, mathDescription.getSubDomains().nextElement());
    }
    // get Jacobian and RateSensitivities from MathDescription and create SensStateVariables
    for (int v = 0; v < sensVariables.size(); v++) {
        fieldStateVariables.addElement(new SensStateVariable(sensVariables.elementAt(v), rateSensitivity, jacobian, sensVariables, simSymbolTable));
    }
}

Example 12 with SimulationSymbolTable

use of cbit.vcell.solver.SimulationSymbolTable in project vcell by virtualcell.

the class MovingBoundaryFileWriter method flattenExpression.

private Expression flattenExpression(Expression ex, VariableDomain varDomain) throws ExpressionException, MathException {
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    Variable normalX = new Variable("normalX", null) {

        public boolean compareEqual(Matchable object, boolean bIgnoreMissingDomains) {
            return false;
        }

        public String getVCML() throws MathException {
            return null;
        }
    };
    Variable normalY = new Variable("normalY", null) {

        public boolean compareEqual(Matchable object, boolean bIgnoreMissingDomains) {
            return false;
        }

        public String getVCML() throws MathException {
            return null;
        }
    };
    SymbolTable augmentedSymbolTable = new SymbolTable() {

        @Override
        public SymbolTableEntry getEntry(String identifierString) {
            if (identifierString.equals(normalX.getName())) {
                return normalX;
            }
            if (identifierString.equals(normalY.getName())) {
                return normalY;
            }
            return simSymbolTable.getEntry(identifierString);
        }

        @Override
        public void getEntries(Map<String, SymbolTableEntry> entryMap) {
            simSymbolTable.getEntries(entryMap);
            entryMap.put(normalX.getName(), normalX);
            entryMap.put(normalY.getName(), normalY);
        }
    };
    ex = new Expression(ex);
    ex.bindExpression(augmentedSymbolTable);
    Expression flattended = MathUtilities.substituteFunctions(ex, augmentedSymbolTable).flatten();
    Expression substituted = SolverUtilities.substituteSizeAndNormalFunctions(flattended, varDomain).flatten();
    return substituted;
// return simSymbolTable.substituteFunctions(ex).flatten().infix();
}

Example 13 with SimulationSymbolTable

use of cbit.vcell.solver.SimulationSymbolTable in project vcell by virtualcell.

the class NetCDFWriter method getReactionRateLaws.

/**
 * @param reacs
 * @return ReactionRateLaw[]
 */
private ReactionRateLaw[] getReactionRateLaws(Expression[] probs) throws ExpressionException, MathException {
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    ReactionRateLaw[] results = new ReactionRateLaw[probs.length];
    varInProbOrderHash = new Hashtable[probs.length];
    for (int i = 0; i < probs.length; i++) {
        varInProbOrderHash[i] = new Hashtable<String, Integer>();
        results[i] = new ReactionRateLaw();
        Expression prob = probs[i];
        String[] symbols = prob.getSymbols();
        // get stoch variables involved in the reaction
        String[] varSymbols = getVariableSymbols(symbols);
        // max allowed order for each variable(species)
        int maxOrder = 5;
        // order of the reaction
        int totalOrder = 0;
        // the expression of the rate constant
        Expression coefExp = null;
        if (symbols != null) {
            if ((varSymbols != null && !Compare.isEqual(symbols, varSymbols)) || varSymbols == null) {
                throw new MathException("Unrecognized symbols in propensity function " + prob + ". Propensity function should contain stochastic variable symbols only.");
            }
            if (symbols != null && varSymbols != null) {
                PropensitySolver.PropensityFunction pf = PropensitySolver.solvePropensity(prob, varSymbols, maxOrder);
                // get total order and save each var's order according to each reaction.
                for (int j = 0; j < pf.getSpeciesOrders().length; j++) {
                    varInProbOrderHash[i].put(pf.getSpeceisNames()[j], pf.getSpeciesOrders()[j]);
                    totalOrder = totalOrder + pf.getSpeciesOrders()[j];
                }
                // get rate constant
                coefExp = new Expression(pf.getRateExpression());
            }
        } else {
            coefExp = new Expression(prob);
        }
        if (totalOrder == 0) {
            results[i].setLawType(ReactionRateLaw.order_0);
            coefExp = new Expression(coefExp.flatten().infix() + "/6.02e23");
            try {
                coefExp.bindExpression(simSymbolTable);
                coefExp = simSymbolTable.substituteFunctions(coefExp).flatten();
                double val = coefExp.evaluateConstant();
                results[i].setRateConstant(val);
            } catch (Exception e) {
                e.printStackTrace(System.err);
                throw new ExpressionException(e.getMessage());
            }
        } else if (totalOrder == 1) {
            results[i].setLawType(ReactionRateLaw.order_1);
            coefExp = coefExp.flatten();
            try {
                coefExp.bindExpression(simSymbolTable);
                coefExp = simSymbolTable.substituteFunctions(coefExp).flatten();
                double val = coefExp.evaluateConstant();
                results[i].setRateConstant(val);
            } catch (Exception e) {
                e.printStackTrace(System.err);
                throw new ExpressionException(e.getMessage());
            }
        } else if (totalOrder == 2) {
            if (// second order, two same molecules
            varSymbols.length == 1) {
                // k= c*6.02e23/2, since in VCell, "/2" is already incorporated into c, so we don'y need to take care of this item from the conversion.
                coefExp = new Expression(coefExp.flatten().infix() + "*6.02e23");
                results[i].setLawType(ReactionRateLaw.order_2_1substrate);
                try {
                    coefExp.bindExpression(simSymbolTable);
                    coefExp = simSymbolTable.substituteFunctions(coefExp).flatten();
                    double val = coefExp.evaluateConstant();
                    results[i].setRateConstant(val);
                } catch (Exception e) {
                    e.printStackTrace(System.err);
                    throw new ExpressionException(e.getMessage());
                }
            } else if (// second order, two different molecules
            varSymbols.length == 2) {
                coefExp = new Expression(coefExp.flatten().infix() + "*6.02e23");
                results[i].setLawType(ReactionRateLaw.order_2_2substrate);
                try {
                    coefExp.bindExpression(simSymbolTable);
                    coefExp = simSymbolTable.substituteFunctions(coefExp).flatten();
                    double val = coefExp.evaluateConstant();
                    results[i].setRateConstant(val);
                } catch (Exception e) {
                    e.printStackTrace(System.err);
                    throw new ExpressionException(e.getMessage());
                }
            }
        } else if (totalOrder == 3) {
            if (// third order, three same molecules
            varSymbols.length == 1) {
                results[i].setLawType(ReactionRateLaw.order_3_1substrate);
                // use c directly.  "/3!" is already incorporated into c, so we compensate this item from the conversion.
                coefExp = new Expression(coefExp.flatten().infix() + "*6");
            } else if (// third order, three different molecules
            varSymbols.length == 3) {
                results[i].setLawType(ReactionRateLaw.order_3_3substrate);
            } else if (// third order, two different molecules
            varSymbols.length == 2) {
                results[i].setLawType(ReactionRateLaw.order_3_2substrate);
                // use c directly.  "/(2!*1!)" is already incorporated into c, so we compensate this item from the conversion.
                coefExp = new Expression(coefExp.flatten().infix() + "*2");
            }
            coefExp = coefExp.flatten();
            try {
                coefExp.bindExpression(simSymbolTable);
                coefExp = simSymbolTable.substituteFunctions(coefExp).flatten();
                double val = coefExp.evaluateConstant();
                results[i].setRateConstant(val);
            } catch (Exception e) {
                e.printStackTrace(System.err);
                throw new ExpressionException(e.getMessage());
            }
        } else // order more than 3, throw exception
        {
            throw new RuntimeException("Reaction order is more than 3, we couldn't solve it by Hybrid simulation.");
        }
    }
    return results;
}

Example 14 with SimulationSymbolTable

use of cbit.vcell.solver.SimulationSymbolTable in project vcell by virtualcell.

the class StochFileWriter method write.

/**
 * Write the model to a text file which serves as an input for Stochastic simulation engine.
 * Creation date: (6/22/2006 5:37:26 PM)
 */
public void write(String[] parameterNames) throws Exception, ExpressionException {
    Simulation simulation = simTask.getSimulation();
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    initialize();
    if (bUseMessaging) {
        writeJMSParamters();
    }
    // Write control information
    printWriter.println("<control>");
    cbit.vcell.solver.SolverTaskDescription solverTaskDescription = simulation.getSolverTaskDescription();
    cbit.vcell.solver.TimeBounds timeBounds = solverTaskDescription.getTimeBounds();
    cbit.vcell.solver.OutputTimeSpec outputTimeSpec = solverTaskDescription.getOutputTimeSpec();
    ErrorTolerance errorTolerance = solverTaskDescription.getErrorTolerance();
    NonspatialStochSimOptions stochOpt = solverTaskDescription.getStochOpt();
    printWriter.println("STARTING_TIME" + "\t" + timeBounds.getStartingTime());
    printWriter.println("ENDING_TIME " + "\t" + timeBounds.getEndingTime());
    // pw.println("MAX_ITERATION"+"\t"+outputTimeSpec.getKeepAtMost());
    printWriter.println("TOLERANCE " + "\t" + errorTolerance.getAbsoluteErrorTolerance());
    if (outputTimeSpec.isDefault()) {
        printWriter.println("SAMPLE_INTERVAL" + "\t" + ((DefaultOutputTimeSpec) outputTimeSpec).getKeepEvery());
        printWriter.println("MAX_SAVE_POINTS" + "\t" + ((DefaultOutputTimeSpec) outputTimeSpec).getKeepAtMost());
    } else if (outputTimeSpec.isUniform()) {
        printWriter.println("SAVE_PERIOD" + "\t" + ((UniformOutputTimeSpec) outputTimeSpec).getOutputTimeStep());
    }
    printWriter.println("NUM_TRIAL" + "\t" + solverTaskDescription.getStochOpt().getNumOfTrials());
    if (stochOpt.isUseCustomSeed()) {
        printWriter.println("SEED" + "\t" + stochOpt.getCustomSeed());
    } else {
        // we generate our own random seed
        RandomDataGenerator rdg = new RandomDataGenerator();
        int randomSeed = rdg.nextInt(1, Integer.MAX_VALUE);
        printWriter.println("SEED" + "\t" + randomSeed);
    }
    printWriter.println("</control>");
    printWriter.println();
    // write model information
    // Model info. will be extracted from subDomain of mathDescription
    Enumeration<SubDomain> e = simulation.getMathDescription().getSubDomains();
    SubDomain subDomain = null;
    if (e.hasMoreElements()) {
        subDomain = e.nextElement();
    }
    if (subDomain != null) {
        printWriter.println("<model>");
        // variables
        printWriter.println("<discreteVariables>");
        // Species iniCondition (if in concentration) is sampled from a poisson distribution(which has a mean of the current iniExp value)
        // There is only one subDomain for compartmental model
        List<VarIniCondition> varInis = subDomain.getVarIniConditions();
        if ((varInis != null) && (varInis.size() > 0)) {
            RandomDataGenerator dist = new RandomDataGenerator();
            if (simulation.getSolverTaskDescription().getStochOpt().isUseCustomSeed()) {
                Integer randomSeed = simulation.getSolverTaskDescription().getStochOpt().getCustomSeed();
                if (randomSeed != null) {
                    dist.reSeed(randomSeed);
                }
            }
            printWriter.println("TotalVars" + "\t" + varInis.size());
            for (VarIniCondition varIniCondition : varInis) {
                try {
                    Expression iniExp = varIniCondition.getIniVal();
                    iniExp.bindExpression(simSymbolTable);
                    iniExp = simSymbolTable.substituteFunctions(iniExp).flatten();
                    double expectedCount = iniExp.evaluateConstant();
                    // 1000 mill
                    final Integer limit = 1000000000;
                    if (limit < expectedCount) {
                        String eMessage = "The Initial count for Species '" + varIniCondition.getVar().getName() + "' is " + BigDecimal.valueOf(expectedCount).toBigInteger() + "\n";
                        eMessage += "which is higher than the internal vCell limit of " + limit + ".\n";
                        eMessage += "Please reduce the Initial Condition value for this Species or reduce the compartment size.";
                        throw new MathFormatException(eMessage);
                    }
                    long varCount = 0;
                    if (varIniCondition instanceof VarIniCount) {
                        varCount = (long) expectedCount;
                    } else {
                        if (expectedCount > 0) {
                            varCount = dist.nextPoisson(expectedCount);
                        }
                    }
                    // System.out.println("expectedCount: " + expectedCount + ", varCount: " + varCount);
                    printWriter.println(varIniCondition.getVar().getName() + "\t" + varCount);
                } catch (ExpressionException ex) {
                    ex.printStackTrace();
                    throw new MathFormatException("variable " + varIniCondition.getVar().getName() + "'s initial condition is required to be a constant.");
                }
            }
        } else
            printWriter.println("TotalVars" + "\t" + "0");
        printWriter.println("</discreteVariables>");
        printWriter.println();
        // jump processes
        printWriter.println("<jumpProcesses>");
        List<JumpProcess> jumpProcesses = subDomain.getJumpProcesses();
        if ((jumpProcesses != null) && (jumpProcesses.size() > 0)) {
            printWriter.println("TotalProcesses" + "\t" + jumpProcesses.size());
            for (int i = 0; i < jumpProcesses.size(); i++) {
                printWriter.println(jumpProcesses.get(i).getName());
            }
        } else
            printWriter.println("TotalProcesses" + "\t" + "0");
        printWriter.println("</jumpProcesses>");
        printWriter.println();
        // process description
        printWriter.println("<processDesc>");
        if ((jumpProcesses != null) && (jumpProcesses.size() > 0)) {
            printWriter.println("TotalDescriptions" + "\t" + jumpProcesses.size());
            for (int i = 0; i < jumpProcesses.size(); i++) {
                JumpProcess temProc = (JumpProcess) jumpProcesses.get(i);
                // jump process name
                printWriter.println("JumpProcess" + "\t" + temProc.getName());
                Expression probExp = temProc.getProbabilityRate();
                try {
                    probExp.bindExpression(simSymbolTable);
                    probExp = simSymbolTable.substituteFunctions(probExp).flatten();
                    if (!isValidProbabilityExpression(probExp)) {
                        throw new MathFormatException("probability rate in jump process " + temProc.getName() + " has illegal symbols(should only contain variable names).");
                    }
                } catch (cbit.vcell.parser.ExpressionException ex) {
                    ex.printStackTrace();
                    throw new cbit.vcell.parser.ExpressionException("Binding math description error in probability rate in jump process " + temProc.getName() + ". Some symbols can not be resolved.");
                }
                // Expression temp = replaceVarIniInProbability(probExp);
                // Propensity
                printWriter.println("\t" + "Propensity" + "\t" + probExp.infix());
                // effects
                printWriter.println("\t" + "Effect" + "\t" + temProc.getActions().size());
                for (int j = 0; j < temProc.getActions().size(); j++) {
                    printWriter.print("\t\t" + ((Action) temProc.getActions().get(j)).getVar().getName() + "\t" + ((Action) temProc.getActions().get(j)).getOperation());
                    printWriter.println("\t" + ((Action) temProc.getActions().get(j)).evaluateOperand());
                    printWriter.println();
                }
                // dependencies
                Vector<String> dependencies = getDependencies(temProc, jumpProcesses);
                if ((dependencies != null) && (dependencies.size() > 0)) {
                    printWriter.println("\t" + "DependentProcesses" + "\t" + dependencies.size());
                    for (int j = 0; j < dependencies.size(); j++) printWriter.println("\t\t" + dependencies.elementAt(j));
                } else
                    printWriter.println("\t" + "DependentProcesses" + "\t" + "0");
                printWriter.println();
            }
        } else
            printWriter.println("TotalDescriptions" + "\t" + "0");
        printWriter.println("</processDesc>");
        printWriter.println("</model>");
    }
// if (subDomain != null)
}

Example 15 with SimulationSymbolTable

use of cbit.vcell.solver.SimulationSymbolTable in project vcell by virtualcell.

the class AbstractSolver method getFunctionSensitivity.

/**
 * This method was created by a SmartGuide.
 * @return double[]
 * @param identifier java.lang.String
 */
public Expression getFunctionSensitivity(Expression funcExpr, Constant constant, StateVariable[] stateVariables) throws ExpressionException {
    if (stateVariables == null || stateVariables.length == 0) {
        return null;
    }
    // this uses the chain rule
    // d F(x)   del F(x)        del F(x)   del xi
    // ------ = -------- + Sum (-------- * ------)
    // d k      del k      i    del xi    del k
    // explicit     via state variables
    // dependence
    // 
    // collect the explicit term
    // 
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    Expression sensFuncExp = funcExpr.differentiate(constant.getName());
    sensFuncExp.bindExpression(simSymbolTable);
    sensFuncExp = sensFuncExp.flatten();
    for (int i = 0; i < stateVariables.length; i++) {
        if (stateVariables[i] instanceof ODEStateVariable) {
            ODEStateVariable sv = (ODEStateVariable) stateVariables[i];
            VolVariable volVar = (VolVariable) sv.getVariable();
            // 
            // get corresponding sensitivity state variable
            // 
            SensStateVariable sensStateVariable = null;
            for (int j = 0; j < stateVariables.length; j++) {
                if (stateVariables[j] instanceof SensStateVariable && ((SensVariable) ((SensStateVariable) stateVariables[j]).getVariable()).getVolVariable().getName().equals(volVar.getName()) && ((SensStateVariable) stateVariables[j]).getParameter().getName().equals(constant.getName())) {
                    sensStateVariable = (SensStateVariable) stateVariables[j];
                }
            }
            if (sensStateVariable == null) {
                System.out.println("sens of " + volVar.getName() + " wrt " + constant.getName() + " is not found");
                return null;
            }
            // 
            // get coefficient of proportionality   (e.g.  A = total + b*B + c*C ... gives dA/dK = b*dB/dK + c*dC/dK)
            // 
            Expression coeffExpression = funcExpr.differentiate(volVar.getName());
            coeffExpression.bindExpression(simSymbolTable);
            coeffExpression = MathUtilities.substituteFunctions(coeffExpression, simSymbolTable);
            coeffExpression.bindExpression(simSymbolTable);
            coeffExpression = coeffExpression.flatten();
            sensFuncExp = Expression.add(sensFuncExp, Expression.mult(coeffExpression, new Expression(sensStateVariable.getVariable().getName())));
        }
    }
    return new Expression(sensFuncExp);
}