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

Example 1 with SubDomain

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

the class SmoldynFileWriter method writeMolecules.

private void writeMolecules() throws ExpressionException, MathException {
    // write molecules
    StringBuilder sb = new StringBuilder();
    int max_mol = 0;
    Enumeration<SubDomain> subDomainEnumeration = mathDesc.getSubDomains();
    while (subDomainEnumeration.hasMoreElements()) {
        SubDomain subDomain = subDomainEnumeration.nextElement();
        for (ParticleProperties particleProperties : subDomain.getParticleProperties()) {
            ArrayList<ParticleInitialCondition> particleInitialConditions = particleProperties.getParticleInitialConditions();
            String variableName = getVariableName(particleProperties.getVariable(), subDomain);
            for (ParticleInitialCondition pic : particleInitialConditions) {
                if (pic instanceof ParticleInitialConditionCount) {
                    max_mol += writeInitialCount((ParticleInitialConditionCount) pic, subDomain, variableName, sb);
                } else if (pic instanceof ParticleInitialConditionConcentration) {
                    max_mol += writeInitialConcentration((ParticleInitialConditionConcentration) pic, subDomain, particleProperties.getVariable(), variableName, sb);
                }
            }
            if (lg.isDebugEnabled()) {
                lg.debug("subdomain " + subDomain.getName() + ' ' + variableName + " processed, maximum mol estimate now " + max_mol);
            }
        }
    }
    if (max_mol > MAX_MOLECULE_LIMIT) {
        throw new MathException(VCellErrorMessages.getSmoldynMaxMolReachedErrorMessage((long) max_mol, MAX_MOLECULE_LIMIT));
    }
    int max_adjusted = max_mol * MOLECULE_MAX_COEFFICIENT;
    if (max_adjusted < MIN_MOLECULE_LIMIT) {
        if (lg.isInfoEnabled()) {
            lg.info("adjusting computed max " + max_adjusted + " to minimum " + MIN_MOLECULE_LIMIT);
        }
        max_adjusted = MIN_MOLECULE_LIMIT;
    }
    if (max_adjusted > MAX_MOLECULE_LIMIT) {
        if (lg.isInfoEnabled()) {
            lg.info("adjusting computed max " + max_adjusted + " to maximum " + MAX_MOLECULE_LIMIT);
        }
        max_adjusted = MAX_MOLECULE_LIMIT;
    }
    printWriter.println("# molecules");
    printWriter.println(SmoldynVCellMapper.SmoldynKeyword.max_mol + " " + max_adjusted);
    printWriter.println(sb);
}

Example 2 with SubDomain

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

the class SmoldynFileWriter method writeReactions.

private void writeReactions() throws ExpressionException, MathException {
    printWriter.println("# reactions");
    Enumeration<SubDomain> subdomains = mathDesc.getSubDomains();
    while (subdomains.hasMoreElements()) {
        SubDomain subdomain = subdomains.nextElement();
        for (ParticleJumpProcess pjp : subdomain.getParticleJumpProcesses()) {
            ArrayList<Variable> reactants = new ArrayList<Variable>();
            ArrayList<Variable> products = new ArrayList<Variable>();
            for (Action a : pjp.getActions().toArray(new Action[pjp.getActions().size()])) {
                if (a.getOperation().equals(Action.ACTION_CREATE)) {
                    products.add(a.getVar());
                } else if (a.getOperation().equals(Action.ACTION_DESTROY)) {
                    reactants.add(a.getVar());
                }
            }
            Expression rateDefinition = null;
            JumpProcessRateDefinition jprd = pjp.getParticleRateDefinition();
            if (jprd instanceof MacroscopicRateConstant) {
                rateDefinition = subsituteFlatten(((MacroscopicRateConstant) jprd).getExpression());
            } else if (jprd instanceof InteractionRadius) {
                rateDefinition = subsituteFlatten(((InteractionRadius) jprd).getExpression());
            } else {
                new RuntimeException("The jump process rate definition is not supported");
            }
            if (rateDefinition.isZero()) {
                continue;
            }
            if (mathDesc.isSpatialHybrid()) {
                String[] symbols = rateDefinition.getSymbols();
                if (symbols != null) {
                    if (subdomain instanceof MembraneSubDomain) {
                        rateDefinition = new Expression(FiniteVolumeFileWriter.replaceVolumeVariable(getSimulationTask(), (MembraneSubDomain) subdomain, rateDefinition));
                    }
                }
            } else {
                try {
                    rateDefinition.evaluateConstant();
                } catch (ExpressionException ex) {
                    throw new ExpressionException("reaction rate for jump process " + pjp.getName() + " is not a constant. Constants are required for all reaction rates.");
                }
            }
            // Smoldyn takes maximum 2nd order reaction.
            if (reactants.size() > 2) {
                throw new MathException("VCell spatial stochastic models support up to 2nd order reactions. \n" + "The reaction:" + pjp.getName() + " has more than 2 reactants.");
            }
            if (products.size() > 2) {
                throw new MathException("VCell spatial stochastic models support up to 2nd order reactions. \n" + "The reaction:" + pjp.getName() + " has more than 2 products.");
            }
            String rateDefinitionStr = simulation.getMathDescription().isSpatialHybrid() ? rateDefinition.infix() + ";" : rateDefinition.evaluateConstant() + "";
            if (subdomain instanceof CompartmentSubDomain) {
                // 0th order reaction, product limited to one and we'll let the reaction know where it happens
                if (reactants.size() == 0 && products.size() == 1) {
                    printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_cmpt + " " + subdomain.getName() + " " + pjp.getName() + " ");
                } else {
                    printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction + " " + /* + subdomain.getName() + " "*/
                    pjp.getName() + " ");
                }
                writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
            } else if (subdomain instanceof MembraneSubDomain) {
                // 0th order reaction, product limited to one and it can be on mem or in vol
                if (reactants.size() == 0 && products.size() == 1) {
                    printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
                    writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
                } else // consuming of a species to nothing, limited to one reactant
                if (reactants.size() == 1 && products.size() == 0) {
                    if (// consuming a mem species in mem reaction
                    getMembraneVariableCount(reactants) == 1) {
                        printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
                        writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
                    } else // it equals to adsorption, species A from vol adsorbed to mem as again species A, and then we kill the speceis A on mem.
                    if (getVolumeVariableCount(reactants) == 1) {
                        writeRateTransitionCommand(reactants, products, subdomain, rateDefinitionStr);
                        String speciesName = reactants.get(0).getName();
                        String killMolCmd = "cmd " + SmoldynVCellMapper.SmoldynKeyword.E + " " + SmoldynVCellMapper.SmoldynKeyword.killmol + " " + speciesName + "(" + SmoldynKeyword.up + ")";
                        killMolCommands.add(killMolCmd);
                    }
                } else // Use rate command for any membrane reactions with 1 reactant and 1 product
                if ((reactants.size() == 1) && (products.size() == 1)) {
                    // Membrane reaction (1 react to 1 product).
                    if (getMembraneVariableCount(products) == 1 && getMembraneVariableCount(reactants) == 1) {
                        printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
                        writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
                    } else // Other single molecular reactions
                    {
                        writeRateTransitionCommand(reactants, products, subdomain, rateDefinitionStr);
                    }
                } else // membrane reactions which are not one to one, or 0th order, or consuming species
                {
                    if (// membrane reaction has one membrane bound reactant
                    (getMembraneVariableCount(reactants) == 1)) {
                        printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
                        writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
                    } else if (// bimolecular membrane reaction
                    getMembraneVariableCount(reactants) == 2) {
                        if (jprd instanceof InteractionRadius) {
                            printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
                            writeReactionByInteractionRadius(reactants, products, subdomain, rateDefinitionStr, pjp.getName());
                        } else {
                            // throw new MathException("Error with reaction: " + pjp.getName() + ".\nVCell Spatial stochastic modeling requires macroscopic or microscopic kinetics for bimolecular membrane reactions.");
                            printWriter.print(SmoldynVCellMapper.SmoldynKeyword.reaction_surface + " " + subdomain.getName() + " " + pjp.getName() + " ");
                            writeReactionCommand(reactants, products, subdomain, rateDefinitionStr);
                        }
                    } else if (getMembraneVariableCount(reactants) == 0) {
                        throw new MathException("Error with reaction: " + pjp.getName() + ".\nIn VCell spatial stochastic modeling, the membrane reaction requires at least one membrane bound reactant.");
                    }
                }
            }
        }
    }
    printWriter.println();
}

Example 3 with SubDomain

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

the class SmoldynFileWriter method writeDrifts.

private void writeDrifts() throws ExpressionBindingException, ExpressionException, MathException {
    // writer diffusion properties
    printWriter.println("# drift properties");
    Enumeration<SubDomain> subDomainEnumeration = mathDesc.getSubDomains();
    while (subDomainEnumeration.hasMoreElements()) {
        SubDomain subDomain = subDomainEnumeration.nextElement();
        List<ParticleProperties> particlePropertiesList = subDomain.getParticleProperties();
        for (ParticleProperties pp : particlePropertiesList) {
            String variableName = null;
            if (subDomain instanceof MembraneSubDomain) {
                variableName = getVariableName(pp.getVariable(), subDomain);
            } else {
                variableName = getVariableName(pp.getVariable(), null);
            }
            try {
                double driftX = 0.0;
                if (pp.getDriftX() != null) {
                    driftX = subsituteFlattenToConstant(pp.getDriftX());
                }
                double driftY = 0.0;
                if (pp.getDriftY() != null) {
                    driftY = subsituteFlattenToConstant(pp.getDriftY());
                }
                double driftZ = 0.0;
                if (pp.getDriftZ() != null) {
                    driftZ = subsituteFlattenToConstant(pp.getDriftZ());
                }
                printWriter.println(SmoldynVCellMapper.SmoldynKeyword.drift + " " + variableName + " " + driftX + " " + driftY + " " + driftZ);
            } catch (NotAConstantException ex) {
                throw new ExpressionException("diffusion coefficient for variable " + variableName + " is not a constant. Constants are required for all diffusion coefficients");
            }
        }
    }
    printWriter.println();
}

Example 4 with SubDomain

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

the class Xmlproducer method getXML.

/**
 * This method returns a XML representation of a MathDescription object.
 * Creation date: (3/2/2001 10:57:25 AM)
 * @return Element
 * @param mathdes cbit.vcell.math.MathDescription
 */
Element getXML(MathDescription mathdes) throws XmlParseException {
    Element math = new Element(XMLTags.MathDescriptionTag);
    // Add attributes
    math.setAttribute(XMLTags.NameAttrTag, mangle(mathdes.getName()));
    // Add annotation
    if (mathdes.getDescription() != null && mathdes.getDescription().length() > 0) {
        Element annotationElem = new Element(XMLTags.AnnotationTag);
        annotationElem.setText(mangle(mathdes.getDescription()));
        math.addContent(annotationElem);
    }
    List<ParticleMolecularType> particleMolecularTypes = mathdes.getParticleMolecularTypes();
    for (ParticleMolecularType particleMolecularType : particleMolecularTypes) {
        math.addContent(getXML(particleMolecularType));
    }
    // Add Constant subelements
    Enumeration<Variable> enum1 = mathdes.getVariables();
    /*java.util.Iterator k;
    try {
    	VariableHash varHash = new VariableHash();
    	while (enum1.hasMoreElements()) 
    	 	varHash.addVariable((Variable)enum1.nextElement());
    	Variable vars [] = varHash.getReorderedVariables();
    	k = new ArrayList(java.util.Arrays.asList(vars)).iterator();
    } catch (cbit.vcell.mapping.MappingException e) {
		e.printStackTrace();
		return null;
    }*/
    while (enum1.hasMoreElements()) {
        Variable var = enum1.nextElement();
        Element element = null;
        if (var instanceof Constant) {
            element = getXML((Constant) var);
        } else if (var instanceof FilamentRegionVariable) {
            element = getXML((FilamentRegionVariable) var);
        } else if (var instanceof FilamentVariable) {
            element = getXML((FilamentVariable) var);
        } else if (var instanceof PointVariable) {
            element = getXML((PointVariable) var);
        } else if (var instanceof Function) {
            element = getXML((Function) var);
        } else if (var instanceof RandomVariable) {
            element = getXML((RandomVariable) var);
        } else if (var instanceof InsideVariable) {
            // *** for internal use! Ignore it ***
            continue;
        } else if (var instanceof MembraneRegionVariable) {
            element = getXML((MembraneRegionVariable) var);
        } else if (var instanceof MemVariable) {
            element = getXML((MemVariable) var);
        } else if (var instanceof OutsideVariable) {
            // *** for internal use! Ignore it ****
            continue;
        } else if (var instanceof VolumeRegionVariable) {
            element = getXML((VolumeRegionVariable) var);
        } else if (var instanceof VolVariable) {
            element = getXML((VolVariable) var);
        } else if (var instanceof StochVolVariable) {
            // added for stochastic volumn variables
            element = getXML((StochVolVariable) var);
        } else if (var instanceof ParticleVariable) {
            element = getXML((ParticleVariable) var);
        } else if (var instanceof ParticleObservable) {
            element = getXML((ParticleObservable) var);
        } else {
            throw new XmlParseException("An unknown variable type " + var.getClass().getName() + " was found when parsing the mathdescription " + mathdes.getName() + "!");
        }
        transcribeComments(var, element);
        math.addContent(element);
    }
    // this was moved to the simspec!
    /*	buffer.append("\n");
    	if (geometry != null){
    		buffer.append(geometry.getXML());
    	}	
    	buffer.append("\n");*/
    // Add subdomains
    Enumeration<SubDomain> enum2 = mathdes.getSubDomains();
    while (enum2.hasMoreElements()) {
        SubDomain subDomain = enum2.nextElement();
        math.addContent(getXML(subDomain));
    }
    // Add Metadata (Version) if there is!
    if (mathdes.getVersion() != null) {
        math.addContent(getXML(mathdes.getVersion(), mathdes));
    }
    Iterator<Event> iter = mathdes.getEvents();
    while (iter.hasNext()) {
        math.addContent(getXML(iter.next()));
    }
    PostProcessingBlock postProcessingBlock = mathdes.getPostProcessingBlock();
    if (postProcessingBlock.getNumDataGenerators() > 0) {
        math.addContent(getXML(postProcessingBlock));
    }
    return math;
}

Example 5 with SubDomain

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

the class StochMathMapping method refreshMathDescription.

/**
 * set up a math description based on current simulationContext.
 */
@Override
protected void refreshMathDescription() throws MappingException, MatrixException, MathException, ExpressionException, ModelException {
    // use local variable instead of using getter all the time.
    SimulationContext simContext = getSimulationContext();
    GeometryClass geometryClass = simContext.getGeometry().getGeometrySpec().getSubVolumes()[0];
    Domain domain = new Domain(geometryClass);
    // 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);
    }
    simContext.checkValidity();
    // 
    if (simContext.getGeometry().getDimension() > 0) {
        throw new MappingException("nonspatial stochastic math mapping requires 0-dimensional geometry");
    }
    // 
    for (int i = 0; i < structureMappings.length; i++) {
        if (structureMappings[i] instanceof MembraneMapping) {
            if (((MembraneMapping) structureMappings[i]).getCalculateVoltage()) {
                throw new MappingException("electric potential not yet supported for particle models");
            }
        }
    }
    // 
    // fail if any events
    // 
    BioEvent[] bioEvents = simContext.getBioEvents();
    if (bioEvents != null && bioEvents.length > 0) {
        throw new MappingException("events not yet supported for particle-based models");
    }
    // 
    // 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 && ((FeatureMapping) sm).getGeometryClass() == 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++) {
        Structure[] mappedStructures = simContext.getGeometryContext().getStructuresFromGeometryClass(subVolumes[i]);
        if (mappedStructures == null || mappedStructures.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()) {
            rsList.add(reactionSpecs[i].getReactionStep());
        }
    }
    // 
    for (ReactionStep reactionStep : rsList) {
        Kinetics.UnresolvedParameter[] unresolvedParameters = reactionStep.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("In Application '" + simContext.getName() + "', " + reactionStep.getDisplayType() + " '" + reactionStep.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();
    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.getPI_CONSTANT(), null), getIdentifierSubstitutions(model.getPI_CONSTANT().getExpression(), model.getPI_CONSTANT().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());
        }
    }
    // 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(), geometryClass);
        varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], geometryClass), expr, geometryClass));
    }
    // 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.getGeometryClass()), getIdentifierSubstitutions(memMapping.getInitialVoltageParameter().getExpression(), memMapping.getInitialVoltageParameter().getUnitDefinition(), memMapping.getGeometryClass()), memMapping.getGeometryClass()));
            } catch (ExpressionException e) {
                e.printStackTrace(System.out);
                throw new MappingException("Membrane initial voltage: " + initialVoltageParm.getName() + " cannot be evaluated as constant.");
            }
        }
    }
    // 
    for (ReactionStep rs : rsList) {
        if (rs.getKinetics() instanceof LumpedKinetics) {
            throw new RuntimeException("Lumped Kinetics not yet supported for Stochastic Math Generation");
        }
        Kinetics.KineticsParameter[] parameters = rs.getKinetics().getKineticsParameters();
        for (KineticsParameter parameter : parameters) {
            // 
            if ((parameter.getRole() == Kinetics.ROLE_CurrentDensity) && (parameter.getExpression() == null || parameter.getExpression().isZero())) {
                continue;
            }
            // 
            // don't add rate, we'll do it later when creating the jump processes
            // 
            // if (parameter.getRole() == Kinetics.ROLE_ReactionRate) {
            // continue;
            // }
            // 
            // don't add mass action reverse parameter if irreversible
            // 
            // if (!rs.isReversible() && parameters[i].getRole() == Kinetics.ROLE_KReverse){
            // continue;
            // }
            Expression expr = getSubstitutedExpr(parameter.getExpression(), true, false);
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameter, geometryClass), getIdentifierSubstitutions(expr, parameter.getUnitDefinition(), geometryClass), geometryClass));
        }
    }
    // 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.getGeometryClass()), 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.");
            }
        }
    }
    SpeciesContextSpec[] speciesContextSpecs = getSimulationContext().getReactionContext().getSpeciesContextSpecs();
    addInitialConditions(domain, speciesContextSpecs, varHash);
    // 
    // 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 in Application " + simContext.getName());
    }
    // 
    // create subDomains
    // 
    SubVolume subVolume = simContext.getGeometry().getGeometrySpec().getSubVolumes()[0];
    SubDomain subDomain = new CompartmentSubDomain(subVolume.getName(), 0);
    mathDesc.addSubDomain(subDomain);
    // 
    // 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.getGeometryClass()), getIdentifierSubstitutions(exp, spCountParam.getUnitDefinition(), sm.getGeometryClass()), domain));
        }
    }
    addJumpProcesses(varHash, geometryClass, subDomain);
    // 
    for (int i = 0; i < fieldMathMappingParameters.length; i++) {
        if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass()));
        }
    }
    // 
    // 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.getGeometryClass());
        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 = null;
            if (!scSpecs[i].isConstant() && getSimulationContext().isRandomizeInitCondition()) {
                varIni = new VarIniPoissonExpectedCount(var, new Expression(getMathSymbol(initParm, sm.getGeometryClass())));
            } else {
                varIni = new VarIniCount(var, new Expression(getMathSymbol(initParm, sm.getGeometryClass())));
            }
            subDomain.addVarIniCondition(varIni);
        }
    }
    // 
    for (int i = 0; i < fieldMathMappingParameters.length; i++) {
        if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
            Variable variable = newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass());
            if (mathDesc.getVariable(variable.getName()) == null) {
                mathDesc.addVariable(variable);
            }
        }
        if (fieldMathMappingParameters[i] instanceof ObservableCountParameter) {
            Variable variable = newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass());
            if (mathDesc.getVariable(variable.getName()) == null) {
                mathDesc.addVariable(variable);
            }
        }
    }
    if (!mathDesc.isValid()) {
        System.out.println(mathDesc.getVCML_database());
        throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
    }
}