Search in sources :

Example 6 with VarIniCondition

use of cbit.vcell.math.VarIniCondition 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;
}
Also used : VarIniCondition(cbit.vcell.math.VarIniCondition) Element(org.jdom.Element) ParticleJumpProcess(cbit.vcell.math.ParticleJumpProcess) ComputeNormalComponentEquation(cbit.vcell.math.ComputeNormalComponentEquation) VolumeRegionEquation(cbit.vcell.math.VolumeRegionEquation) PdeEquation(cbit.vcell.math.PdeEquation) ComputeMembraneMetricEquation(cbit.vcell.math.ComputeMembraneMetricEquation) OdeEquation(cbit.vcell.math.OdeEquation) ComputeCentroidComponentEquation(cbit.vcell.math.ComputeCentroidComponentEquation) MembraneRegionEquation(cbit.vcell.math.MembraneRegionEquation) Equation(cbit.vcell.math.Equation) JumpProcess(cbit.vcell.math.JumpProcess) ParticleJumpProcess(cbit.vcell.math.ParticleJumpProcess) ParticleProperties(cbit.vcell.math.ParticleProperties) BoundaryConditionSpec(cbit.vcell.math.SubDomain.BoundaryConditionSpec)

Example 7 with VarIniCondition

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

the class MathDescPanel method refreshEquations.

/**
 * This method was created by a SmartGuide.
 */
private void refreshEquations() throws ExpressionException {
    TreePath selectionPath = getselectionModel1().getSelectionPath();
    if (selectionPath != null) {
        BioModelNode node = (BioModelNode) selectionPath.getLastPathComponent();
        Object[] objectPath = node.getUserObjectPath();
        Object lastPathObject = objectPath[objectPath.length - 1];
        // for non-stochastic
        Expression[] expArray = null;
        // for stochastic
        String[] expString = null;
        getExpressionCanvas1().setStrings(null);
        getExpressionCanvas1().setExpressions(null);
        if (lastPathObject instanceof Function) {
            Function function = (Function) lastPathObject;
            Expression exp = function.getExpression();
            if (getJRadioButtonValue().isSelected()) {
                exp.bindExpression(getMathDescription());
            } else {
                exp.bindExpression(null);
            }
            exp = exp.flatten();
            // stochastic variable initial value must be an int, therefore we must make it here.
            if (getMathDescription().isNonSpatialStoch() && function.getName().startsWith(DiffEquMathMapping.MATH_FUNC_SUFFIX_SPECIES_INIT_CONC_UNIT_PREFIX)) {
                try {
                    double value = exp.evaluateConstant();
                    expString = new String[] { function.getName() + " = " + Math.round(value) };
                    getExpressionCanvas1().setStrings(expString);
                } catch (Exception e) {
                    expArray = new Expression[] { Expression.assign(new Expression(function.getName()), exp) };
                    getExpressionCanvas1().setExpressions(expArray);
                }
            } else {
                expArray = new Expression[] { Expression.assign(new Expression(function.getName()), exp) };
                getExpressionCanvas1().setExpressions(expArray);
            }
        } else if (lastPathObject instanceof Constant) {
            Constant constant = (Constant) lastPathObject;
            Expression exp = constant.getExpression();
            if (getJRadioButtonValue().isSelected()) {
                exp.bindExpression(getMathDescription());
            } else {
                exp.bindExpression(null);
            }
            exp = exp.flatten();
            // stochastic variable initial value must be an int, therefore we must make it here.
            if (getMathDescription().isNonSpatialStoch() && constant.getName().startsWith(DiffEquMathMapping.MATH_FUNC_SUFFIX_SPECIES_INIT_CONC_UNIT_PREFIX)) {
                try {
                    double value = exp.evaluateConstant();
                    expString = new String[] { constant.getName() + " = " + Math.round(value) };
                    getExpressionCanvas1().setStrings(expString);
                } catch (Exception e) {
                    expArray = new Expression[] { Expression.assign(new Expression(constant.getName()), exp) };
                    getExpressionCanvas1().setExpressions(expArray);
                }
            } else {
                expArray = new Expression[] { Expression.assign(new Expression(constant.getName()), exp) };
                getExpressionCanvas1().setExpressions(expArray);
            }
        } else if (lastPathObject instanceof Equation) {
            Equation equ = (Equation) lastPathObject;
            Enumeration<Expression> enum_equ = equ.getTotalExpressions();
            Vector<Expression> expList = new Vector<Expression>();
            while (enum_equ.hasMoreElements()) {
                Expression exp = new Expression(enum_equ.nextElement());
                if (getJRadioButtonValue().isSelected()) {
                    exp.bindExpression(getMathDescription());
                } else {
                    exp.bindExpression(null);
                }
                expList.addElement(exp.flatten());
            }
            expArray = new Expression[expList.size()];
            expList.copyInto(expArray);
            getExpressionCanvas1().setExpressions(expArray);
        } else if (// added Sept 29, 2006
        lastPathObject instanceof VarIniCondition) {
            VarIniCondition varIni = (VarIniCondition) lastPathObject;
            expString = new String[] { varIni.toString() };
            getExpressionCanvas1().setStrings(expString);
        } else if (// added Sept 29, 2006
        lastPathObject instanceof Action) {
            Action action = (Action) lastPathObject;
            expString = new String[] { action.toString() };
            getExpressionCanvas1().setStrings(expString);
        } else if (// added Sept 29, 2006
        lastPathObject instanceof String) {
            // to check if the string is the probability string
            int index1 = ((String) lastPathObject).indexOf("probability_rate");
            int index2 = ((String) lastPathObject).indexOf("=");
            if (// has probability eqution,but don't calculate prob string, since it is not a constant value.
            (index1 == 0) && ((index2 + 2) < ((String) lastPathObject).length())) {
                expString = new String[] { (String) lastPathObject };
                getExpressionCanvas1().setStrings(expString);
            }
        }
    } else {
        getExpressionCanvas1().setExpressions(null);
        getExpressionCanvas1().setStrings(null);
        getExpressionCanvas1().repaint();
    }
// Variable var = getMathDescription().getVariable(varName);
// SubDomain subDomain = getMathDescription().getSubDomain(subDomainName);
// Equation equ = null;
// if (var != null && subDomain != null){
// equ = subDomain.getEquation(var);
// }
// if (equ == null){
// if (subDomain instanceof MembraneSubDomain && var instanceof VolVariable){
// equ = ((MembraneSubDomain)subDomain).getJumpCondition((VolVariable)var);
// }
// }
// if (equ == null){
// Expression exp = null;
// String title = "hello";
// if (var instanceof Constant || var instanceof Function){
// exp = var.getExpression();
// title = "Expression for "+var.getName();
// }
// if (getJRadioButtonValue().isSelected()){
// exp.bindExpression(getMathDescription());
// }else{
// exp.bindExpression(null);
// }
// getExpressionCanvas1().setExpression(exp.flatten());
// getEquationTitleLabel().setText(title);
// }else{
// Enumeration enum_equ = equ.getTotalExpressions();
// Vector expList = new Vector();
// while (enum_equ.hasMoreElements()){
// Expression exp = new Expression((Expression)enum_equ.nextElement());
// if (getJRadioButtonValue().isSelected()){
// exp.bindExpression(getMathDescription());
// }else{
// exp.bindExpression(null);
// }
// expList.addElement(exp.flatten());
// }
// Expression expArray[] = new Expression[expList.size()];
// expList.copyInto(expArray);
// getExpressionCanvas1().setExpressions(expArray);
// getEquationTitleLabel().setText("Equations for "+var.getName()+" within "+subDomain.getName());
// }
// } catch (java.lang.Throwable ivjExc) {
// handleException(ivjExc);
// getEquationTitleLabel().setText("Equations for "+varName+" within "+subDomainName+" *** "+ivjExc.getMessage());
// }
// }
// }else{
// getEquationTitleLabel().setText("no variable selected");
// getExpressionCanvas1().setExpression((Expression)null);
// }
}
Also used : VarIniCondition(cbit.vcell.math.VarIniCondition) Action(cbit.vcell.math.Action) Constant(cbit.vcell.math.Constant) Equation(cbit.vcell.math.Equation) BioModelNode(cbit.vcell.desktop.BioModelNode) ExpressionException(cbit.vcell.parser.ExpressionException) Function(cbit.vcell.math.Function) TreePath(javax.swing.tree.TreePath) Expression(cbit.vcell.parser.Expression) Vector(java.util.Vector)

Example 8 with VarIniCondition

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

the class MathDescriptionTreeModel method createBaseTree.

/**
 * Insert the method's description here.
 * Creation date: (11/28/00 1:06:51 PM)
 * @return cbit.vcell.desktop.BioModelNode
 * @param docManager cbit.vcell.clientdb.DocumentManager
 */
private BioModelNode createBaseTree() {
    if (getMathDescription() == null) {
        return new BioModelNode(" ", false);
    }
    // 
    // make root node
    // 
    BioModelNode rootNode = new BioModelNode("math description", true);
    // 
    // create subTree of Constants
    // 
    BioModelNode constantRootNode = new BioModelNode("constants", true);
    Enumeration<Constant> enum1 = getMathDescription().getConstants();
    while (enum1.hasMoreElements()) {
        Constant constant = enum1.nextElement();
        BioModelNode constantNode = new BioModelNode(constant, false);
        constantRootNode.add(constantNode);
    }
    if (constantRootNode.getChildCount() > 0) {
        rootNode.add(constantRootNode);
    }
    // 
    // create subTree of Functions
    // 
    BioModelNode functionRootNode = new BioModelNode("functions", true);
    Enumeration<Function> enum2 = getMathDescription().getFunctions();
    while (enum2.hasMoreElements()) {
        Function function = enum2.nextElement();
        BioModelNode functionNode = new BioModelNode(function, false);
        functionRootNode.add(functionNode);
    }
    if (functionRootNode.getChildCount() > 0) {
        rootNode.add(functionRootNode);
    }
    // 
    // create subTree of VolumeSubDomains and MembraneSubDomains
    // 
    BioModelNode volumeRootNode = new BioModelNode("volume domains", true);
    BioModelNode membraneRootNode = new BioModelNode("membrane domains", true);
    Enumeration<SubDomain> enum3 = getMathDescription().getSubDomains();
    while (enum3.hasMoreElements()) {
        SubDomain subDomain = enum3.nextElement();
        if (subDomain instanceof cbit.vcell.math.CompartmentSubDomain) {
            CompartmentSubDomain volumeSubDomain = (CompartmentSubDomain) subDomain;
            BioModelNode volumeSubDomainNode = new BioModelNode(volumeSubDomain, true);
            if (// stochastic subtree
            getMathDescription().isNonSpatialStoch()) {
                // add stoch variable initial conditions
                BioModelNode varIniConditionNode = new BioModelNode("variable_initial_conditions", true);
                for (VarIniCondition varIni : volumeSubDomain.getVarIniConditions()) {
                    BioModelNode varIniNode = new BioModelNode(varIni, false);
                    varIniConditionNode.add(varIniNode);
                }
                volumeSubDomainNode.add(varIniConditionNode);
                // add jump processes
                for (JumpProcess jp : volumeSubDomain.getJumpProcesses()) {
                    BioModelNode jpNode = new BioModelNode(jp, true);
                    // add probability rate.
                    String probRate = "P_" + jp.getName();
                    BioModelNode prNode = new BioModelNode("probability_rate = " + probRate, false);
                    jpNode.add(prNode);
                    // add Actions
                    Enumeration<Action> actions = Collections.enumeration(jp.getActions());
                    while (actions.hasMoreElements()) {
                        Action action = actions.nextElement();
                        BioModelNode actionNode = new BioModelNode(action, false);
                        jpNode.add(actionNode);
                    }
                    volumeSubDomainNode.add(jpNode);
                }
            } else // non-stochastic subtree
            {
                // 
                // add equation children
                // 
                Enumeration<Equation> eqnEnum = volumeSubDomain.getEquations();
                while (eqnEnum.hasMoreElements()) {
                    Equation equation = eqnEnum.nextElement();
                    BioModelNode equationNode = new BioModelNode(equation, false);
                    volumeSubDomainNode.add(equationNode);
                }
                // 
                // add fast system
                // 
                FastSystem fastSystem = volumeSubDomain.getFastSystem();
                if (fastSystem != null) {
                    BioModelNode fsNode = new BioModelNode(fastSystem, true);
                    Enumeration<FastInvariant> enumFI = fastSystem.getFastInvariants();
                    while (enumFI.hasMoreElements()) {
                        FastInvariant fi = enumFI.nextElement();
                        fsNode.add(new BioModelNode(fi, false));
                    }
                    Enumeration<FastRate> enumFR = fastSystem.getFastRates();
                    while (enumFR.hasMoreElements()) {
                        FastRate fr = enumFR.nextElement();
                        fsNode.add(new BioModelNode(fr, false));
                    }
                    volumeSubDomainNode.add(fsNode);
                }
            }
            volumeRootNode.add(volumeSubDomainNode);
        } else if (subDomain instanceof MembraneSubDomain) {
            MembraneSubDomain membraneSubDomain = (MembraneSubDomain) subDomain;
            BioModelNode membraneSubDomainNode = new BioModelNode(membraneSubDomain, true);
            // 
            // add equation children
            // 
            Enumeration<Equation> eqnEnum = membraneSubDomain.getEquations();
            while (eqnEnum.hasMoreElements()) {
                Equation equation = eqnEnum.nextElement();
                BioModelNode equationNode = new BioModelNode(equation, false);
                membraneSubDomainNode.add(equationNode);
            }
            // 
            // add jump condition children
            // 
            Enumeration<JumpCondition> jcEnum = membraneSubDomain.getJumpConditions();
            while (jcEnum.hasMoreElements()) {
                JumpCondition jumpCondition = jcEnum.nextElement();
                BioModelNode jcNode = new BioModelNode(jumpCondition, false);
                membraneSubDomainNode.add(jcNode);
            }
            // 
            // add fast system
            // 
            FastSystem fastSystem = membraneSubDomain.getFastSystem();
            if (fastSystem != null) {
                BioModelNode fsNode = new BioModelNode(fastSystem, true);
                Enumeration<FastInvariant> enumFI = fastSystem.getFastInvariants();
                while (enumFI.hasMoreElements()) {
                    FastInvariant fi = enumFI.nextElement();
                    fsNode.add(new BioModelNode(fi, false));
                }
                Enumeration<FastRate> enumFR = fastSystem.getFastRates();
                while (enumFR.hasMoreElements()) {
                    FastRate fr = enumFR.nextElement();
                    fsNode.add(new BioModelNode(fr, false));
                }
                membraneSubDomainNode.add(fsNode);
            }
            membraneRootNode.add(membraneSubDomainNode);
        }
    }
    if (volumeRootNode.getChildCount() > 0) {
        rootNode.add(volumeRootNode);
    }
    if (membraneRootNode.getChildCount() > 0) {
        rootNode.add(membraneRootNode);
    }
    return rootNode;
}
Also used : JumpCondition(cbit.vcell.math.JumpCondition) VarIniCondition(cbit.vcell.math.VarIniCondition) Action(cbit.vcell.math.Action) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) Enumeration(java.util.Enumeration) Constant(cbit.vcell.math.Constant) Equation(cbit.vcell.math.Equation) FastRate(cbit.vcell.math.FastRate) BioModelNode(cbit.vcell.desktop.BioModelNode) FastInvariant(cbit.vcell.math.FastInvariant) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) Function(cbit.vcell.math.Function) FastSystem(cbit.vcell.math.FastSystem) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) JumpProcess(cbit.vcell.math.JumpProcess)

Example 9 with VarIniCondition

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

the class StochMathMapping_4_8 method refreshMathDescription.

/**
 * set up a math description based on current simulationContext.
 */
private void refreshMathDescription() throws MappingException, MatrixException, MathException, ExpressionException, ModelException {
    // use local variable instead of using getter all the time.
    SimulationContext simContext = getSimulationContext();
    // local structure mapping list
    StructureMapping[] structureMappings = simContext.getGeometryContext().getStructureMappings();
    // We have to check if all the reactions are able to tranform to stochastic jump processes before generating the math.
    String stochChkMsg = simContext.getModel().isValidForStochApp();
    if (!(stochChkMsg.equals(""))) {
        throw new ModelException("Problem updating math description: " + simContext.getName() + "\n" + stochChkMsg);
    }
    // All sizes must be set for new ODE models and ratios must be set for old ones.
    simContext.checkValidity();
    // 
    // verify that all structures are mapped to subvolumes and all subvolumes are mapped to a structure
    // 
    Structure[] structures = simContext.getGeometryContext().getModel().getStructures();
    for (int i = 0; i < structures.length; i++) {
        StructureMapping sm = simContext.getGeometryContext().getStructureMapping(structures[i]);
        if (sm == null || (sm instanceof FeatureMapping && getSubVolume(((FeatureMapping) sm)) == null)) {
            throw new MappingException("model structure '" + structures[i].getName() + "' not mapped to a geometry subVolume");
        }
        if (sm != null && (sm instanceof MembraneMapping) && ((MembraneMapping) sm).getVolumeFractionParameter() != null) {
            Expression volFractExp = ((MembraneMapping) sm).getVolumeFractionParameter().getExpression();
            try {
                if (volFractExp != null) {
                    double volFract = volFractExp.evaluateConstant();
                    if (volFract >= 1.0) {
                        throw new MappingException("model structure '" + (getSimulationContext().getModel().getStructureTopology().getInsideFeature(((MembraneMapping) sm).getMembrane()).getName() + "' has volume fraction >= 1.0"));
                    }
                }
            } catch (ExpressionException e) {
                e.printStackTrace(System.out);
            }
        }
    }
    SubVolume[] subVolumes = simContext.getGeometryContext().getGeometry().getGeometrySpec().getSubVolumes();
    for (int i = 0; i < subVolumes.length; i++) {
        if (getStructures(subVolumes[i]) == null || getStructures(subVolumes[i]).length == 0) {
            throw new MappingException("geometry subVolume '" + subVolumes[i].getName() + "' not mapped from a model structure");
        }
    }
    // 
    // gather only those reactionSteps that are not "excluded"
    // 
    ReactionSpec[] reactionSpecs = simContext.getReactionContext().getReactionSpecs();
    Vector<ReactionStep> rsList = new Vector<ReactionStep>();
    for (int i = 0; i < reactionSpecs.length; i++) {
        if (reactionSpecs[i].isExcluded() == false) {
            rsList.add(reactionSpecs[i].getReactionStep());
        }
    }
    ReactionStep[] reactionSteps = new ReactionStep[rsList.size()];
    rsList.copyInto(reactionSteps);
    // 
    for (int i = 0; i < reactionSteps.length; i++) {
        Kinetics.UnresolvedParameter[] unresolvedParameters = reactionSteps[i].getKinetics().getUnresolvedParameters();
        if (unresolvedParameters != null && unresolvedParameters.length > 0) {
            StringBuffer buffer = new StringBuffer();
            for (int j = 0; j < unresolvedParameters.length; j++) {
                if (j > 0) {
                    buffer.append(", ");
                }
                buffer.append(unresolvedParameters[j].getName());
            }
            throw new MappingException(reactionSteps[i].getDisplayType() + " '" + reactionSteps[i].getName() + "' contains unresolved identifier(s): " + buffer);
        }
    }
    // 
    // create new MathDescription (based on simContext's previous MathDescription if possible)
    // 
    MathDescription oldMathDesc = simContext.getMathDescription();
    mathDesc = null;
    if (oldMathDesc != null) {
        if (oldMathDesc.getVersion() != null) {
            mathDesc = new MathDescription(oldMathDesc.getVersion());
        } else {
            mathDesc = new MathDescription(oldMathDesc.getName());
        }
    } else {
        mathDesc = new MathDescription(simContext.getName() + "_generated");
    }
    // 
    // temporarily place all variables in a hashtable (before binding) and discarding duplicates
    // 
    VariableHash varHash = new VariableHash();
    // 
    // conversion factors
    // 
    Model model = simContext.getModel();
    ModelUnitSystem modelUnitSystem = model.getUnitSystem();
    varHash.addVariable(new Constant(getMathSymbol(model.getKMOLE(), null), getIdentifierSubstitutions(model.getKMOLE().getExpression(), model.getKMOLE().getUnitDefinition(), null)));
    varHash.addVariable(new Constant(getMathSymbol(model.getN_PMOLE(), null), getIdentifierSubstitutions(model.getN_PMOLE().getExpression(), model.getN_PMOLE().getUnitDefinition(), null)));
    varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT().getExpression(), model.getFARADAY_CONSTANT().getUnitDefinition(), null)));
    varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT_NMOLE(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT_NMOLE().getExpression(), model.getFARADAY_CONSTANT_NMOLE().getUnitDefinition(), null)));
    varHash.addVariable(new Constant(getMathSymbol(model.getGAS_CONSTANT(), null), getIdentifierSubstitutions(model.getGAS_CONSTANT().getExpression(), model.getGAS_CONSTANT().getUnitDefinition(), null)));
    varHash.addVariable(new Constant(getMathSymbol(model.getTEMPERATURE(), null), getIdentifierSubstitutions(new Expression(simContext.getTemperatureKelvin()), model.getTEMPERATURE().getUnitDefinition(), null)));
    Enumeration<SpeciesContextMapping> enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = enum1.nextElement();
        if (scm.getVariable() instanceof StochVolVariable) {
            varHash.addVariable(scm.getVariable());
        }
    }
    // 
    // add rate term for all reactions
    // add current source terms for each reaction step in a membrane
    // 
    /*for (int i = 0; i < reactionSteps.length; i++){
			boolean bAllReactionParticipantsFixed = true;
			ReactionParticipant rp_Array[] = reactionSteps[i].getReactionParticipants();
			for (int j = 0; j < rp_Array.length; j++) {
				SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(rp_Array[j].getSpeciesContext());
				if (!(rp_Array[j] instanceof Catalyst) && !scs.isConstant()){
					bAllReactionParticipantsFixed = false;  // found at least one reactionParticipant that is not fixed and needs this rate
				}
			}
			StructureMapping sm = simContext.getGeometryContext().getStructureMapping(reactionSteps[i].getStructure());
		}---don't think it's useful, isn't it?*/
    // deals with model parameters
    ModelParameter[] modelParameters = simContext.getModel().getModelParameters();
    for (int j = 0; j < modelParameters.length; j++) {
        Expression expr = getSubstitutedExpr(modelParameters[j].getExpression(), true, false);
        expr = getIdentifierSubstitutions(expr, modelParameters[j].getUnitDefinition(), null);
        varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], null), expr));
    }
    // added July 2009, ElectricalStimulusParameter electric mapping tab
    ElectricalStimulus[] elecStimulus = simContext.getElectricalStimuli();
    if (elecStimulus.length > 0) {
        throw new MappingException("Modles with electrophysiology are not supported for stochastic applications.");
    }
    for (int j = 0; j < structureMappings.length; j++) {
        if (structureMappings[j] instanceof MembraneMapping) {
            MembraneMapping memMapping = (MembraneMapping) structureMappings[j];
            Parameter initialVoltageParm = memMapping.getInitialVoltageParameter();
            try {
                Expression exp = initialVoltageParm.getExpression();
                exp.evaluateConstant();
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(memMapping.getMembrane().getMembraneVoltage(), memMapping), getIdentifierSubstitutions(memMapping.getInitialVoltageParameter().getExpression(), memMapping.getInitialVoltageParameter().getUnitDefinition(), memMapping)));
            } catch (ExpressionException e) {
                e.printStackTrace(System.out);
                throw new MappingException("Membrane initial voltage: " + initialVoltageParm.getName() + " cannot be evaluated as constant.");
            }
        }
    }
    // 
    for (int j = 0; j < reactionSteps.length; j++) {
        ReactionStep rs = reactionSteps[j];
        if (simContext.getReactionContext().getReactionSpec(rs).isExcluded()) {
            continue;
        }
        if (rs.getKinetics() instanceof LumpedKinetics) {
            throw new RuntimeException("Lumped Kinetics not yet supported for Stochastic Math Generation");
        }
        Kinetics.KineticsParameter[] parameters = rs.getKinetics().getKineticsParameters();
        StructureMapping sm = simContext.getGeometryContext().getStructureMapping(rs.getStructure());
        if (parameters != null) {
            for (int i = 0; i < parameters.length; i++) {
                if ((parameters[i].getRole() == Kinetics.ROLE_CurrentDensity) && (parameters[i].getExpression() == null || parameters[i].getExpression().isZero())) {
                    continue;
                }
                // don't add rate, we'll do it later when creating the jump processes
                if (parameters[i].getRole() != Kinetics.ROLE_ReactionRate) {
                    Expression expr = getSubstitutedExpr(parameters[i].getExpression(), true, false);
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[i], sm), getIdentifierSubstitutions(expr, parameters[i].getUnitDefinition(), sm)));
                }
            }
        }
    }
    // the parameter "Size" is already put into mathsymbolmapping in refreshSpeciesContextMapping()
    for (int i = 0; i < structureMappings.length; i++) {
        StructureMapping sm = structureMappings[i];
        StructureMapping.StructureMappingParameter parm = sm.getParameterFromRole(StructureMapping.ROLE_Size);
        if (parm.getExpression() != null) {
            try {
                double value = parm.getExpression().evaluateConstant();
                varHash.addVariable(new Constant(getMathSymbol(parm, sm), new Expression(value)));
            } catch (ExpressionException e) {
                // varHash.addVariable(new Function(getMathSymbol0(parm,sm),getIdentifierSubstitutions(parm.getExpression(),parm.getUnitDefinition(),sm)));
                e.printStackTrace(System.out);
                throw new MappingException("Size of structure:" + sm.getNameScope().getName() + " cannot be evaluated as constant.");
            }
        }
    }
    // 
    // species initial values (either function or constant)
    // 
    SpeciesContextSpec[] speciesContextSpecs = simContext.getReactionContext().getSpeciesContextSpecs();
    for (int i = 0; i < speciesContextSpecs.length; i++) {
        // can be concentration or amount
        SpeciesContextSpec.SpeciesContextSpecParameter initParam = null;
        Expression iniExp = null;
        StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
        if (speciesContextSpecs[i].getInitialConcentrationParameter() != null && speciesContextSpecs[i].getInitialConcentrationParameter().getExpression() != null) {
            // use concentration, need to set up amount functions
            initParam = speciesContextSpecs[i].getInitialConcentrationParameter();
            iniExp = initParam.getExpression();
            iniExp = getSubstitutedExpr(iniExp, true, !speciesContextSpecs[i].isConstant());
            // now create the appropriate function or Constant for the speciesContextSpec.
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(initParam, sm), getIdentifierSubstitutions(iniExp, initParam.getUnitDefinition(), sm)));
            // add function for initial amount
            SpeciesContextSpec.SpeciesContextSpecParameter initAmountParam = speciesContextSpecs[i].getInitialCountParameter();
            Expression iniAmountExp = getExpressionConcToAmt(new Expression(initParam, getNameScope()), speciesContextSpecs[i].getSpeciesContext());
            // iniAmountExp.bindExpression(this);
            varHash.addVariable(new Function(getMathSymbol(initAmountParam, sm), getIdentifierSubstitutions(iniAmountExp, initAmountParam.getUnitDefinition(), sm), nullDomain));
        } else if (speciesContextSpecs[i].getInitialCountParameter() != null && speciesContextSpecs[i].getInitialCountParameter().getExpression() != null) {
            // use amount
            initParam = speciesContextSpecs[i].getInitialCountParameter();
            iniExp = initParam.getExpression();
            iniExp = getSubstitutedExpr(iniExp, false, !speciesContextSpecs[i].isConstant());
            // now create the appropriate function or Constant for the speciesContextSpec.
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(initParam, sm), getIdentifierSubstitutions(iniExp, initParam.getUnitDefinition(), sm)));
        }
        // add spConcentration (concentration of species) to varHash as function or constant
        SpeciesConcentrationParameter spConcParam = getSpeciesConcentrationParameter(speciesContextSpecs[i].getSpeciesContext());
        varHash.addVariable(newFunctionOrConstant(getMathSymbol(spConcParam, sm), getIdentifierSubstitutions(spConcParam.getExpression(), spConcParam.getUnitDefinition(), sm)));
    }
    // 
    // constant species (either function or constant)
    // 
    enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
        if (scm.getVariable() instanceof Constant) {
            varHash.addVariable(scm.getVariable());
        }
    }
    // 
    if (simContext.getGeometryContext().getGeometry() != null) {
        try {
            mathDesc.setGeometry(simContext.getGeometryContext().getGeometry());
        } catch (java.beans.PropertyVetoException e) {
            e.printStackTrace(System.out);
            throw new MappingException("failure setting geometry " + e.getMessage());
        }
    } else {
        throw new MappingException("geometry must be defined");
    }
    // 
    // functions: species which is not a variable, but has dependency expression
    // 
    enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
        if (scm.getVariable() == null && scm.getDependencyExpression() != null) {
            StructureMapping sm = simContext.getGeometryContext().getStructureMapping(scm.getSpeciesContext().getStructure());
            Expression exp = scm.getDependencyExpression();
            exp.bindExpression(this);
            SpeciesCountParameter spCountParam = getSpeciesCountParameter(scm.getSpeciesContext());
            varHash.addVariable(new Function(getMathSymbol(spCountParam, sm), getIdentifierSubstitutions(exp, spCountParam.getUnitDefinition(), sm), nullDomain));
        }
    }
    // 
    // create subDomains
    // 
    SubDomain subDomain = null;
    subVolumes = simContext.getGeometryContext().getGeometry().getGeometrySpec().getSubVolumes();
    for (int j = 0; j < subVolumes.length; j++) {
        SubVolume subVolume = (SubVolume) subVolumes[j];
        // 
        // get priority of subDomain
        // 
        int priority;
        Feature spatialFeature = getResolvedFeature(subVolume);
        if (spatialFeature == null) {
            if (simContext.getGeometryContext().getGeometry().getDimension() > 0) {
                throw new MappingException("no compartment (in Physiology) is mapped to subdomain '" + subVolume.getName() + "' (in Geometry)");
            } else {
                priority = CompartmentSubDomain.NON_SPATIAL_PRIORITY;
            }
        } else {
            // now does not have to match spatial feature, *BUT* needs to be unique
            priority = j;
        }
        subDomain = new CompartmentSubDomain(subVolume.getName(), priority);
        mathDesc.addSubDomain(subDomain);
    }
    // ReactionSpec[] reactionSpecs = simContext.getReactionContext().getReactionSpecs();---need to take a look here!
    for (int i = 0; i < reactionSpecs.length; i++) {
        if (reactionSpecs[i].isExcluded()) {
            continue;
        }
        // get the reaction
        ReactionStep reactionStep = reactionSpecs[i].getReactionStep();
        Kinetics kinetics = reactionStep.getKinetics();
        // the structure where reaction happens
        StructureMapping sm = simContext.getGeometryContext().getStructureMapping(reactionStep.getStructure());
        // create symbol table for jump process based on reactionStep and structure mapping
        // final ReactionStep finalRS = reactionStep;
        // final StructureMapping finalSM = sm;
        // SymbolTable symTable = new SymbolTable(){
        // public SymbolTableEntry getEntry(String identifierString) throws ExpressionBindingException {
        // SymbolTableEntry ste = finalRS.getEntry(identifierString);
        // if(ste == null)
        // {
        // ste = finalSM.getEntry(identifierString);
        // }
        // return ste;
        // }
        // };
        // Different ways to deal with simple reactions and flux reactions
        // probability parameter from modelUnitSystem
        VCUnitDefinition probabilityParamUnit = modelUnitSystem.getStochasticSubstanceUnit().divideBy(modelUnitSystem.getTimeUnit());
        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)) {
                forwardRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KForward).getExpression();
                reverseRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KReverse).getExpression();
            } else if (kinetics.getKineticsDescription().equals(KineticsDescription.General)) {
                Expression rateExp = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate).getExpression();
                MassActionSolver.MassActionFunction maFunc = MassActionSolver.solveMassAction(null, null, 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 (kinetics.getKineticsDescription().getName().compareTo(KineticsDescription.HMM_irreversible.getName())==0)
			    {
				    forwardRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Km).getExpression();
				}
			    else if (kinetics.getKineticsDescription().getName().compareTo(KineticsDescription.HMM_reversible.getName())==0)
			    {
					forwardRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KmFwd).getExpression();
					reverseRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KmRev).getExpression();
				}*/
            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 mass actions
                exp = getProbabilityRate(reactionStep, true);
                // bind symbol table before substitute identifiers in the reaction step
                exp.bindExpression(this);
                MathMapping_4_8.ProbabilityParameter probParm = null;
                try {
                    probParm = addProbabilityParameter("P_" + jpName, exp, MathMapping_4_8.PARAMETER_ROLE_P, probabilityParamUnit, reactionSpecs[i]);
                } catch (PropertyVetoException pve) {
                    pve.printStackTrace();
                    throw new MappingException(pve.getMessage());
                }
                // add probability to function or constant
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(probParm, sm), getIdentifierSubstitutions(exp, probabilityParamUnit, sm)));
                JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probParm, sm)));
                // 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 = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression("-" + String.valueOf(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 = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", 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()) + "_reverse";
                Expression exp = null;
                // reactions are mass actions
                exp = getProbabilityRate(reactionStep, false);
                // bind symbol table before substitute identifiers in the reaction step
                exp.bindExpression(this);
                MathMapping_4_8.ProbabilityParameter probRevParm = null;
                try {
                    probRevParm = addProbabilityParameter("P_" + jpName, exp, MathMapping_4_8.PARAMETER_ROLE_P_reverse, probabilityParamUnit, reactionSpecs[i]);
                } catch (PropertyVetoException pve) {
                    pve.printStackTrace();
                    throw new MappingException(pve.getMessage());
                }
                // add probability to function or constant
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(probRevParm, sm), getIdentifierSubstitutions(exp, probabilityParamUnit, sm)));
                JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probRevParm, sm)));
                // 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 = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", 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 = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression("-" + String.valueOf(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)) {
                Expression fluxRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate).getExpression();
                // we have to pass the math description para to flux solver, coz somehow math description in simulation context is not updated.
                MassActionSolver.MassActionFunction fluxFunc = MassActionSolver.solveMassAction(null, null, fluxRate, (FluxReaction) reactionStep);
                // create jump process for forward flux if it exists.
                if (fluxFunc.getForwardRate() != null && !fluxFunc.getForwardRate().isZero()) {
                    // jump process name
                    // +"_reverse";
                    String jpName = TokenMangler.mangleToSName(reactionStep.getName());
                    // we do it here instead of fluxsolver, coz we need to use getMathSymbol0(), structuremapping...etc.
                    Expression rate = fluxFunc.getForwardRate();
                    // get species expression (depend on structure, if mem: Species/mem_Size, if vol: species*KMOLE/vol_size)
                    SpeciesContext scOut = fluxFunc.getReactants().get(0).getSpeciesContext();
                    Expression speciesFactor = null;
                    if (scOut.getStructure() instanceof Feature) {
                        Expression exp1 = new Expression(Model.reservedConstantsMap.get(ReservedSymbolRole.KMOLE));
                        Expression exp2 = new Expression(scOut.getStructure().getStructureSize(), getNameScope());
                        speciesFactor = Expression.div(Expression.invert(exp1), exp2);
                    } else {
                        throw new MappingException("Species involved in a flux have to be volume species.");
                    }
                    Expression speciesExp = Expression.mult(speciesFactor, new Expression(scOut, getNameScope()));
                    // get probability expression by adding factor to rate (rate: rate*size_mem/KMOLE)
                    Expression expr1 = Expression.mult(rate, speciesExp);
                    Expression numeratorExpr = Expression.mult(expr1, new Expression(sm.getStructure().getStructureSize(), getNameScope()));
                    Expression exp = new Expression(Model.reservedConstantsMap.get(ReservedSymbolRole.KMOLE));
                    Expression probExp = Expression.mult(numeratorExpr, exp);
                    // bind symbol table before substitute identifiers in the reaction step
                    probExp.bindExpression(reactionStep);
                    MathMapping_4_8.ProbabilityParameter probParm = null;
                    try {
                        probParm = addProbabilityParameter("P_" + jpName, probExp, MathMapping_4_8.PARAMETER_ROLE_P, probabilityParamUnit, reactionSpecs[i]);
                    } catch (PropertyVetoException pve) {
                        pve.printStackTrace();
                        throw new MappingException(pve.getMessage());
                    }
                    // add probability to function or constant
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(probParm, sm), getIdentifierSubstitutions(probExp, probabilityParamUnit, sm)));
                    JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probParm, sm)));
                    // actions
                    Action action = null;
                    SpeciesContext sc = fluxFunc.getReactants().get(0).getSpeciesContext();
                    if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
                        action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression(-1));
                        jp.addAction(action);
                    }
                    sc = fluxFunc.getProducts().get(0).getSpeciesContext();
                    if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
                        action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression(1));
                        jp.addAction(action);
                    }
                    subDomain.addJumpProcess(jp);
                }
                if (fluxFunc.getReverseRate() != null && !fluxFunc.getReverseRate().isZero()) {
                    // jump process name
                    String jpName = TokenMangler.mangleToSName(reactionStep.getName()) + "_reverse";
                    Expression rate = fluxFunc.getReverseRate();
                    // get species expression (depend on structure, if mem: Species/mem_Size, if vol: species*KMOLE/vol_size)
                    SpeciesContext scIn = fluxFunc.getProducts().get(0).getSpeciesContext();
                    Expression speciesFactor = null;
                    if (scIn.getStructure() instanceof Feature) {
                        Expression exp1 = new Expression(Model.reservedConstantsMap.get(ReservedSymbolRole.KMOLE));
                        Expression exp2 = new Expression(scIn.getStructure().getStructureSize(), getNameScope());
                        speciesFactor = Expression.div(Expression.invert(exp1), exp2);
                    } else {
                        throw new MappingException("Species involved in a flux have to be volume species.");
                    }
                    Expression speciesExp = Expression.mult(speciesFactor, new Expression(scIn, getNameScope()));
                    // get probability expression by adding factor to rate (rate: rate*size_mem/KMOLE)
                    Expression expr1 = Expression.mult(rate, speciesExp);
                    Expression numeratorExpr = Expression.mult(expr1, new Expression(sm.getStructure().getStructureSize(), getNameScope()));
                    Expression exp = new Expression(Model.reservedConstantsMap.get(ReservedSymbolRole.KMOLE));
                    Expression probRevExp = Expression.mult(numeratorExpr, exp);
                    // bind symbol table before substitute identifiers in the reaction step
                    probRevExp.bindExpression(reactionStep);
                    MathMapping_4_8.ProbabilityParameter probRevParm = null;
                    try {
                        probRevParm = addProbabilityParameter("P_" + jpName, probRevExp, MathMapping_4_8.PARAMETER_ROLE_P_reverse, probabilityParamUnit, reactionSpecs[i]);
                    } catch (PropertyVetoException pve) {
                        pve.printStackTrace();
                        throw new MappingException(pve.getMessage());
                    }
                    // add probability to function or constant
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(probRevParm, sm), getIdentifierSubstitutions(probRevExp, probabilityParamUnit, sm)));
                    JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probRevParm, sm)));
                    // actions
                    Action action = null;
                    SpeciesContext sc = fluxFunc.getReactants().get(0).getSpeciesContext();
                    if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
                        action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression(1));
                        jp.addAction(action);
                    }
                    sc = fluxFunc.getProducts().get(0).getSpeciesContext();
                    if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
                        action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression(-1));
                        jp.addAction(action);
                    }
                    subDomain.addJumpProcess(jp);
                }
            }
        }
    // end of if (simplereaction)...else if(fluxreaction)
    }
    // end of reaction step loop
    // 
    // set Variables to MathDescription all at once with the order resolved by "VariableHash"
    // 
    mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
    // set up variable initial conditions in subDomain
    SpeciesContextSpec[] scSpecs = simContext.getReactionContext().getSpeciesContextSpecs();
    for (int i = 0; i < speciesContextSpecs.length; i++) {
        // get stochastic variable by name
        SpeciesCountParameter spCountParam = getSpeciesCountParameter(speciesContextSpecs[i].getSpeciesContext());
        StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
        String varName = getMathSymbol(spCountParam, sm);
        if (scSpecs[i].isConstant()) {
            continue;
        }
        StochVolVariable var = (StochVolVariable) mathDesc.getVariable(varName);
        // stochastic use initial number of particles
        SpeciesContextSpec.SpeciesContextSpecParameter initParm = scSpecs[i].getInitialCountParameter();
        // stochastic variables initial expression.
        if (initParm != null) {
            VarIniCondition varIni = new VarIniCount(var, new Expression(getMathSymbol(initParm, sm)));
            subDomain.addVarIniCondition(varIni);
        }
    }
    if (!mathDesc.isValid()) {
        throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
    }
}
Also used : VarIniCondition(cbit.vcell.math.VarIniCondition) MembraneMapping(cbit.vcell.mapping.MembraneMapping) LumpedKinetics(cbit.vcell.model.LumpedKinetics) MathDescription(cbit.vcell.math.MathDescription) SpeciesContextMapping(cbit.vcell.mapping.SpeciesContextMapping) Product(cbit.vcell.model.Product) FluxReaction(cbit.vcell.model.FluxReaction) SpeciesContext(cbit.vcell.model.SpeciesContext) SpeciesContextSpec(cbit.vcell.mapping.SpeciesContextSpec) Feature(cbit.vcell.model.Feature) Reactant(cbit.vcell.model.Reactant) ExpressionException(cbit.vcell.parser.ExpressionException) MappingException(cbit.vcell.mapping.MappingException) PropertyVetoException(java.beans.PropertyVetoException) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter) SubVolume(cbit.vcell.geometry.SubVolume) Vector(java.util.Vector) ModelException(cbit.vcell.model.ModelException) ReactionSpec(cbit.vcell.mapping.ReactionSpec) PropertyVetoException(java.beans.PropertyVetoException) ModelParameter(cbit.vcell.model.Model.ModelParameter) VCUnitDefinition(cbit.vcell.units.VCUnitDefinition) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) ReactionStep(cbit.vcell.model.ReactionStep) Kinetics(cbit.vcell.model.Kinetics) LumpedKinetics(cbit.vcell.model.LumpedKinetics) ReactionParticipant(cbit.vcell.model.ReactionParticipant) Action(cbit.vcell.math.Action) VariableHash(cbit.vcell.math.VariableHash) Constant(cbit.vcell.math.Constant) StructureMapping(cbit.vcell.mapping.StructureMapping) Function(cbit.vcell.math.Function) FeatureMapping(cbit.vcell.mapping.FeatureMapping) JumpProcess(cbit.vcell.math.JumpProcess) Structure(cbit.vcell.model.Structure) StochVolVariable(cbit.vcell.math.StochVolVariable) ModelUnitSystem(cbit.vcell.model.ModelUnitSystem) SimpleReaction(cbit.vcell.model.SimpleReaction) VarIniCount(cbit.vcell.math.VarIniCount) SimulationContext(cbit.vcell.mapping.SimulationContext) ElectricalStimulus(cbit.vcell.mapping.ElectricalStimulus) Expression(cbit.vcell.parser.Expression) Model(cbit.vcell.model.Model) ProxyParameter(cbit.vcell.model.ProxyParameter) Parameter(cbit.vcell.model.Parameter) KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter) ModelParameter(cbit.vcell.model.Model.ModelParameter)

Example 10 with VarIniCondition

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

the class NetCDFWriter method writeHybridInputFile.

/**
 * Write the model to a NetCDF file which serves as an input for stoch hybrid simulator.
 * To write to a NetCDF file is a bit complicated. First, we have to create a NetCDF-3
 * file. And then feed in the data.
 * Creation date: (5/22/2007 5:36:03 PM)
 */
public void writeHybridInputFile(String[] parameterNames) throws Exception, cbit.vcell.parser.ExpressionException, IOException, MathException, InvalidRangeException {
    Simulation simulation = simTask.getSimulation();
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    if (initialize()) {
        // we need to get model and control information first
        NetcdfFileWriteable ncfile = NetcdfFileWriteable.createNew(filename, false);
        // Model info. will be extracted from subDomain of mathDescription
        java.util.Enumeration<SubDomain> e = simulation.getMathDescription().getSubDomains();
        // remember we are dealing with compartmental model here. only 1 subdomain.
        SubDomain subDomain = e.nextElement();
        JumpProcess[] reactions = (JumpProcess[]) subDomain.getJumpProcesses().toArray(new JumpProcess[subDomain.getJumpProcesses().size()]);
        // get species variable names
        Variable[] variables = simSymbolTable.getVariables();
        String[] speciesNames = new String[variables.length];
        for (int i = 0; i < variables.length; i++) speciesNames[i] = variables[i].getName();
        // the probabilities for reactions
        Expression[] probs = new Expression[reactions.length];
        for (int i = 0; i < reactions.length; i++) {
            probs[i] = simSymbolTable.substituteFunctions(reactions[i].getProbabilityRate());
            probs[i] = probs[i].flatten();
        }
        VarIniCondition[] varInis = (VarIniCondition[]) subDomain.getVarIniConditions().toArray(new VarIniCondition[subDomain.getVarIniConditions().size()]);
        // the non-constant stoch variables
        Vector<Variable> vars = new Vector<Variable>();
        for (int i = 0; i < varInis.length; i++) {
            if (varInis[i].getVar() instanceof StochVolVariable) {
                vars.addElement(varInis[i].getVar());
            }
        }
        // get reaction rate law types and rate constants
        ReactionRateLaw[] reactionRateLaws = getReactionRateLaws(probs);
        SolverTaskDescription solverTaskDescription = simulation.getSolverTaskDescription();
        TimeBounds timeBounds = solverTaskDescription.getTimeBounds();
        UniformOutputTimeSpec timeSpec = (UniformOutputTimeSpec) solverTaskDescription.getOutputTimeSpec();
        UniformOutputTimeSpec outputTimeSpec = ((UniformOutputTimeSpec) solverTaskDescription.getOutputTimeSpec());
        NonspatialStochSimOptions stochOpt = solverTaskDescription.getStochOpt();
        // create an empty NetCDF-3 file
        // define dimensions
        /* these sizes must match the buffers allocated in corresponding Fortran code -- see globalvariables.f90
			in numerics Hy3S/src directory */
        Dimension numTrial = ncfile.addDimension("NumTrials", (int) stochOpt.getNumOfTrials());
        Dimension numSpecies = ncfile.addDimension("NumSpecies", vars.size());
        Dimension numReactions = ncfile.addDimension("NumReactions", subDomain.getJumpProcesses().size());
        int outPoints = ((int) ((timeBounds.getEndingTime() - timeBounds.getStartingTime()) / outputTimeSpec.getOutputTimeStep())) + 1;
        Dimension numTimePoints = ncfile.addDimension("NumTimePoints", outPoints);
        Dimension numModels = ncfile.addDimension("NumModels", 1);
        Dimension numMaxDepList = ncfile.addDimension("NumMaxDepList", 6);
        Dimension numMaxStoichList = ncfile.addDimension("NumMaxStoichList", 25);
        Dimension stringLen = ncfile.addDimension("StringLen", 72);
        // define variables
        // jms info
        ArrayList<Dimension> dims = new ArrayList<Dimension>();
        dims.add(stringLen);
        if (bMessaging) {
            ncfile.addVariable("JMS_BROKER", DataType.CHAR, dims);
            ncfile.addVariable("JMS_USER", DataType.CHAR, dims);
            ncfile.addVariable("JMS_PASSWORD", DataType.CHAR, dims);
            ncfile.addVariable("JMS_QUEUE", DataType.CHAR, dims);
            ncfile.addVariable("JMS_TOPIC", DataType.CHAR, dims);
            ncfile.addVariable("VCELL_USER", DataType.CHAR, dims);
            ncfile.addVariable("SIMULATION_KEY", DataType.INT, new ArrayList<Dimension>());
            ncfile.addVariable("JOB_INDEX", DataType.INT, new ArrayList<Dimension>());
        }
        // scalars
        ncfile.addVariable("TStart", DataType.DOUBLE, new ArrayList<Dimension>());
        ncfile.addVariable("TEnd", DataType.DOUBLE, new ArrayList<Dimension>());
        ncfile.addVariable("SaveTime", DataType.DOUBLE, new ArrayList<Dimension>());
        ncfile.addVariable("Volume", DataType.DOUBLE, new ArrayList<Dimension>());
        ncfile.addVariable("CellGrowthTime", DataType.DOUBLE, new ArrayList<Dimension>());
        ncfile.addVariable("CellGrowthTimeSD", DataType.DOUBLE, new ArrayList<Dimension>());
        ncfile.addVariable("ExpType", DataType.INT, new ArrayList<Dimension>());
        ncfile.addVariable("LastTrial", DataType.INT, new ArrayList<Dimension>());
        ncfile.addVariable("LastModel", DataType.INT, new ArrayList<Dimension>());
        ncfile.addVariable("MaxNumModels", DataType.INT, new ArrayList<Dimension>());
        ncfile.addVariable("NumModels", DataType.INT, new ArrayList<Dimension>());
        // variables with at least 1 dimension
        ArrayList<Dimension> dimspecies = new ArrayList<Dimension>();
        dimspecies.add(numSpecies);
        ArrayList<Dimension> dimreactions = new ArrayList<Dimension>();
        dimreactions.add(numReactions);
        ncfile.addVariable("SpeciesSplitOnDivision", DataType.INT, dimspecies);
        ncfile.addVariable("SaveSpeciesData", DataType.INT, dimspecies);
        ncfile.addVariable("Reaction_Rate_Laws", DataType.INT, dimreactions);
        ncfile.addVariable("Reaction_DListLen", DataType.INT, dimreactions);
        ncfile.addVariable("Reaction_StoichListLen", DataType.INT, dimreactions);
        ncfile.addVariable("Reaction_OptionalData", DataType.INT, dimreactions);
        dims.clear();
        dims.add(numReactions);
        dims.add(numMaxStoichList);
        ncfile.addVariable("Reaction_StoichCoeff", DataType.INT, dims);
        ncfile.addVariable("Reaction_StoichSpecies", DataType.INT, dims);
        dims.clear();
        dims.add(numReactions);
        dims.add(numMaxDepList);
        ncfile.addVariable("Reaction_DepList", DataType.INT, dims);
        dims.clear();
        dims.add(numReactions);
        dims.add(stringLen);
        ncfile.addVariable("Reaction_names", DataType.CHAR, dims);
        dims.clear();
        dims.add(numSpecies);
        dims.add(stringLen);
        ncfile.addVariable("Species_names", DataType.CHAR, dims);
        ncfile.addVariable("SpeciesIC", DataType.INT, dimspecies);
        dims.clear();
        dims.add(numReactions);
        dims.add(numMaxDepList);
        ncfile.addVariable("Reaction_Rate_Constants", DataType.DOUBLE, dims);
        // create the file
        try {
            ncfile.create();
        } catch (IOException ioe) {
            ioe.printStackTrace(System.err);
            throw new IOException("Error creating hybrid file " + filename + ": " + ioe.getMessage());
        }
        // write data to the NetCDF file
        try {
            // write jms info
            if (bMessaging) {
                ArrayChar.D1 jmsString = new ArrayChar.D1(stringLen.getLength());
                String jmshost = PropertyLoader.getRequiredProperty(PropertyLoader.jmsSimHostExternal);
                // 
                // Used for new REST HTTP messaging api (USE THIS WHEN Hyrbid Solvers are compiled).
                // 
                // String jmsrestport = PropertyLoader.getRequiredProperty(PropertyLoader.jmsRestPortExternal);
                // String jmsurl = jmshost+":"+jmsrestport;
                // 
                // connect to messaging using legacy AMQP protocol instead of new REST api.  Needed for legacy pre-compiled solvers.
                // 
                String jmsport = PropertyLoader.getRequiredProperty(PropertyLoader.jmsSimPortExternal);
                String jmsurl = "failover:(tcp://" + jmshost + ":" + jmsport + ")";
                jmsString.setString(jmsurl);
                ncfile.write("JMS_BROKER", jmsString);
                jmsString.setString(PropertyLoader.getRequiredProperty(PropertyLoader.jmsUser));
                ncfile.write("JMS_USER", jmsString);
                String jmsPassword = PropertyLoader.getSecretValue(PropertyLoader.jmsPasswordValue, PropertyLoader.jmsPasswordFile);
                jmsString.setString(jmsPassword);
                ncfile.write("JMS_PASSWORD", jmsString);
                jmsString.setString(VCellQueue.WorkerEventQueue.getName());
                ncfile.write("JMS_QUEUE", jmsString);
                jmsString.setString(VCellTopic.ServiceControlTopic.getName());
                ncfile.write("JMS_TOPIC", jmsString);
                jmsString.setString(simulation.getVersion().getOwner().getName());
                ncfile.write("VCELL_USER", jmsString);
                ArrayInt.D0 scalarJMS = new ArrayInt.D0();
                scalarJMS.set(Integer.parseInt(simulation.getVersion().getVersionKey() + ""));
                ncfile.write("SIMULATION_KEY", scalarJMS);
                scalarJMS.set(simTask.getSimulationJob().getJobIndex());
                ncfile.write("JOB_INDEX", scalarJMS);
            }
            ArrayDouble.D0 scalarDouble = new ArrayDouble.D0();
            // TStart, TEnd, SaveTime
            if ((timeBounds.getEndingTime() > timeBounds.getStartingTime()) && (outputTimeSpec.getOutputTimeStep() > 0)) {
                scalarDouble.set(timeBounds.getStartingTime());
                ncfile.write("TStart", scalarDouble);
                scalarDouble.set(timeBounds.getEndingTime());
                ncfile.write("TEnd", scalarDouble);
                scalarDouble.set(outputTimeSpec.getOutputTimeStep());
                ncfile.write("SaveTime", scalarDouble);
            } else {
                System.err.println("Time setting error. Ending time smaller than starting time or save interval is not a positive value.");
                throw new RuntimeException("Time setting error. Ending time smaller than starting time or save interval is not a positive value.");
            }
            // Volume
            // we set volume to 1. This model file cannot support multi-compartmental sizes.
            // When writting the rate constants, we must take the volume into account according to the reaction type.
            scalarDouble.set(1);
            ncfile.write("Volume", scalarDouble);
            // CellGrowthTime, CellGrowthTimeSD,
            scalarDouble.set(0);
            ncfile.write("CellGrowthTime", scalarDouble);
            ncfile.write("CellGrowthTimeSD", scalarDouble);
            // ExpType, Last Trial, Last Model, MaxNumModels, NumModels
            ArrayInt.D0 scalarInt = new ArrayInt.D0();
            scalarInt.set(0);
            ncfile.write("LastTrial", scalarInt);
            ncfile.write("LastModel", scalarInt);
            scalarInt.set(1);
            ncfile.write("ExpType", scalarInt);
            ncfile.write("MaxNumModels", scalarInt);
            ncfile.write("NumModels", scalarInt);
            // SpeciesSplitOnDivision
            ArrayInt A1 = new ArrayInt.D1(numSpecies.getLength());
            Index idx = A1.getIndex();
            for (int i = 0; i < numSpecies.getLength(); i++) {
                A1.setInt(idx.set(i), 0);
            }
            ncfile.write("SpeciesSplitOnDivision", new int[1], A1);
            // SaveSpeciesData
            ArrayInt A2 = new ArrayInt.D1(numSpecies.getLength());
            idx = A2.getIndex();
            for (int i = 0; i < numSpecies.getLength(); i++) {
                A2.setInt(idx.set(i), 1);
            }
            ncfile.write("SaveSpeciesData", new int[1], A2);
            // Reaction_Rate_Laws
            ArrayInt A3 = new ArrayInt.D1(numReactions.getLength());
            idx = A3.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                A3.setInt(idx.set(i), reactionRateLaws[i].getLawType());
            }
            ncfile.write("Reaction_Rate_Laws", new int[1], A3);
            // Reaction_DListLen
            ArrayInt A4 = new ArrayInt.D1(numReactions.getLength());
            idx = A4.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                if (reactionRateLaws[i].getLawType() == ReactionRateLaw.order_0)
                    A4.setInt(idx.set(i), 0);
                else if ((reactionRateLaws[i].getLawType() == ReactionRateLaw.order_1) || (reactionRateLaws[i].getLawType() == ReactionRateLaw.order_2_1substrate) || (reactionRateLaws[i].getLawType() == ReactionRateLaw.order_3_1substrate))
                    A4.setInt(idx.set(i), 1);
                else if ((reactionRateLaws[i].getLawType() == ReactionRateLaw.order_2_2substrate) || (reactionRateLaws[i].getLawType() == ReactionRateLaw.order_3_2substrate))
                    A4.setInt(idx.set(i), 2);
                else if (reactionRateLaws[i].getLawType() == ReactionRateLaw.order_3_3substrate)
                    A4.setInt(idx.set(i), 3);
            }
            ncfile.write("Reaction_DListLen", new int[1], A4);
            // Reaction_StoichListLen
            ArrayInt A5 = new ArrayInt.D1(numReactions.getLength());
            idx = A5.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                A5.setInt(idx.set(i), reactions[i].getActions().size());
            }
            ncfile.write("Reaction_StoichListLen", new int[1], A5);
            // Reaction_OptionalData
            ArrayInt A6 = new ArrayInt.D1(numReactions.getLength());
            idx = A6.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                A6.setInt(idx.set(i), 0);
            }
            ncfile.write("Reaction_OptionalData", new int[1], A6);
            // Reaction_StoichCoeff
            ArrayInt A7 = new ArrayInt.D2(numReactions.getLength(), numMaxStoichList.getLength());
            idx = A7.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                Action[] actions = (Action[]) reactions[i].getActions().toArray(new Action[reactions[i].getActions().size()]);
                for (int j = 0; j < actions.length; j++) {
                    try {
                        actions[j].getOperand().evaluateConstant();
                        int coeff = (int) Math.round(actions[j].getOperand().evaluateConstant());
                        A7.setInt(idx.set(i, j), coeff);
                    } catch (ExpressionException ex) {
                        ex.printStackTrace(System.err);
                        throw new ExpressionException(ex.getMessage());
                    }
                }
            }
            ncfile.write("Reaction_StoichCoeff", new int[2], A7);
            // Reaction_StoichSpecies
            ArrayInt A8 = new ArrayInt.D2(numReactions.getLength(), numMaxStoichList.getLength());
            idx = A8.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                ArrayList<Action> actions = reactions[i].getActions();
                for (int j = 0; j < actions.size(); j++) {
                    A8.setInt(idx.set(i, j), getVariableIndex(((Action) actions.get(j)).getVar().getName(), vars));
                }
            }
            ncfile.write("Reaction_StoichSpecies", new int[2], A8);
            // Reaction_DepList
            ArrayInt A9 = new ArrayInt.D2(numReactions.getLength(), numMaxDepList.getLength());
            idx = A9.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                ReactionRateLaw rl = reactionRateLaws[i];
                Hashtable<String, Integer> tem = varInProbOrderHash[i];
                Enumeration<String> varnames = tem.keys();
                if (rl.getLawType() == ReactionRateLaw.order_0) {
                // don't do anything here.
                } else if ((rl.getLawType() == ReactionRateLaw.order_1) || (rl.getLawType() == ReactionRateLaw.order_2_1substrate) || (rl.getLawType() == ReactionRateLaw.order_3_1substrate) || (rl.getLawType() == ReactionRateLaw.order_2_2substrate) || (rl.getLawType() == ReactionRateLaw.order_3_3substrate)) {
                    int j = 0;
                    while (varnames.hasMoreElements()) {
                        String name = varnames.nextElement();
                        A9.setInt(idx.set(i, j), getVariableIndex(name, vars));
                        j++;
                    }
                } else if (rl.getLawType() == ReactionRateLaw.order_3_2substrate) {
                    int order = 0;
                    String highOrderName = "";
                    String lowOrderName = "";
                    // we must make sure to put the higher order species first.
                    while (varnames.hasMoreElements()) {
                        lowOrderName = varnames.nextElement();
                        if (tem.get(lowOrderName) > order) {
                            String s = highOrderName;
                            highOrderName = lowOrderName;
                            lowOrderName = s;
                            order = tem.get(highOrderName);
                        }
                    }
                    A9.setInt(idx.set(i, 0), getVariableIndex(highOrderName, vars));
                    A9.setInt(idx.set(i, 1), getVariableIndex(lowOrderName, vars));
                }
            }
            ncfile.write("Reaction_DepList", new int[2], A9);
            // Reaction_names
            ArrayChar A10 = new ArrayChar.D2(numReactions.getLength(), stringLen.getLength());
            for (int i = 0; i < numReactions.getLength(); i++) {
                String name = reactions[i].getName();
                int diff = stringLen.getLength() - name.length();
                if (diff >= 0) {
                    for (int j = 0; j < diff; j++) {
                        name = name + " ";
                    }
                    A10.setString(i, name);
                } else
                    throw new RuntimeException("Name of Reaction:" + name + " is too long. Please shorten to " + stringLen.getLength() + " chars.");
            }
            ncfile.write("Reaction_names", A10);
            // Species_names
            ArrayChar A11 = new ArrayChar.D2(numSpecies.getLength(), stringLen.getLength());
            for (int i = 0; i < numSpecies.getLength(); i++) {
                String name = vars.elementAt(i).getName();
                int diff = stringLen.getLength() - name.length();
                if (diff >= 0) {
                    for (int j = 0; j < diff; j++) {
                        name = name + " ";
                    }
                    A11.setString(i, name);
                } else
                    throw new RuntimeException("Name of Species:" + name + " is too long. Please shorten to " + stringLen.getLength() + " chars.");
            }
            ncfile.write("Species_names", A11);
            // Species Initial Condition (in number of molecules).
            // Species iniCondition are sampled from a poisson distribution(which has a mean of the current iniExp value)
            RandomDataGenerator dist = new RandomDataGenerator();
            if (stochOpt.isUseCustomSeed()) {
                Integer randomSeed = stochOpt.getCustomSeed();
                if (randomSeed != null) {
                    dist.reSeed(randomSeed);
                }
            }
            ArrayLong A12 = new ArrayLong.D1(numSpecies.getLength());
            idx = A12.getIndex();
            for (int i = 0; i < numSpecies.getLength(); i++) {
                try {
                    VarIniCondition varIniCondition = subDomain.getVarIniCondition(vars.elementAt(i));
                    Expression varIniExp = varIniCondition.getIniVal();
                    varIniExp.bindExpression(simSymbolTable);
                    varIniExp = simSymbolTable.substituteFunctions(varIniExp).flatten();
                    double expectedCount = varIniExp.evaluateConstant();
                    long varCount = 0;
                    if (varIniCondition instanceof VarIniCount) {
                        varCount = (long) expectedCount;
                    } else {
                        if (expectedCount > 0) {
                            varCount = dist.nextPoisson(expectedCount);
                        }
                    }
                    A12.setLong(idx.set(i), varCount);
                } catch (ExpressionException ex) {
                    ex.printStackTrace(System.err);
                    throw new ExpressionException(ex.getMessage());
                }
            }
            ncfile.write("SpeciesIC", new int[1], A12);
            // Reaction_Rate_Constants(NumReactions, NumMaxDepList) ;
            ArrayDouble A13 = new ArrayDouble.D2(numReactions.getLength(), numMaxDepList.getLength());
            idx = A13.getIndex();
            for (int i = 0; i < numReactions.getLength(); i++) {
                ReactionRateLaw rl = reactionRateLaws[i];
                A13.setDouble(idx.set(i, 0), rl.getRateConstant());
            }
            ncfile.write("Reaction_Rate_Constants", A13);
        } catch (IOException ioe) {
            ioe.printStackTrace(System.err);
            throw new IOException("Error writing hybrid input file " + filename + ": " + ioe.getMessage());
        } catch (InvalidRangeException ire) {
            ire.printStackTrace(System.err);
            throw new InvalidRangeException("Error writing hybrid input file " + filename + ": " + ire.getMessage());
        }
        try {
            ncfile.close();
        } catch (IOException ioe) {
            throw new IOException("Error closing file " + filename + ". " + ioe.getMessage());
        }
    }
}
Also used : VarIniCondition(cbit.vcell.math.VarIniCondition) NonspatialStochSimOptions(cbit.vcell.solver.NonspatialStochSimOptions) ArrayList(java.util.ArrayList) Index(ucar.ma2.Index) ExpressionException(cbit.vcell.parser.ExpressionException) SubDomain(cbit.vcell.math.SubDomain) SolverTaskDescription(cbit.vcell.solver.SolverTaskDescription) Vector(java.util.Vector) UniformOutputTimeSpec(cbit.vcell.solver.UniformOutputTimeSpec) RandomDataGenerator(org.apache.commons.math3.random.RandomDataGenerator) InvalidRangeException(ucar.ma2.InvalidRangeException) SimulationSymbolTable(cbit.vcell.solver.SimulationSymbolTable) ArrayLong(ucar.ma2.ArrayLong) ArrayInt(ucar.ma2.ArrayInt) NetcdfFileWriteable(ucar.nc2.NetcdfFileWriteable) Action(cbit.vcell.math.Action) StochVolVariable(cbit.vcell.math.StochVolVariable) Variable(cbit.vcell.math.Variable) TimeBounds(cbit.vcell.solver.TimeBounds) ArrayChar(ucar.ma2.ArrayChar) ArrayDouble(ucar.ma2.ArrayDouble) JumpProcess(cbit.vcell.math.JumpProcess) StochVolVariable(cbit.vcell.math.StochVolVariable) VarIniCount(cbit.vcell.math.VarIniCount) Dimension(ucar.nc2.Dimension) IOException(java.io.IOException) Simulation(cbit.vcell.solver.Simulation) Expression(cbit.vcell.parser.Expression)

Aggregations

VarIniCondition (cbit.vcell.math.VarIniCondition)11 Expression (cbit.vcell.parser.Expression)9 ExpressionException (cbit.vcell.parser.ExpressionException)8 JumpProcess (cbit.vcell.math.JumpProcess)7 SubDomain (cbit.vcell.math.SubDomain)7 Action (cbit.vcell.math.Action)5 Function (cbit.vcell.math.Function)5 StochVolVariable (cbit.vcell.math.StochVolVariable)5 VarIniCount (cbit.vcell.math.VarIniCount)5 Variable (cbit.vcell.math.Variable)5 Vector (java.util.Vector)5 Constant (cbit.vcell.math.Constant)4 Equation (cbit.vcell.math.Equation)4 ModelException (cbit.vcell.model.ModelException)4 PropertyVetoException (java.beans.PropertyVetoException)4 MappingException (cbit.vcell.mapping.MappingException)3 CompartmentSubDomain (cbit.vcell.math.CompartmentSubDomain)3 MathException (cbit.vcell.math.MathException)3 MathFormatException (cbit.vcell.math.MathFormatException)3 ImageException (cbit.image.ImageException)2