Examples with JumpProcess cbit.vcell.math.JumpProcess used on opensource projects

Example 1 with JumpProcess

use of cbit.vcell.math.JumpProcess in project vcell by virtualcell.

the class SimulationWorkspace method getEstimatedNumTimePointsForStoch.

private static long getEstimatedNumTimePointsForStoch(SimulationSymbolTable simSymbolTable) {
    Simulation sim = simSymbolTable.getSimulation();
    SolverTaskDescription solverTaskDescription = sim.getSolverTaskDescription();
    TimeBounds tb = solverTaskDescription.getTimeBounds();
    double startTime = tb.getStartingTime();
    double endTime = tb.getEndingTime();
    OutputTimeSpec tSpec = solverTaskDescription.getOutputTimeSpec();
    // hybrid G_E and G_M are fixed time step methods using uniform output time spec
    if (tSpec.isUniform()) {
        double outputTimeStep = ((UniformOutputTimeSpec) tSpec).getOutputTimeStep();
        return (long) ((endTime - startTime) / outputTimeStep);
    }
    double maxProbability = 0;
    SubDomain subDomain = sim.getMathDescription().getSubDomains().nextElement();
    List<VarIniCondition> varInis = subDomain.getVarIniConditions();
    // get all the probability expressions
    ArrayList<Expression> probList = new ArrayList<Expression>();
    for (JumpProcess jp : subDomain.getJumpProcesses()) {
        probList.add(jp.getProbabilityRate());
    }
    // loop through probability expressions
    for (int i = 0; i < probList.size(); i++) {
        try {
            Expression pExp = new Expression(probList.get(i));
            pExp.bindExpression(simSymbolTable);
            pExp = simSymbolTable.substituteFunctions(pExp);
            pExp = pExp.flatten();
            String[] symbols = pExp.getSymbols();
            // substitute stoch vars with it's initial condition expressions
            if (symbols != null) {
                for (int j = 0; symbols != null && j < symbols.length; j++) {
                    for (int k = 0; k < varInis.size(); k++) {
                        if (symbols[j].equals(varInis.get(k).getVar().getName())) {
                            pExp.substituteInPlace(new Expression(symbols[j]), new Expression(varInis.get(k).getIniVal()));
                            break;
                        }
                    }
                }
            }
            pExp = simSymbolTable.substituteFunctions(pExp);
            pExp = pExp.flatten();
            double val = pExp.evaluateConstant();
            if (maxProbability < val) {
                maxProbability = val;
            }
        } catch (ExpressionBindingException e) {
            System.out.println("Cannot estimate the total time points for stochastic simulation!! Due to the reason below...");
            e.printStackTrace();
        } catch (ExpressionException ex) {
            System.out.println("Cannot estimate the total time points for stochastic simulation!! Due to the reason below...");
            ex.printStackTrace();
        } catch (MathException e) {
            System.out.println("Cannot estimate the total time points for stochastic simulation!! Due to the reason below...");
            e.printStackTrace();
        }
    }
    int keepEvery = 1;
    if (tSpec.isDefault()) {
        keepEvery = ((DefaultOutputTimeSpec) tSpec).getKeepEvery();
    }
    // points = (endt-startt)/(t*keepEvery) = (endt - startt)/(keepEvery*1/prob)
    long estimatedPoints = Math.round((tb.getEndingTime() - tb.getStartingTime()) * maxProbability / keepEvery) + 1;
    return estimatedPoints;
}

Example 2 with JumpProcess

use of cbit.vcell.math.JumpProcess in project vcell by virtualcell.

the class FieldUtilities method collectFieldFuncAndExpressions.

private static Hashtable<FieldFunctionArguments, Vector<Expression>> collectFieldFuncAndExpressions(MathDescription mathDesc) throws MathException, ExpressionException {
    // make sure each field only added once
    Hashtable<FieldFunctionArguments, Vector<Expression>> fieldFuncArgsExpHash = new Hashtable<FieldFunctionArguments, Vector<Expression>>();
    Enumeration<SubDomain> enum1 = mathDesc.getSubDomains();
    Enumeration<Variable> mathDescEnumeration = mathDesc.getVariables();
    while (mathDescEnumeration.hasMoreElements()) {
        Variable variable = mathDescEnumeration.nextElement();
        if (variable instanceof Function) {
            Function function = (Function) variable;
            FieldUtilities.addFieldFuncArgsAndExpToCollection(fieldFuncArgsExpHash, function.getExpression());
        }
    }
    // go through each subdomain
    while (enum1.hasMoreElements()) {
        SubDomain subDomain = enum1.nextElement();
        // go through each equation
        Enumeration<Equation> enum_equ = subDomain.getEquations();
        while (enum_equ.hasMoreElements()) {
            Equation equation = (Equation) enum_equ.nextElement();
            Vector<Expression> exs = equation.getExpressions(mathDesc);
            // go through each expresson
            for (int i = 0; i < exs.size(); i++) {
                Expression exp = (Expression) exs.get(i);
                FieldUtilities.addFieldFuncArgsAndExpToCollection(fieldFuncArgsExpHash, exp);
            }
        }
        // go through each Fast System
        FastSystem fastSystem = subDomain.getFastSystem();
        if (fastSystem != null) {
            Expression[] fsExprArr = fastSystem.getExpressions();
            for (int i = 0; i < fsExprArr.length; i++) {
                FieldUtilities.addFieldFuncArgsAndExpToCollection(fieldFuncArgsExpHash, fsExprArr[i]);
            }
        }
        // go through each Jump Process
        for (JumpProcess jumpProcess : subDomain.getJumpProcesses()) {
            Expression[] jpExprArr = jumpProcess.getExpressions();
            for (int i = 0; i < jpExprArr.length; i++) {
                FieldUtilities.addFieldFuncArgsAndExpToCollection(fieldFuncArgsExpHash, jpExprArr[i]);
            }
        }
        // go through VarInitConditions
        for (VarIniCondition varInitCond : subDomain.getVarIniConditions()) {
            Expression[] vicExprArr = new Expression[] { varInitCond.getIniVal(), varInitCond.getVar().getExpression() };
            for (int i = 0; i < vicExprArr.length; i++) {
                FieldUtilities.addFieldFuncArgsAndExpToCollection(fieldFuncArgsExpHash, vicExprArr[i]);
            }
        }
    }
    return fieldFuncArgsExpHash;
}

Example 3 with JumpProcess

use of cbit.vcell.math.JumpProcess in project vcell by virtualcell.

the class StochMathMapping method addJumpProcesses.

private void addJumpProcesses(VariableHash varHash, GeometryClass geometryClass, SubDomain subDomain) throws ExpressionException, ModelException, MappingException, MathException {
    // set up jump processes
    // get all the reactions from simulation context
    // ReactionSpec[] reactionSpecs = simContext.getReactionContext().getReactionSpecs();---need to take a look here!
    ModelUnitSystem modelUnitSystem = getSimulationContext().getModel().getUnitSystem();
    ReactionSpec[] reactionSpecs = getSimulationContext().getReactionContext().getReactionSpecs();
    for (ReactionSpec reactionSpec : reactionSpecs) {
        if (reactionSpec.isExcluded()) {
            continue;
        }
        // get the reaction
        ReactionStep reactionStep = reactionSpec.getReactionStep();
        Kinetics kinetics = reactionStep.getKinetics();
        // probability parameter from modelUnitSystem
        VCUnitDefinition probabilityParamUnit = modelUnitSystem.getStochasticSubstanceUnit().divideBy(modelUnitSystem.getTimeUnit());
        // Different ways to deal with simple reactions and flux reactions
        if (// simple reactions
        reactionStep instanceof SimpleReaction) {
            // check the reaction rate law to see if we need to decompose a reaction(reversible) into two jump processes.
            // rate constants are important in calculating the probability rate.
            // for Mass Action, we use KForward and KReverse,
            // for General Kinetics we parse reaction rate J to see if it is in Mass Action form.
            Expression forwardRate = null;
            Expression reverseRate = null;
            if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction) || kinetics.getKineticsDescription().equals(KineticsDescription.General)) {
                Expression rateExp = new Expression(kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate), reactionStep.getNameScope());
                Parameter forwardRateParameter = null;
                Parameter reverseRateParameter = null;
                if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
                    forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KForward);
                    reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KReverse);
                }
                MassActionSolver.MassActionFunction maFunc = MassActionSolver.solveMassAction(forwardRateParameter, reverseRateParameter, rateExp, reactionStep);
                if (maFunc.getForwardRate() == null && maFunc.getReverseRate() == null) {
                    throw new MappingException("Cannot generate stochastic math mapping for the reaction:" + reactionStep.getName() + "\nLooking for the rate function according to the form of k1*Reactant1^Stoir1*Reactant2^Stoir2...-k2*Product1^Stoip1*Product2^Stoip2.");
                } else {
                    if (maFunc.getForwardRate() != null) {
                        forwardRate = maFunc.getForwardRate();
                    }
                    if (maFunc.getReverseRate() != null) {
                        reverseRate = maFunc.getReverseRate();
                    }
                }
            } else // if it's macro/microscopic kinetics, we'll have them set up as reactions with only forward rate.
            if (kinetics.getKineticsDescription().equals(KineticsDescription.Macroscopic_irreversible) || kinetics.getKineticsDescription().equals(KineticsDescription.Microscopic_irreversible)) {
                Expression Kon = getIdentifierSubstitutions(new Expression(reactionStep.getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_KOn), getNameScope()), reactionStep.getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_Binding_Radius).getUnitDefinition(), geometryClass);
                if (Kon != null) {
                    Expression KonCopy = new Expression(Kon);
                    try {
                        MassActionSolver.substituteParameters(KonCopy, true).evaluateConstant();
                        forwardRate = new Expression(Kon);
                    } catch (ExpressionException e) {
                        throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Problem with Kon parameter in " + reactionStep.getName() + ":  '" + KonCopy.infix() + "', " + e.getMessage()));
                    }
                } else {
                    throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Kon parameter of " + reactionStep.getName() + " is null."));
                }
            }
            boolean isForwardRatePresent = false;
            boolean isReverseRatePresent = false;
            if (forwardRate != null) {
                isForwardRatePresent = true;
            }
            if (reverseRate != null) {
                isReverseRatePresent = true;
            }
            // we process it as forward reaction
            if ((isForwardRatePresent)) /*|| ((forwardRate == null) && (reverseRate == null))*/
            {
                // get jump process name
                String jpName = TokenMangler.mangleToSName(reactionStep.getName());
                // get probability
                Expression exp = null;
                // reactions are of mass action form
                exp = getProbabilityRate(reactionStep, forwardRate, true);
                ProbabilityParameter probParm = null;
                try {
                    probParm = addProbabilityParameter(PARAMETER_PROBABILITYRATE_PREFIX + jpName, exp, PARAMETER_ROLE_P, probabilityParamUnit, reactionStep);
                } catch (PropertyVetoException pve) {
                    pve.printStackTrace();
                    throw new MappingException(pve.getMessage());
                }
                // add probability to function or constant
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(probParm, geometryClass), getIdentifierSubstitutions(exp, probabilityParamUnit, geometryClass), geometryClass));
                JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probParm, geometryClass)));
                // actions
                ReactionParticipant[] reacPart = reactionStep.getReactionParticipants();
                for (int j = 0; j < reacPart.length; j++) {
                    Action action = null;
                    SpeciesCountParameter spCountParam = getSpeciesCountParameter(reacPart[j].getSpeciesContext());
                    if (reacPart[j] instanceof Reactant) {
                        // check if the reactant is a constant. If the species is a constant, there will be no action taken on this species
                        if (// not a constant
                        !simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
                            int stoi = ((Reactant) reacPart[j]).getStoichiometry();
                            action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(-stoi));
                            jp.addAction(action);
                        }
                    } else if (reacPart[j] instanceof Product) {
                        // check if the product is a constant. If the product is a constant, there will be no action taken on this species
                        if (// not a constant
                        !simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
                            int stoi = ((Product) reacPart[j]).getStoichiometry();
                            action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(stoi));
                            jp.addAction(action);
                        }
                    }
                }
                // add jump process to compartment subDomain
                subDomain.addJumpProcess(jp);
            }
            if (// one more jump process for a reversible reaction
            isReverseRatePresent) {
                // get jump process name
                String jpName = TokenMangler.mangleToSName(reactionStep.getName()) + PARAMETER_PROBABILITY_RATE_REVERSE_SUFFIX;
                Expression exp = null;
                // reactions are mass actions
                exp = getProbabilityRate(reactionStep, reverseRate, false);
                ProbabilityParameter probRevParm = null;
                try {
                    probRevParm = addProbabilityParameter(PARAMETER_PROBABILITYRATE_PREFIX + jpName, exp, PARAMETER_ROLE_P_reverse, probabilityParamUnit, reactionStep);
                } catch (PropertyVetoException pve) {
                    pve.printStackTrace();
                    throw new MappingException(pve.getMessage());
                }
                // add probability to function or constant
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(probRevParm, geometryClass), getIdentifierSubstitutions(exp, probabilityParamUnit, geometryClass), geometryClass));
                JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probRevParm, geometryClass)));
                // actions
                ReactionParticipant[] reacPart = reactionStep.getReactionParticipants();
                for (int j = 0; j < reacPart.length; j++) {
                    Action action = null;
                    SpeciesCountParameter spCountParam = getSpeciesCountParameter(reacPart[j].getSpeciesContext());
                    if (reacPart[j] instanceof Reactant) {
                        // check if the reactant is a constant. If the species is a constant, there will be no action taken on this species
                        if (// not a constant
                        !simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
                            int stoi = ((Reactant) reacPart[j]).getStoichiometry();
                            action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(stoi));
                            jp.addAction(action);
                        }
                    } else if (reacPart[j] instanceof Product) {
                        // check if the product is a constant. If the product is a constant, there will be no action taken on this species
                        if (// not a constant
                        !simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
                            int stoi = ((Product) reacPart[j]).getStoichiometry();
                            action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(-stoi));
                            jp.addAction(action);
                        }
                    }
                }
                // add jump process to compartment subDomain
                subDomain.addJumpProcess(jp);
            }
        // end of if(isForwardRateNonZero), if(isReverseRateNonRate)
        } else if (// flux reactions
        reactionStep instanceof FluxReaction) {
            // we could set jump processes for general flux rate in forms of p1*Sout + p2*Sin
            if (kinetics.getKineticsDescription().equals(KineticsDescription.General) || kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
                Expression fluxRate = new Expression(kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate), reactionStep.getNameScope());
                // we have to pass the math description para to flux solver, coz somehow math description in simulation context is not updated.
                // forward and reverse rate parameters may be null
                Parameter forwardRateParameter = null;
                Parameter reverseRateParameter = null;
                if (kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
                    forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
                    reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
                }
                MassActionSolver.MassActionFunction fluxFunc = MassActionSolver.solveMassAction(forwardRateParameter, reverseRateParameter, fluxRate, (FluxReaction) reactionStep);
                // create jump process for forward flux if it exists.
                Expression rsStructureSize = new Expression(reactionStep.getStructure().getStructureSize(), getNameScope());
                VCUnitDefinition probRateUnit = modelUnitSystem.getStochasticSubstanceUnit().divideBy(modelUnitSystem.getAreaUnit()).divideBy(modelUnitSystem.getTimeUnit());
                Expression rsRateUnitFactor = getUnitFactor(probRateUnit.divideBy(modelUnitSystem.getFluxReactionUnit()));
                if (fluxFunc.getForwardRate() != null && !fluxFunc.getForwardRate().isZero()) {
                    Expression rate = fluxFunc.getForwardRate();
                    // get species expression (depend on structure, if mem: Species/mem_Size, if vol: species*KMOLE/vol_size)
                    if (fluxFunc.getReactants().size() != 1) {
                        throw new MappingException("Flux " + reactionStep.getName() + " should have only one reactant.");
                    }
                    SpeciesContext scReactant = fluxFunc.getReactants().get(0).getSpeciesContext();
                    Expression scConcExpr = new Expression(getSpeciesConcentrationParameter(scReactant), getNameScope());
                    Expression probExp = Expression.mult(rate, rsRateUnitFactor, rsStructureSize, scConcExpr);
                    // jump process name
                    // +"_reverse";
                    String jpName = TokenMangler.mangleToSName(reactionStep.getName());
                    ProbabilityParameter probParm = null;
                    try {
                        probParm = addProbabilityParameter(PARAMETER_PROBABILITYRATE_PREFIX + jpName, probExp, PARAMETER_ROLE_P, probabilityParamUnit, reactionStep);
                    } catch (PropertyVetoException pve) {
                        pve.printStackTrace();
                        throw new MappingException(pve.getMessage());
                    }
                    // add probability to function or constant
                    String ms = getMathSymbol(probParm, geometryClass);
                    Expression is = getIdentifierSubstitutions(probExp, probabilityParamUnit, geometryClass);
                    Variable nfoc = newFunctionOrConstant(ms, is, geometryClass);
                    varHash.addVariable(nfoc);
                    JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probParm, geometryClass)));
                    // actions
                    Action action = null;
                    SpeciesContext sc = fluxFunc.getReactants().get(0).getSpeciesContext();
                    if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
                        action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(-1));
                        jp.addAction(action);
                    }
                    sc = fluxFunc.getProducts().get(0).getSpeciesContext();
                    if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
                        action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(1));
                        jp.addAction(action);
                    }
                    subDomain.addJumpProcess(jp);
                }
                // create jump process for reverse flux if it exists.
                if (fluxFunc.getReverseRate() != null && !fluxFunc.getReverseRate().isZero()) {
                    // jump process name
                    String jpName = TokenMangler.mangleToSName(reactionStep.getName()) + PARAMETER_PROBABILITY_RATE_REVERSE_SUFFIX;
                    Expression rate = fluxFunc.getReverseRate();
                    // get species expression (depend on structure, if mem: Species/mem_Size, if vol: species*KMOLE/vol_size)
                    if (fluxFunc.getProducts().size() != 1) {
                        throw new MappingException("Flux " + reactionStep.getName() + " should have only one product.");
                    }
                    SpeciesContext scProduct = fluxFunc.getProducts().get(0).getSpeciesContext();
                    Expression scConcExpr = new Expression(getSpeciesConcentrationParameter(scProduct), getNameScope());
                    Expression probExp = Expression.mult(rate, rsRateUnitFactor, rsStructureSize, scConcExpr);
                    ProbabilityParameter probRevParm = null;
                    try {
                        probRevParm = addProbabilityParameter(PARAMETER_PROBABILITYRATE_PREFIX + jpName, probExp, PARAMETER_ROLE_P_reverse, probabilityParamUnit, reactionStep);
                    } catch (PropertyVetoException pve) {
                        pve.printStackTrace();
                        throw new MappingException(pve.getMessage());
                    }
                    // add probability to function or constant
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(probRevParm, geometryClass), getIdentifierSubstitutions(probExp, probabilityParamUnit, geometryClass), geometryClass));
                    JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probRevParm, geometryClass)));
                    // actions
                    Action action = null;
                    SpeciesContext sc = fluxFunc.getReactants().get(0).getSpeciesContext();
                    if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
                        action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(1));
                        jp.addAction(action);
                    }
                    sc = fluxFunc.getProducts().get(0).getSpeciesContext();
                    if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
                        action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(-1));
                        jp.addAction(action);
                    }
                    subDomain.addJumpProcess(jp);
                }
            }
        }
    // end of if (simplereaction)...else if(fluxreaction)
    }
// end of reaction step loop
}

Example 4 with JumpProcess

use of cbit.vcell.math.JumpProcess 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());
        // need to overwrite limit hardcoded in C++
        double savePoints = (timeBounds.getEndingTime() - timeBounds.getStartingTime()) / ((UniformOutputTimeSpec) outputTimeSpec).getOutputTimeStep();
        printWriter.println("MAX_SAVE_POINTS" + "\t" + (Math.ceil(savePoints) + 1));
    }
    // boolean isMultiTrial = !solverTaskDescription.getStochOpt().isHistogram() &&
    // solverTaskDescription.getStochOpt().getNumOfTrials() > 1;
    // //Multi-trial 'NUM_TRIAL' handled by slurm array within .slurm.sh script
    // printWriter.println("NUM_TRIAL"+"\t"+(isMultiTrial?1:solverTaskDescription.getStochOpt().getNumOfTrials()));
    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);
    }
    if (isMultiTrialNonHisto) {
        printWriter.println("BMULTIBUTNOTHISTO" + "\t" + "1");
    }
    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) Math.round(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 5 with JumpProcess

use of cbit.vcell.math.JumpProcess in project vcell by virtualcell.

the class Xmlproducer method getXML.

/**
 * This method returns a XML representation of a CompartmentSubDomain object.
 * Creation date: (3/2/2001 1:18:55 PM)
 * @return Element
 * @param param cbit.vcell.math.CompartmentSubDomain
 */
private Element getXML(CompartmentSubDomain param) throws XmlParseException {
    Element compartment = new Element(XMLTags.CompartmentSubDomainTag);
    compartment.setAttribute(XMLTags.NameAttrTag, mangle(param.getName()));
    // Add boundaryType subelements
    Element boundary;
    // Xm
    boundary = new Element(XMLTags.BoundaryTypeTag);
    boundary.setAttribute(XMLTags.BoundaryAttrTag, XMLTags.BoundaryAttrValueXm);
    boundary.setAttribute(XMLTags.BoundaryTypeAttrTag, param.getBoundaryConditionXm().boundaryTypeStringValue());
    compartment.addContent(boundary);
    // Xp
    boundary = new Element(XMLTags.BoundaryTypeTag);
    boundary.setAttribute(XMLTags.BoundaryAttrTag, XMLTags.BoundaryAttrValueXp);
    boundary.setAttribute(XMLTags.BoundaryTypeAttrTag, param.getBoundaryConditionXp().boundaryTypeStringValue());
    compartment.addContent(boundary);
    // Ym
    boundary = new Element(XMLTags.BoundaryTypeTag);
    boundary.setAttribute(XMLTags.BoundaryAttrTag, XMLTags.BoundaryAttrValueYm);
    boundary.setAttribute(XMLTags.BoundaryTypeAttrTag, param.getBoundaryConditionYm().boundaryTypeStringValue());
    compartment.addContent(boundary);
    // Yp
    boundary = new Element(XMLTags.BoundaryTypeTag);
    boundary.setAttribute(XMLTags.BoundaryAttrTag, XMLTags.BoundaryAttrValueYp);
    boundary.setAttribute(XMLTags.BoundaryTypeAttrTag, param.getBoundaryConditionYp().boundaryTypeStringValue());
    compartment.addContent(boundary);
    // Zm
    boundary = new Element(XMLTags.BoundaryTypeTag);
    boundary.setAttribute(XMLTags.BoundaryAttrTag, XMLTags.BoundaryAttrValueZm);
    boundary.setAttribute(XMLTags.BoundaryTypeAttrTag, param.getBoundaryConditionZm().boundaryTypeStringValue());
    compartment.addContent(boundary);
    // Zp
    boundary = new Element(XMLTags.BoundaryTypeTag);
    boundary.setAttribute(XMLTags.BoundaryAttrTag, XMLTags.BoundaryAttrValueZp);
    boundary.setAttribute(XMLTags.BoundaryTypeAttrTag, param.getBoundaryConditionZp().boundaryTypeStringValue());
    compartment.addContent(boundary);
    // add BoundaryConditionSpecs
    for (BoundaryConditionSpec bcs : param.getBoundaryconditionSpecs()) {
        compartment.addContent(getXML(bcs));
    }
    // Add Equations
    Enumeration<Equation> enum1 = param.getEquations();
    while (enum1.hasMoreElements()) {
        Equation equ = enum1.nextElement();
        compartment.addContent(getXML(equ));
    }
    // Add FastSystem
    if (param.getFastSystem() != null) {
        compartment.addContent(getXML(param.getFastSystem()));
    }
    // Add Variable Initial Condition
    for (VarIniCondition varIni : param.getVarIniConditions()) {
        compartment.addContent(getXML(varIni));
    }
    // Add JumpProcesses
    for (JumpProcess jp : param.getJumpProcesses()) {
        compartment.addContent(getXML(jp));
    }
    for (ParticleJumpProcess pjp : param.getParticleJumpProcesses()) {
        compartment.addContent(getXML(pjp));
    }
    for (ParticleProperties pp : param.getParticleProperties()) {
        compartment.addContent(getXML(pp));
    }
    return compartment;
}