Examples with StructureMapping cbit.vcell.mapping.StructureMapping used on opensource projects

Example 31 with StructureMapping

use of cbit.vcell.mapping.StructureMapping in project vcell by virtualcell.

the class SbmlVcmlConverter method main.

/**
 * @param args -import or -export
 */
public static void main(String[] args) {
    ResourceUtil.setNativeLibraryDirectory();
    if (args.length < 2 || args.length > 3) {
        System.out.println("Usage:\n\t -export path_of_input_file\n\tOR\n\t -import path_of_input_file [-simulate]");
        System.exit(1);
    }
    if (args.length > 2 && args[0].equals("-export")) {
        System.out.println("Export cannot have arguments other than input file.");
        System.out.println("Usage:\n\t -export path_of_input_file\n\tOR\n\t-import path_of_input_file [-simulate]");
        System.exit(1);
    }
    try {
        String pathName = args[1];
        // Read xml file (Sbml or Vcml) into stringbuffer, pass the string into VCell's importer/exporter
        String xmlString = getXMLString(pathName);
        if (args[0].equals("-import")) {
            if (args.length > 3 || (args.length == 3 && !args[2].equals("-simulate"))) {
                System.out.println(args[2] + " : Unknown arguments for import; please check and re-run import command");
                System.out.println("Usage:\n\t -import path_of_input_file [-simulate]");
                System.exit(1);
            }
            boolean bSimulate = false;
            if (args.length == 3 && args[2].equals("-simulate")) {
                bSimulate = true;
            }
            // Create a default VCLogger - SBMLImporter needs it
            cbit.util.xml.VCLogger logger = new cbit.util.xml.VCLogger() {

                @Override
                public void sendMessage(Priority p, ErrorType et, String message) {
                    System.err.println("LOGGER: msgLevel=" + p + ", msgType=" + et + ", " + message);
                    if (p == VCLogger.Priority.HighPriority) {
                        throw new RuntimeException("Import failed : " + message);
                    }
                }

                public void sendAllMessages() {
                }

                public boolean hasMessages() {
                    return false;
                }
            };
            // invoke SBMLImporter, which returns a Biomodel, convert that to VCML using XMLHelper
            try {
                // import SBML model into VCML, store VCML string in 'fileName.vcml'
                // Instantiate an SBMLImporter to get the speciesUnitsHash - to compute the conversion factor from VC->SB species units.
                // and import SBML  (sbml->bioModel)
                org.vcell.sbml.vcell.SBMLImporter sbmlImporter = new org.vcell.sbml.vcell.SBMLImporter(pathName, logger, false);
                BioModel bioModel = sbmlImporter.getBioModel();
                // Hashtable<String, SBMLImporter.SBVCConcentrationUnits> speciesUnitsHash = sbmlImporter.getSpeciesUnitsHash();
                // double timeFactor = sbmlImporter.getSBMLTimeUnitsFactor();
                // convert biomodel to vcml and save to file.
                String vcmlString = XmlHelper.bioModelToXML(bioModel);
                String fileExtensionStr;
                if (pathName.endsWith(".xml")) {
                    fileExtensionStr = ".xml";
                } else if (pathName.endsWith(".sbml")) {
                    fileExtensionStr = ".sbml";
                } else {
                    throw new RuntimeException("Unknown file extension for SBML file name; Exiting SbmlConverter.");
                }
                String vcmlFileName = pathName.replace(fileExtensionStr, ".vcml");
                File vcmlFile = new File(vcmlFileName);
                XmlUtil.writeXMLStringToFile(vcmlString, vcmlFile.getAbsolutePath(), true);
                // If user doesn't choose to simulate, you are done.
                if (!bSimulate) {
                    return;
                }
                // Generate math for lone simContext
                SimulationContext simContext = (SimulationContext) bioModel.getSimulationContext(0);
                MathDescription mathDesc = simContext.createNewMathMapping().getMathDescription();
                simContext.setMathDescription(mathDesc);
                // Create basic simulation, with IDA solver (set in solve method) and other defaults, and end time 'Te'
                org.vcell.util.document.SimulationVersion simVersion = new org.vcell.util.document.SimulationVersion(new KeyValue("12345"), "name", new org.vcell.util.document.User("user", new KeyValue("123")), new org.vcell.util.document.GroupAccessNone(), // versionBranchPointRef
                null, // branchID
                new java.math.BigDecimal(1.0), new java.util.Date(), org.vcell.util.document.VersionFlag.Archived, "", null);
                Simulation sim1 = new Simulation(simVersion, simContext.getMathDescription());
                simContext.addSimulation(sim1);
                sim1.setName("sim1");
                // using a default end time of 10.0 secs, and a forcing default output time step of 0.01.
                // If user wants anything different, user can modify the .idaInput file and re-run simulation separately.
                // double newEndTime = endTime * timeFactor;
                double newEndTime = endTime;
                sim1.getSolverTaskDescription().setTimeBounds(new cbit.vcell.solver.TimeBounds(0, newEndTime));
                TimeStep timeStep_1 = new TimeStep();
                sim1.getSolverTaskDescription().setTimeStep(new TimeStep(timeStep_1.getMinimumTimeStep(), timeStep_1.getDefaultTimeStep(), newEndTime / 10000));
                sim1.getSolverTaskDescription().setOutputTimeSpec(new UniformOutputTimeSpec((newEndTime - 0) / numTimeSteps));
                sim1.getSolverTaskDescription().setErrorTolerance(new ErrorTolerance(1e-10, 1e-12));
                // save the new vcml with generated math and new sim in file named : fileName_IDA_simulation.vcml
                String newVcmlString = XmlHelper.bioModelToXML(bioModel);
                String newVcmlFileName = vcmlFile.getPath().replace(".vcml", "_IDA_simulation.vcml");
                File newVcmlFile = new File(newVcmlFileName);
                XmlUtil.writeXMLStringToFile(newVcmlString, newVcmlFile.getAbsolutePath(), true);
                // Solve simulation, which also generates and saves the .idainput file and .csv
                // SimSpec to get vars to solve and output in .csv
                SimSpec simSpec = SimSpec.fromSBML(xmlString);
                solveSimulation(new SimulationJob(sim1, 0, null), vcmlFile.getPath(), simSpec);
            } catch (Exception e) {
                e.printStackTrace(System.err);
                // and an annotation saying why it failed.
                try {
                    String fileExtensionStr;
                    if (pathName.endsWith(".xml")) {
                        fileExtensionStr = ".xml";
                    } else if (pathName.endsWith(".sbml")) {
                        fileExtensionStr = ".sbml";
                    } else {
                        throw new RuntimeException("Unknown file extension for SBML file name; Exiting SbmlConverter.");
                    }
                    String dummyVcmlFileName = pathName.replace(fileExtensionStr, ".vcml");
                    File dummyVcmlFile = new File(dummyVcmlFileName);
                    BioModel dummy_biomodel = new BioModel(null);
                    String dummyName = dummyVcmlFile.getName().substring(0, dummyVcmlFile.getName().indexOf(".vcml"));
                    dummy_biomodel.setName(dummyName);
                    dummy_biomodel.setDescription("SBML Model could not be automatically converted to VCML : " + e.getMessage());
                    String vcmlString = XmlHelper.bioModelToXML(dummy_biomodel);
                    XmlUtil.writeXMLStringToFile(vcmlString, dummyVcmlFile.getAbsolutePath(), true);
                } catch (Exception e1) {
                    e.printStackTrace(System.err);
                }
            }
        } else if (args[0].equals("-export")) {
            try {
                BioModel bioModel = XmlHelper.XMLToBioModel(new XMLSource(xmlString));
                for (int i = 0; i < bioModel.getSimulationContexts().length; i++) {
                    SimulationContext simContext = bioModel.getSimulationContext(i);
                    if (simContext.getGeometry().getDimension() == 0 && !simContext.isStoch()) {
                        if (!simContext.getGeometryContext().isAllSizeSpecifiedPositive()) {
                            Structure structure = simContext.getModel().getStructure(0);
                            double structureSize = 1.0;
                            StructureMapping structMapping = simContext.getGeometryContext().getStructureMapping(structure);
                            StructureSizeSolver.updateAbsoluteStructureSizes(simContext, structure, structureSize, structMapping.getSizeParameter().getUnitDefinition());
                        }
                        // Export the application itself, with default overrides
                        String sbmlString = XmlHelper.exportSBML(bioModel, 2, 4, 0, false, simContext, null);
                        String filePath = pathName.substring(0, pathName.lastIndexOf("\\") + 1);
                        String sbmlFileName = TokenMangler.mangleToSName(bioModel.getName() + "_" + i);
                        File sbmlFile = new File(filePath + sbmlFileName + ".xml");
                        XmlUtil.writeXMLStringToFile(sbmlString, sbmlFile.getAbsolutePath(), true);
                        // Now export each simulation in the application that has overrides
                        Simulation[] simulations = bioModel.getSimulations(simContext);
                        for (int j = 0; j < simulations.length; j++) {
                            if (simulations[j].getMathOverrides().hasOverrides()) {
                                // Check for parameter scan and create simJob to pass into exporter
                                for (int k = 0; k < simulations[j].getScanCount(); k++) {
                                    SimulationJob simJob = new SimulationJob(simulations[j], k, null);
                                    sbmlString = XmlHelper.exportSBML(bioModel, 2, 4, 0, false, simContext, simJob);
                                    String fileName = TokenMangler.mangleToSName(sbmlFileName + "_" + j + "_" + k);
                                    sbmlFile = new File(filePath + fileName + ".xml");
                                    XmlUtil.writeXMLStringToFile(sbmlString, sbmlFile.getAbsolutePath(), true);
                                }
                            }
                        }
                    }
                }
                System.out.println("Successfully translated model : " + pathName);
            } catch (Exception e) {
                e.printStackTrace(System.err);
            }
        }
        System.exit(0);
    } catch (IOException e) {
        e.printStackTrace(System.err);
        System.exit(1);
    }
}

Example 32 with StructureMapping

use of cbit.vcell.mapping.StructureMapping in project vcell by virtualcell.

the class StructureMappingCartoonTool method menuAction.

protected void menuAction(Shape shape, String menuAction) {
    // 
    if (shape == null) {
        return;
    }
    // 
    if (menuAction.equals(CartoonToolEditActions.Delete.MENU_ACTION)) {
        if (shape instanceof StructureMappingShape) {
            try {
                StructureMapping sm = (StructureMapping) ((StructureMappingShape) shape).getModelObject();
                getStructureMappingCartoon().getGeometryContext().assignStructure(sm.getStructure(), null);
                getStructureMappingCartoon().refreshAll();
            } catch (IllegalMappingException e) {
                e.printStackTrace(System.out);
                PopupGenerator.showErrorDialog(getGraphPane(), e.getMessage());
            } catch (java.beans.PropertyVetoException e) {
                e.printStackTrace(System.out);
                PopupGenerator.showErrorDialog(getGraphPane(), e.getMessage());
            } catch (MappingException e) {
                e.printStackTrace(System.out);
                PopupGenerator.showErrorDialog(getGraphPane(), e.getMessage());
            }
        }
    } else {
        // default action is to ignore
        System.out.println("unsupported menu action '" + menuAction + "' on shape '" + shape + "'");
    }
}

Example 33 with StructureMapping

use of cbit.vcell.mapping.StructureMapping in project vcell by virtualcell.

the class SpeciesContextSpecParameterTableModel method propertyChange.

/**
 * This method gets called when a bound property is changed.
 * @param evt A PropertyChangeEvent object describing the event source
 *   and the property that has changed.
 */
public void propertyChange(java.beans.PropertyChangeEvent evt) {
    // 
    try {
        // 
        if (fieldSpeciesContextSpec != null && evt.getSource() == fieldSpeciesContextSpec.getSimulationContext().getGeometryContext() && evt.getPropertyName().equals("geometry")) {
            refreshData();
        }
        // 
        if (fieldSpeciesContextSpec != null && evt.getSource() == fieldSpeciesContextSpec.getSimulationContext().getGeometryContext() && evt.getPropertyName().equals("structureMappings")) {
            StructureMapping[] oldStructureMappings = (StructureMapping[]) evt.getOldValue();
            for (int i = 0; oldStructureMappings != null && i < oldStructureMappings.length; i++) {
                oldStructureMappings[i].removePropertyChangeListener(this);
            }
            StructureMapping[] newStructureMappings = (StructureMapping[]) evt.getNewValue();
            for (int i = 0; newStructureMappings != null && i < newStructureMappings.length; i++) {
                newStructureMappings[i].addPropertyChangeListener(this);
            }
            refreshData();
        }
        // 
        if (evt.getSource() instanceof StructureMapping) {
            refreshData();
        }
        if (evt.getSource() == this && evt.getPropertyName().equals("speciesContextSpec")) {
            SpeciesContextSpec oldValue = (SpeciesContextSpec) evt.getOldValue();
            if (oldValue != null) {
                oldValue.removePropertyChangeListener(this);
                Parameter[] oldParameters = oldValue.getParameters();
                if (oldParameters != null) {
                    for (int i = 0; i < oldParameters.length; i++) {
                        oldParameters[i].removePropertyChangeListener(this);
                    }
                }
                SimulationContext oldSimContext = oldValue.getSimulationContext();
                if (oldSimContext != null) {
                    oldSimContext.getGeometryContext().removePropertyChangeListener(this);
                    StructureMapping[] oldStructureMappings = oldSimContext.getGeometryContext().getStructureMappings();
                    if (oldStructureMappings != null) {
                        for (int i = 0; i < oldStructureMappings.length; i++) {
                            oldStructureMappings[i].removePropertyChangeListener(this);
                        }
                    }
                }
            }
            SpeciesContextSpec newValue = (SpeciesContextSpec) evt.getNewValue();
            if (newValue != null) {
                newValue.addPropertyChangeListener(this);
                Parameter[] newParameters = newValue.getParameters();
                if (newParameters != null) {
                    for (int i = 0; i < newParameters.length; i++) {
                        newParameters[i].addPropertyChangeListener(this);
                    }
                }
                SimulationContext newSimContext = newValue.getSimulationContext();
                if (newSimContext != null) {
                    newSimContext.getGeometryContext().addPropertyChangeListener(this);
                    StructureMapping[] newStructureMappings = newSimContext.getGeometryContext().getStructureMappings();
                    if (newStructureMappings != null) {
                        for (int i = 0; i < newStructureMappings.length; i++) {
                            newStructureMappings[i].addPropertyChangeListener(this);
                        }
                    }
                }
            }
            refreshData();
        }
        if (evt.getSource() instanceof SpeciesContextSpec) {
            // if parameters changed must update listeners
            if (evt.getPropertyName().equals("parameters")) {
                Parameter[] oldParameters = (Parameter[]) evt.getOldValue();
                if (oldParameters != null) {
                    for (int i = 0; i < oldParameters.length; i++) {
                        oldParameters[i].removePropertyChangeListener(this);
                    }
                }
                Parameter[] newParameters = (Parameter[]) evt.getNewValue();
                if (newParameters != null) {
                    for (int i = 0; i < newParameters.length; i++) {
                        newParameters[i].addPropertyChangeListener(this);
                    }
                }
            }
            if (!evt.getPropertyName().equals(SpeciesContextSpec.PARAMETER_NAME_PROXY_PARAMETERS)) {
                // for any change to the SpeciesContextSpec, want to update all.
                // proxy parameters don't affect table
                refreshData();
            }
        }
        if (evt.getSource() instanceof Parameter) {
            refreshData();
        }
    } catch (Exception e) {
        e.printStackTrace(System.out);
    }
}

Example 34 with StructureMapping

use of cbit.vcell.mapping.StructureMapping 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());
    }
}

Example 35 with StructureMapping

use of cbit.vcell.mapping.StructureMapping in project vcell by virtualcell.

the class StochMathMapping_4_8 method getProbabilityRate.

/**
 * Get probability expression for the specific elementary reaction.
 * Input: ReactionStep, the reaction. isForwardDirection, if the elementary reaction is forward from the reactionstep.
 * Output: Expression. the probability expression.
 * Creation date: (9/14/2006 3:22:58 PM)
 */
Expression getProbabilityRate(ReactionStep rs, boolean isForwardDirection) throws MappingException {
    ReactionStep reactionStep = rs;
    Expression probExp = null;
    // get kinetics of the reaction step
    Kinetics kinetics = reactionStep.getKinetics();
    // to compose the rate constant expression e.g. Kf, Kr
    Expression rateConstantExpr = null;
    // to compose the stochastic variable(species) expression, e.g. s*(s-1)*(s-2)* speciesFactor.
    Expression rxnProbabilityExpr = null;
    // to compose the factor that the probability expression multiplies with, which convert the rate expression under stochastic context
    Expression factorExpr = null;
    // the structure where reaction happens
    StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(rs.getStructure());
    Model model = getSimulationContext().getModel();
    try {
        if (// forward reaction
        isForwardDirection) {
            // for HMMs, it's a bit complicated. Vmax/(Km+s)-->Vmax*Size_s/(Km*Size_s+Ns)
            if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
                KineticsParameter kfp = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KForward);
                rateConstantExpr = new Expression(kfp, getNameScope());
            // rateConstantExpr.bindExpression(this);
            }
            // get convert factor for rate constant( membrane:rateConstant*membrane_Size (factor is membrane_size), feature : rateConstant*(feature_size/KMole)(factor is feature_size/KMOLE)) )
            if (sm.getStructure() instanceof Membrane) {
                factorExpr = new Expression(sm.getStructure().getStructureSize(), getNameScope());
            } else {
                factorExpr = new Expression(sm.getStructure().getStructureSize(), getNameScope());
                Expression kmoleExpr = new Expression(Model.reservedConstantsMap.get(ReservedSymbolRole.KMOLE));
                factorExpr = Expression.mult(factorExpr, kmoleExpr);
            }
            // complete the probability expression by the reactants' stoichiometries if it is Mass Action rate law
            if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
                ReactionParticipant[] reacPart = reactionStep.getReactionParticipants();
                for (int i = 0; i < reacPart.length; i++) {
                    int stoichiometry = 0;
                    if (reacPart[i] instanceof Reactant) {
                        stoichiometry = ((Reactant) reacPart[i]).getStoichiometry();
                        // ******the following part is to form the s*(s-1)(s-2)..(s-stoi+1).portion of the probability rate.
                        StructureMapping reactSM = getSimulationContext().getGeometryContext().getStructureMapping(reacPart[i].getStructure());
                        // factor expression for species
                        Expression speciesFactor = null;
                        // convert speceis' unit from moles/liter to molecules.
                        if (reactSM.getStructure() instanceof Membrane) {
                            speciesFactor = Expression.invert(new Expression(reactSM.getStructure().getStructureSize(), getNameScope()));
                        } else {
                            Expression exp1 = new Expression(Model.reservedConstantsMap.get(ReservedSymbolRole.KMOLE));
                            Expression exp2 = new Expression(reactSM.getStructure().getStructureSize(), getNameScope());
                            speciesFactor = Expression.div(Expression.invert(exp1), exp2);
                        }
                        // s*(s-1)(s-2)..(s-stoi+1)
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(reacPart[i].getSpeciesContext());
                        Expression spCount_exp = new Expression(spCountParam, getNameScope());
                        // species from uM to No. of Particles, form s*(s-1)*(s-2)
                        Expression tempExpr = new Expression(spCount_exp);
                        for (int j = 1; j < stoichiometry; j++) {
                            tempExpr = Expression.mult(tempExpr, Expression.add(spCount_exp, new Expression(-j)));
                        }
                        // update total factor with speceies factor
                        if (stoichiometry == 1) {
                            factorExpr = Expression.mult(factorExpr, speciesFactor);
                        } else if (stoichiometry > 1) {
                            // rxnProbExpr * (structSize^stoichiometry)
                            Expression powerExpr = Expression.power(speciesFactor, new Expression(stoichiometry));
                            factorExpr = Expression.mult(factorExpr, powerExpr);
                        }
                        if (rxnProbabilityExpr == null) {
                            rxnProbabilityExpr = new Expression(tempExpr);
                        } else {
                            // for more than one reactant
                            rxnProbabilityExpr = Expression.mult(rxnProbabilityExpr, tempExpr);
                        }
                    }
                }
            }
        } else // reverse reaction
        {
            if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
                KineticsParameter krp = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KReverse);
                rateConstantExpr = new Expression(krp, getNameScope());
            // rateConstantExpr.bindExpression(this);
            }
            // get convert factor for rate constant( membrane:rateConstant*membrane_Size (factor is membrane_size), feature : rateConstant*(feature_size/KMole)(factor is feature_size/KMOLE)) )
            if (sm.getStructure() instanceof Membrane) {
                factorExpr = new Expression(sm.getStructure().getStructureSize(), getNameScope());
            } else {
                factorExpr = new Expression(sm.getStructure().getStructureSize(), getNameScope());
                Expression exp = new Expression(Model.reservedConstantsMap.get(ReservedSymbolRole.KMOLE));
                factorExpr = Expression.mult(factorExpr, exp);
            }
            // complete the remaining part of the probability expression by the products' stoichiometries.
            if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
                ReactionParticipant[] reacPart = reactionStep.getReactionParticipants();
                for (int i = 0; i < reacPart.length; i++) {
                    int stoichiometry = 0;
                    if (reacPart[i] instanceof Product) {
                        stoichiometry = ((Product) reacPart[i]).getStoichiometry();
                        // ******the following part is to form the s*(s-1)*(s-2)...(s-stoi+1).portion of the probability rate.
                        StructureMapping reactSM = getSimulationContext().getGeometryContext().getStructureMapping(reacPart[i].getStructure());
                        // factor expression for species
                        Expression speciesFactor = null;
                        // convert speceis' unit from moles/liter to molecules.
                        if (reactSM.getStructure() instanceof Membrane) {
                            speciesFactor = Expression.invert(new Expression(reactSM.getStructure().getStructureSize(), getNameScope()));
                        } else {
                            Expression exp1 = new Expression(Model.reservedConstantsMap.get(ReservedSymbolRole.KMOLE));
                            Expression exp2 = new Expression(reactSM.getStructure().getStructureSize(), getNameScope());
                            speciesFactor = Expression.div(Expression.invert(exp1), exp2);
                        }
                        // s*(s-1)*(s-2)...(s-stoi+1)
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(reacPart[i].getSpeciesContext());
                        Expression spCount_exp = new Expression(spCountParam, getNameScope());
                        // species from uM to No. of Particles, form s*(s-1)*(s-2)
                        Expression tempExpr = new Expression(spCount_exp);
                        for (int j = 1; j < stoichiometry; j++) {
                            tempExpr = Expression.mult(tempExpr, Expression.add(spCount_exp, new Expression(-j)));
                        }
                        // update total factor with speceies factor
                        if (stoichiometry == 1) {
                            factorExpr = Expression.mult(factorExpr, speciesFactor);
                        } else if (stoichiometry > 1) {
                            // rxnProbExpr * (structSize^stoichiometry)
                            Expression powerExpr = Expression.power(speciesFactor, new Expression(stoichiometry));
                            factorExpr = Expression.mult(factorExpr, powerExpr);
                        }
                        if (rxnProbabilityExpr == null) {
                            rxnProbabilityExpr = new Expression(tempExpr);
                        } else {
                            rxnProbabilityExpr = Expression.mult(rxnProbabilityExpr, tempExpr);
                        }
                    }
                }
            }
        }
        // Now construct the probability expression.
        if (rateConstantExpr == null) {
            throw new MappingException("Can not find reaction rate constant in reaction: " + reactionStep.getName());
        } else if (rxnProbabilityExpr == null) {
            probExp = new Expression(rateConstantExpr);
        } else if ((rateConstantExpr != null) && (rxnProbabilityExpr != null)) {
            probExp = Expression.mult(rateConstantExpr, rxnProbabilityExpr);
        }
        // simplify the factor
        RationalExp factorRatExp = RationalExpUtils.getRationalExp(factorExpr);
        factorExpr = new Expression(factorRatExp.infixString());
        factorExpr.bindExpression(this);
        // get probability rate with converting factor
        probExp = Expression.mult(probExp, factorExpr);
        probExp = probExp.flatten();
    // //
    // // round trip to rational expression for simplifying terms like KMOLE/KMOLE ...
    // // we don't want to loose the symbol binding ... so we make a temporary symbolTable from the original binding.
    // //
    // final Expression finalExp = new Expression(probExp);
    // SymbolTable symbolTable = new SymbolTable(){
    // public void getEntries(Map<String, SymbolTableEntry> entryMap) {
    // throw new RuntimeException("should not be called");
    // }
    // public SymbolTableEntry getEntry(String identifierString) throws ExpressionBindingException {
    // return finalExp.getSymbolBinding(identifierString);
    // }
    // };
    // cbit.vcell.matrix.RationalExp ratExp = cbit.vcell.parser.RationalExpUtils.getRationalExp(probExp);
    // probExp = new Expression(ratExp.infixString());
    // probExp.bindExpression(symbolTable);
    } catch (ExpressionException e) {
        e.printStackTrace();
    }
    return probExp;
}