Search in sources :

Example 51 with Structure

use of cbit.vcell.model.Structure in project vcell by virtualcell.

the class ReactionCartoonFull method applyDefaults.

public void applyDefaults(Diagram diagram) {
    List<NodeReference> nodeList = diagram.getNodeFullList();
    List<NodeReference> orphansList = new ArrayList<NodeReference>();
    for (int i = 0; i < nodeList.size(); i++) {
        NodeReference node = nodeList.get(i);
        Object obj = null;
        Structure struct = diagram.getStructure();
        boolean found = false;
        switch(node.nodeType) {
            case NodeReference.SIMPLE_REACTION_NODE:
                obj = getModel().getReactionStep(node.name);
                if (!(obj instanceof SimpleReaction)) {
                    System.out.println("ReactionCartoon.applyDefaults(), diagram reaction " + node.name + " type mismatch in model, using location anyway");
                }
                break;
            case NodeReference.FLUX_REACTION_NODE:
                obj = getModel().getReactionStep(node.name);
                if (!(obj instanceof FluxReaction)) {
                    System.out.println("ReactionCartoon.applyDefaults(), diagram flux " + node.name + " type mismatch in model, using location anyway");
                }
                break;
            case NodeReference.SPECIES_CONTEXT_NODE:
                obj = getModel().getSpeciesContext(node.name);
                break;
            case NodeReference.REACTION_RULE_NODE:
                obj = getModel().getRbmModelContainer().getReactionRule(node.name);
                break;
            case // obj is a RuleParticipantSignature
            NodeReference.RULE_PARTICIPANT_SIGNATURE_FULL_NODE:
                for (RuleParticipantSignature signature : ruleParticipantSignatures) {
                    if (signature instanceof RuleParticipantLongSignature && signature.getStructure() == struct && signature.compareByCriteria(node.getName(), GroupingCriteria.full)) {
                        obj = signature;
                        found = true;
                        break;
                    }
                }
                if (!found) {
                    orphansList.add(node);
                }
                break;
            case NodeReference.RULE_PARTICIPANT_SIGNATURE_SHORT_NODE:
                System.out.println("ReactionCartoonFull, RULE_PARTICIPANT_SIGNATURE_SHORT_NODE detected");
                for (RuleParticipantSignature signature : ruleParticipantSignatures) {
                    if (signature instanceof RuleParticipantShortSignature && signature.getStructure() == struct && signature.compareByCriteria(node.getName(), GroupingCriteria.full)) {
                        obj = signature;
                        found = true;
                        break;
                    }
                }
                if (!found) {
                    orphansList.add(node);
                }
                break;
        }
        // -- switch
        Shape shape = getShapeFromModelObject(obj);
        if (shape != null) {
            Point relPosOld = shape.getRelPos();
            Point relPosNew = node.location;
            // In old models, the same node can appear in multiple diagrams.
            // Now, we have only one diagram, so if a node has multiple positions,
            // some would overwrite others.
            // This attempts to prevent overwriting a position with a worse one.
            // if(relPosOld.x + relPosOld.y < relPosNew.x + relPosNew.y) {
            shape.setRelPos(relPosNew);
        // }
        }
    }
    if (!orphansList.isEmpty()) {
        diagram.removeNodeReferences(NodeReference.Mode.full, orphansList);
    }
}
Also used : RuleParticipantSignature(cbit.vcell.model.RuleParticipantSignature) RuleParticipantLongSignature(cbit.vcell.model.RuleParticipantLongSignature) SimpleReaction(cbit.vcell.model.SimpleReaction) Shape(cbit.gui.graph.Shape) NodeReference(cbit.vcell.model.NodeReference) ArrayList(java.util.ArrayList) FluxReaction(cbit.vcell.model.FluxReaction) Point(java.awt.Point) Point(java.awt.Point) RuleParticipantShortSignature(cbit.vcell.model.RuleParticipantShortSignature) Structure(cbit.vcell.model.Structure)

Example 52 with Structure

use of cbit.vcell.model.Structure in project vcell by virtualcell.

the class DiffEquMathMapping method refreshMathDescription.

/**
 * This method was created in VisualAge.
 */
@SuppressWarnings("deprecation")
private void refreshMathDescription() throws MappingException, MatrixException, MathException, ExpressionException, ModelException {
    // All sizes must be set for new ODE models and ratios must be set for old ones.
    simContext.checkValidity();
    // 
    // temporarily place all variables in a hashtable (before binding) and discarding duplicates (check for equality)
    // 
    VariableHash varHash = new VariableHash();
    StructureMapping[] structureMappings = simContext.getGeometryContext().getStructureMappings();
    // 
    // 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.getGeometryClass() == null)){
    // localIssueList.add(new Issue(structures[i], IssueCategory.StructureNotMapped,"In Application '" + simContext.getName() + "', model structure '"+structures[i].getName()+"' not mapped to a geometry subdomain",Issue.SEVERITY_WARNING));
    // }
    // }
    // 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){
    // localIssueList.add(new Issue(subVolumes[i], IssueCategory.GeometryClassNotMapped,"In Application '" + simContext.getName() + "', geometry subVolume '"+subVolumes[i].getName()+"' not mapped from a model structure",Issue.SEVERITY_WARNING));
    // }
    // }
    // deals with model parameters
    HashMap<VolVariable, EventAssignmentOrRateRuleInitParameter> eventOrRateRuleVolVarHash = new HashMap<VolVariable, EventAssignmentOrRateRuleInitParameter>();
    HashMap<VolVariable, RateRuleRateParameter> rateRuleRateParamHash = new HashMap<VolVariable, RateRuleRateParameter>();
    ArrayList<SymbolTableEntry> rateRuleVarTargets = new ArrayList<SymbolTableEntry>();
    Model model = simContext.getModel();
    ModelUnitSystem modelUnitSystem = model.getUnitSystem();
    VCUnitDefinition timeUnit = modelUnitSystem.getTimeUnit();
    ModelParameter[] modelParameters = model.getModelParameters();
    if (simContext.getGeometry().getDimension() == 0) {
        // 
        // global parameters from model (that presently are constants)
        // 
        BioEvent[] bioEvents = simContext.getBioEvents();
        ArrayList<SymbolTableEntry> eventAssignTargets = new ArrayList<SymbolTableEntry>();
        if (bioEvents != null && bioEvents.length > 0) {
            for (BioEvent be : bioEvents) {
                ArrayList<EventAssignment> eventAssignments = be.getEventAssignments();
                if (eventAssignments != null) {
                    for (EventAssignment ea : eventAssignments) {
                        if (!eventAssignTargets.contains(ea.getTarget())) {
                            eventAssignTargets.add(ea.getTarget());
                        }
                    }
                }
            }
        }
        /**
         * @author anu : RATE RULES
         */
        RateRule[] rateRules = simContext.getRateRules();
        if (rateRules != null && rateRules.length > 0) {
            for (RateRule rr : rateRules) {
                SymbolTableEntry rateRuleVar = rr.getRateRuleVar();
                if (!rateRuleVarTargets.contains(rateRuleVar)) {
                    rateRuleVarTargets.add(rateRuleVar);
                }
            }
        }
        for (int j = 0; j < modelParameters.length; j++) {
            Expression modelParamExpr = modelParameters[j].getExpression();
            GeometryClass geometryClass = getDefaultGeometryClass(modelParamExpr);
            VCUnitDefinition paramUnit = modelParameters[j].getUnitDefinition();
            modelParamExpr = getIdentifierSubstitutions(modelParamExpr, paramUnit, geometryClass);
            if (eventAssignTargets.contains(modelParameters[j]) || rateRuleVarTargets.contains(modelParameters[j])) {
                EventAssignmentOrRateRuleInitParameter eap = null;
                try {
                    eap = addEventAssignmentOrRateRuleInitParameter(modelParameters[j], modelParamExpr, PARAMETER_ROLE_EVENTASSIGN_OR_RATERULE_INITCONDN, paramUnit);
                } catch (PropertyVetoException e) {
                    e.printStackTrace(System.out);
                    throw new MappingException(e.getMessage());
                }
                VolVariable volVar = new VolVariable(modelParameters[j].getName(), null);
                varHash.addVariable(volVar);
                eventOrRateRuleVolVarHash.put(volVar, eap);
                /**
                 * RATE RULES
                 */
                if (rateRuleVarTargets.contains(modelParameters[j])) {
                    RateRuleRateParameter rateParam = null;
                    try {
                        Expression origExp = simContext.getRateRule(modelParameters[j]).getRateRuleExpression();
                        VCUnitDefinition rateUnit = modelUnitSystem.getInstance_TBD();
                        if (paramUnit != null && !paramUnit.equals(modelUnitSystem.getInstance_TBD())) {
                            rateUnit = paramUnit.divideBy(timeUnit);
                        }
                        Expression rateExpr = getIdentifierSubstitutions(origExp, rateUnit, geometryClass);
                        rateParam = addRateRuleRateParameter(modelParameters[j], rateExpr, PARAMETER_ROLE_RATERULE_RATE, rateUnit);
                    } catch (PropertyVetoException e) {
                        e.printStackTrace(System.out);
                        throw new MappingException(e.getMessage());
                    }
                    rateRuleRateParamHash.put(volVar, rateParam);
                }
            } else {
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], geometryClass), modelParamExpr, geometryClass));
            }
        }
    } else {
        for (int j = 0; j < modelParameters.length; j++) {
            Expression modelParamExpr = modelParameters[j].getExpression();
            GeometryClass geometryClass = getDefaultGeometryClass(modelParamExpr);
            modelParamExpr = getIdentifierSubstitutions(modelParamExpr, modelParameters[j].getUnitDefinition(), geometryClass);
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], geometryClass), modelParamExpr, geometryClass));
        }
    }
    // 
    for (SimulationContextParameter scParameter : simContext.getSimulationContextParameters()) {
        Expression scParameterExpression = scParameter.getExpression();
        GeometryClass gc = getDefaultGeometryClass(scParameterExpression);
        varHash.addVariable(newFunctionOrConstant(getMathSymbol(scParameter, gc), getIdentifierSubstitutions(scParameter.getExpression(), scParameter.getUnitDefinition(), gc), gc));
    }
    // 
    for (DataSymbol dataSymbol : simContext.getDataContext().getDataSymbols()) {
        if (dataSymbol instanceof FieldDataSymbol) {
            FieldDataSymbol fieldDataSymbol = (FieldDataSymbol) dataSymbol;
            GeometryClass geometryClass = null;
            FieldFunctionArguments ffs = new FieldFunctionArguments(fieldDataSymbol.getExternalDataIdentifier().getName(), fieldDataSymbol.getFieldDataVarName(), new Expression(fieldDataSymbol.getFieldDataVarTime()), VariableType.getVariableTypeFromVariableTypeName(fieldDataSymbol.getFieldDataVarType()));
            Expression exp = new Expression(ffs.infix());
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(dataSymbol, geometryClass), getIdentifierSubstitutions(exp, dataSymbol.getUnitDefinition(), geometryClass), geometryClass));
        } else {
            throw new RuntimeException("In Application '" + simContext.getName() + "', dataSymbol type '" + dataSymbol.getClass().getName() + "' not yet supported for math generation");
        }
    }
    // 
    // 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("In Application '" + simContext.getName() + "', " + 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");
    }
    // 
    // volume variables
    // 
    Enumeration<SpeciesContextMapping> enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = enum1.nextElement();
        if (scm.getVariable() instanceof VolVariable) {
            if (!(mathDesc.getVariable(scm.getVariable().getName()) instanceof VolVariable)) {
                varHash.addVariable(scm.getVariable());
            }
        }
    }
    // 
    // membrane variables
    // 
    enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
        if (scm.getVariable() instanceof MemVariable) {
            varHash.addVariable(scm.getVariable());
        }
    }
    // 
    // volume region variables
    // 
    enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
        if (scm.getVariable() instanceof VolumeRegionVariable) {
            varHash.addVariable(scm.getVariable());
        }
    }
    // 
    // membrane region variables
    // 
    enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
        if (scm.getVariable() instanceof MembraneRegionVariable) {
            varHash.addVariable(scm.getVariable());
        }
    }
    // 
    // add compartment and membrane subdomains
    // 
    ArrayList<CompartmentSubdomainContext> compartmentSubdomainContexts = new ArrayList<CompartmentSubdomainContext>();
    ArrayList<MembraneSubdomainContext> membraneSubdomainContexts = new ArrayList<MembraneSubdomainContext>();
    addSubdomains(model, compartmentSubdomainContexts, membraneSubdomainContexts);
    // membrane velocities set on MembraneSubdomains later.
    addSpatialProcesses(varHash, compartmentSubdomainContexts, membraneSubdomainContexts);
    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)));
    // 
    // only calculate potential if at least one MembraneMapping has CalculateVoltage == true
    // 
    boolean bCalculatePotential = false;
    for (int i = 0; i < structureMappings.length; i++) {
        if (structureMappings[i] instanceof MembraneMapping) {
            if (((MembraneMapping) structureMappings[i]).getCalculateVoltage()) {
                bCalculatePotential = true;
            }
        }
    }
    potentialMapping = new PotentialMapping(simContext, this);
    if (bCalculatePotential) {
        potentialMapping.computeMath();
        // 
        // copy functions for currents and constants for capacitances
        // 
        ElectricalDevice[] devices = potentialMapping.getElectricalDevices();
        for (int j = 0; j < devices.length; j++) {
            if (devices[j] instanceof MembraneElectricalDevice) {
                MembraneElectricalDevice membraneElectricalDevice = (MembraneElectricalDevice) devices[j];
                MembraneMapping memMapping = membraneElectricalDevice.getMembraneMapping();
                Parameter specificCapacitanceParm = memMapping.getParameterFromRole(MembraneMapping.ROLE_SpecificCapacitance);
                varHash.addVariable(new Constant(getMathSymbol(specificCapacitanceParm, memMapping.getGeometryClass()), getIdentifierSubstitutions(specificCapacitanceParm.getExpression(), specificCapacitanceParm.getUnitDefinition(), memMapping.getGeometryClass())));
                ElectricalDevice.ElectricalDeviceParameter transmembraneCurrentParm = membraneElectricalDevice.getParameterFromRole(ElectricalDevice.ROLE_TransmembraneCurrent);
                ElectricalDevice.ElectricalDeviceParameter totalCurrentParm = membraneElectricalDevice.getParameterFromRole(ElectricalDevice.ROLE_TotalCurrent);
                ElectricalDevice.ElectricalDeviceParameter capacitanceParm = membraneElectricalDevice.getParameterFromRole(ElectricalDevice.ROLE_Capacitance);
                GeometryClass geometryClass = membraneElectricalDevice.getMembraneMapping().getGeometryClass();
                if (totalCurrentParm != null && /* totalCurrentDensityParm.getExpression()!=null && */
                memMapping.getCalculateVoltage()) {
                    Expression totalCurrentDensityExp = (totalCurrentParm.getExpression() != null) ? (totalCurrentParm.getExpression()) : (new Expression(0.0));
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrentParm, geometryClass), getIdentifierSubstitutions(totalCurrentDensityExp, totalCurrentParm.getUnitDefinition(), geometryClass), geometryClass));
                }
                if (transmembraneCurrentParm != null && transmembraneCurrentParm.getExpression() != null && memMapping.getCalculateVoltage()) {
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(transmembraneCurrentParm, geometryClass), getIdentifierSubstitutions(transmembraneCurrentParm.getExpression(), transmembraneCurrentParm.getUnitDefinition(), geometryClass), geometryClass));
                }
                if (capacitanceParm != null && capacitanceParm.getExpression() != null && memMapping.getCalculateVoltage()) {
                    StructureMappingParameter sizeParameter = membraneElectricalDevice.getMembraneMapping().getSizeParameter();
                    if (simContext.getGeometry().getDimension() == 0 && (sizeParameter.getExpression() == null || sizeParameter.getExpression().isZero())) {
                        varHash.addVariable(newFunctionOrConstant(getMathSymbol(capacitanceParm, geometryClass), getIdentifierSubstitutions(Expression.mult(memMapping.getNullSizeParameterValue(), specificCapacitanceParm.getExpression()), capacitanceParm.getUnitDefinition(), geometryClass), geometryClass));
                    } else {
                        varHash.addVariable(newFunctionOrConstant(getMathSymbol(capacitanceParm, geometryClass), getIdentifierSubstitutions(capacitanceParm.getExpression(), capacitanceParm.getUnitDefinition(), geometryClass), geometryClass));
                    }
                }
                // 
                if (membraneElectricalDevice.getDependentVoltageExpression() == null) {
                    // is Voltage Independent?
                    StructureMapping.StructureMappingParameter initialVoltageParm = memMapping.getInitialVoltageParameter();
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(initialVoltageParm, memMapping.getGeometryClass()), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), memMapping.getGeometryClass()), memMapping.getGeometryClass()));
                } else // 
                // membrane forced potential
                // 
                {
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(memMapping.getMembrane().getMembraneVoltage(), memMapping.getGeometryClass()), getIdentifierSubstitutions(membraneElectricalDevice.getDependentVoltageExpression(), memMapping.getMembrane().getMembraneVoltage().getUnitDefinition(), memMapping.getGeometryClass()), memMapping.getGeometryClass()));
                }
            } else if (devices[j] instanceof CurrentClampElectricalDevice) {
                CurrentClampElectricalDevice currentClampDevice = (CurrentClampElectricalDevice) devices[j];
                // total current = current source (no capacitance)
                Parameter totalCurrentParm = currentClampDevice.getParameterFromRole(CurrentClampElectricalDevice.ROLE_TotalCurrent);
                Parameter currentParm = currentClampDevice.getParameterFromRole(CurrentClampElectricalDevice.ROLE_TransmembraneCurrent);
                // Parameter dependentVoltage = currentClampDevice.getCurrentClampStimulus().getVoltageParameter();
                Feature deviceElectrodeFeature = currentClampDevice.getCurrentClampStimulus().getElectrode().getFeature();
                Feature groundElectrodeFeature = simContext.getGroundElectrode().getFeature();
                Membrane membrane = model.getStructureTopology().getMembrane(deviceElectrodeFeature, groundElectrodeFeature);
                GeometryClass geometryClass = null;
                if (membrane != null) {
                    StructureMapping membraneStructureMapping = simContext.getGeometryContext().getStructureMapping(membrane);
                    geometryClass = membraneStructureMapping.getGeometryClass();
                }
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrentParm, geometryClass), getIdentifierSubstitutions(totalCurrentParm.getExpression(), totalCurrentParm.getUnitDefinition(), geometryClass), geometryClass));
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(currentParm, geometryClass), getIdentifierSubstitutions(currentParm.getExpression(), currentParm.getUnitDefinition(), geometryClass), geometryClass));
                // varHash.addVariable(newFunctionOrConstant(getMathSymbol(dependentVoltage,null),getIdentifierSubstitutions(currentClampDevice.getDependentVoltageExpression(),dependentVoltage.getUnitDefinition(),null)));
                // 
                // add user-defined parameters
                // 
                ElectricalDevice.ElectricalDeviceParameter[] parameters = currentClampDevice.getParameters();
                for (int k = 0; k < parameters.length; k++) {
                    if (parameters[k].getExpression() != null) {
                        // guards against voltage parameters that are "variable".
                        varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[k], null), getIdentifierSubstitutions(parameters[k].getExpression(), parameters[k].getUnitDefinition(), geometryClass), geometryClass));
                    }
                }
            } else if (devices[j] instanceof VoltageClampElectricalDevice) {
                VoltageClampElectricalDevice voltageClampDevice = (VoltageClampElectricalDevice) devices[j];
                Feature deviceElectrodeFeature = voltageClampDevice.getVoltageClampStimulus().getElectrode().getFeature();
                Feature groundElectrodeFeature = simContext.getGroundElectrode().getFeature();
                Membrane membrane = model.getStructureTopology().getMembrane(deviceElectrodeFeature, groundElectrodeFeature);
                GeometryClass geometryClass = null;
                if (membrane != null) {
                    StructureMapping membraneStructureMapping = simContext.getGeometryContext().getStructureMapping(membrane);
                    geometryClass = membraneStructureMapping.getGeometryClass();
                }
                // total current = current source (no capacitance)
                Parameter totalCurrent = voltageClampDevice.getParameterFromRole(VoltageClampElectricalDevice.ROLE_TotalCurrent);
                Parameter totalCurrentParm = voltageClampDevice.getParameterFromRole(VoltageClampElectricalDevice.ROLE_TotalCurrent);
                Parameter voltageParm = voltageClampDevice.getParameterFromRole(VoltageClampElectricalDevice.ROLE_Voltage);
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrent, geometryClass), getIdentifierSubstitutions(totalCurrent.getExpression(), totalCurrent.getUnitDefinition(), geometryClass), geometryClass));
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrentParm, geometryClass), getIdentifierSubstitutions(totalCurrentParm.getExpression(), totalCurrentParm.getUnitDefinition(), geometryClass), geometryClass));
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(voltageParm, geometryClass), getIdentifierSubstitutions(voltageParm.getExpression(), voltageParm.getUnitDefinition(), geometryClass), geometryClass));
                // 
                // add user-defined parameters
                // 
                ElectricalDevice.ElectricalDeviceParameter[] parameters = voltageClampDevice.getParameters();
                for (int k = 0; k < parameters.length; k++) {
                    if (parameters[k].getRole() == ElectricalDevice.ROLE_UserDefined) {
                        varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[k], geometryClass), getIdentifierSubstitutions(parameters[k].getExpression(), parameters[k].getUnitDefinition(), geometryClass), geometryClass));
                    }
                }
            }
        }
    } else {
        // 
        for (int j = 0; j < structureMappings.length; j++) {
            if (structureMappings[j] instanceof MembraneMapping) {
                MembraneMapping memMapping = (MembraneMapping) structureMappings[j];
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(memMapping.getMembrane().getMembraneVoltage(), memMapping.getGeometryClass()), getIdentifierSubstitutions(memMapping.getInitialVoltageParameter().getExpression(), memMapping.getInitialVoltageParameter().getUnitDefinition(), memMapping.getGeometryClass()), memMapping.getGeometryClass()));
            }
        }
    }
    // 
    for (int j = 0; j < structureMappings.length; j++) {
        if (structureMappings[j] instanceof MembraneMapping) {
            MembraneMapping membraneMapping = (MembraneMapping) structureMappings[j];
            Membrane.MembraneVoltage membraneVoltage = membraneMapping.getMembrane().getMembraneVoltage();
            ElectricalDevice[] membraneDevices = potentialMapping.getElectricalDevices(membraneMapping.getMembrane());
            // ElectricalDevice membraneDevice = null;
            for (int i = 0; i < membraneDevices.length; i++) {
                if (membraneDevices[i].hasCapacitance() && membraneDevices[i].getDependentVoltageExpression() == null) {
                    GeometryClass geometryClass = membraneMapping.getGeometryClass();
                    if (geometryClass == null) {
                        throw new MappingException("Application '" + getSimulationContext().getName() + "'\nGeometry->StructureMapping->(" + structureMappings[j].getStructure().getTypeName() + ")'" + structureMappings[j].getStructure().getName() + "' must be mapped to geometry domain.\n(see 'Problems' tab)");
                    }
                    Domain domain = new Domain(geometryClass);
                    if (membraneMapping.getCalculateVoltage() && bCalculatePotential) {
                        if (geometryClass instanceof SurfaceClass) {
                            // 
                            if (mathDesc.getVariable(Membrane.MEMBRANE_VOLTAGE_REGION_NAME) == null) {
                                // varHash.addVariable(new MembraneRegionVariable(MembraneVoltage.MEMBRANE_VOLTAGE_REGION_NAME));
                                varHash.addVariable(new MembraneRegionVariable(getMathSymbol(membraneVoltage, geometryClass), domain));
                            }
                        } else {
                            // 
                            // spatially unresolved membrane, and must solve for potential ... make VolVariable for this compartment
                            // 
                            varHash.addVariable(new VolVariable(getMathSymbol(membraneVoltage, geometryClass), domain));
                        }
                        Parameter initialVoltageParm = membraneMapping.getInitialVoltageParameter();
                        Variable initVoltageFunction = newFunctionOrConstant(getMathSymbol(initialVoltageParm, geometryClass), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), geometryClass), geometryClass);
                        varHash.addVariable(initVoltageFunction);
                    } else {
                        // 
                        // don't calculate voltage, still may need it though
                        // 
                        Parameter initialVoltageParm = membraneMapping.getInitialVoltageParameter();
                        Variable voltageFunction = newFunctionOrConstant(getMathSymbol(membraneMapping.getMembrane().getMembraneVoltage(), geometryClass), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), geometryClass), geometryClass);
                        varHash.addVariable(voltageFunction);
                    }
                }
            }
        }
    }
    // 
    for (int j = 0; j < reactionSteps.length; j++) {
        ReactionStep rs = reactionSteps[j];
        if (simContext.getReactionContext().getReactionSpec(rs).isExcluded()) {
            continue;
        }
        Kinetics.KineticsParameter[] parameters = rs.getKinetics().getKineticsParameters();
        GeometryClass geometryClass = null;
        if (rs.getStructure() != null) {
            geometryClass = simContext.getGeometryContext().getStructureMapping(rs.getStructure()).getGeometryClass();
        }
        if (parameters != null) {
            for (int i = 0; i < parameters.length; i++) {
                if (((parameters[i].getRole() == Kinetics.ROLE_CurrentDensity) || (parameters[i].getRole() == Kinetics.ROLE_LumpedCurrent)) && (parameters[i].getExpression() == null || parameters[i].getExpression().isZero())) {
                    continue;
                }
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[i], geometryClass), getIdentifierSubstitutions(parameters[i].getExpression(), parameters[i].getUnitDefinition(), geometryClass), geometryClass));
            }
        }
    }
    // 
    // initial constants (either function or constant)
    // 
    SpeciesContextSpec[] speciesContextSpecs = simContext.getReactionContext().getSpeciesContextSpecs();
    for (int i = 0; i < speciesContextSpecs.length; i++) {
        // add initial count if present (!= null)
        SpeciesContextSpecParameter initCountParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialCount);
        SpeciesContext speciesContext = speciesContextSpecs[i].getSpeciesContext();
        if (initCountParm != null && initCountParm.getExpression() != null) {
            Expression initCountExpr = new Expression(initCountParm.getExpression());
            StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContext.getStructure());
            String[] symbols = initCountExpr.getSymbols();
            // Check if 'initExpr' has other speciesContexts in its expression, need to replace it with 'spContext_init'
            for (int j = 0; symbols != null && j < symbols.length; j++) {
                // if symbol is a speciesContext, replacing it with a reference to initial condition for that speciesContext.
                SpeciesContext spC = null;
                SymbolTableEntry ste = initCountExpr.getSymbolBinding(symbols[j]);
                if (ste instanceof SpeciesContextSpecProxyParameter) {
                    SpeciesContextSpecProxyParameter spspp = (SpeciesContextSpecProxyParameter) ste;
                    if (spspp.getTarget() instanceof SpeciesContext) {
                        spC = (SpeciesContext) spspp.getTarget();
                        SpeciesContextSpec spcspec = simContext.getReactionContext().getSpeciesContextSpec(spC);
                        SpeciesContextSpecParameter spCInitParm = spcspec.getParameterFromRole(SpeciesContextSpec.ROLE_InitialCount);
                        // need to get init condn expression, but can't get it from getMathSymbol() (mapping between bio and math), hence get it as below.
                        Expression scsInitExpr = new Expression(spCInitParm, getNameScope());
                        // scsInitExpr.bindExpression(this);
                        initCountExpr.substituteInPlace(new Expression(spC.getName()), scsInitExpr);
                    }
                }
            }
            // now create the appropriate function for the current speciesContextSpec.
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(initCountParm, sm.getGeometryClass()), getIdentifierSubstitutions(initCountExpr, initCountParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
        }
        // add initial concentration (may be derived from initial count if necessary)
        SpeciesContextSpecParameter initConcParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
        if (initConcParm != null) {
            Expression initConcExpr = null;
            if (initConcParm.getExpression() != null) {
                initConcExpr = new Expression(initConcParm.getExpression());
            } else if (initCountParm != null && initCountParm.getExpression() != null) {
                Expression structureSizeExpr = new Expression(speciesContext.getStructure().getStructureSize(), getNameScope());
                VCUnitDefinition concUnit = initConcParm.getUnitDefinition();
                VCUnitDefinition countDensityUnit = initCountParm.getUnitDefinition().divideBy(speciesContext.getStructure().getStructureSize().getUnitDefinition());
                Expression unitFactor = getUnitFactor(concUnit.divideBy(countDensityUnit));
                initConcExpr = Expression.mult(Expression.div(new Expression(initCountParm, getNameScope()), structureSizeExpr), unitFactor);
            }
            StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContext.getStructure());
            String[] symbols = initConcExpr.getSymbols();
            // Check if 'initExpr' has other speciesContexts in its expression, need to replace it with 'spContext_init'
            for (int j = 0; symbols != null && j < symbols.length; j++) {
                // if symbol is a speciesContext, replacing it with a reference to initial condition for that speciesContext.
                SpeciesContext spC = null;
                SymbolTableEntry ste = initConcExpr.getSymbolBinding(symbols[j]);
                if (ste instanceof SpeciesContextSpecProxyParameter) {
                    SpeciesContextSpecProxyParameter spspp = (SpeciesContextSpecProxyParameter) ste;
                    if (spspp.getTarget() instanceof SpeciesContext) {
                        spC = (SpeciesContext) spspp.getTarget();
                        SpeciesContextSpec spcspec = simContext.getReactionContext().getSpeciesContextSpec(spC);
                        SpeciesContextSpecParameter spCInitParm = spcspec.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
                        // if initConc param expression is null, try initCount
                        if (spCInitParm.getExpression() == null) {
                            spCInitParm = spcspec.getParameterFromRole(SpeciesContextSpec.ROLE_InitialCount);
                        }
                        // need to get init condn expression, but can't get it from getMathSymbol() (mapping between bio and math), hence get it as below.
                        Expression scsInitExpr = new Expression(spCInitParm, getNameScope());
                        // scsInitExpr.bindExpression(this);
                        initConcExpr.substituteInPlace(new Expression(spC.getName()), scsInitExpr);
                    }
                }
            }
            // now create the appropriate function for the current speciesContextSpec.
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(initConcParm, sm.getGeometryClass()), getIdentifierSubstitutions(initConcExpr, initConcParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
        }
    }
    // 
    for (int i = 0; i < speciesContextSpecs.length; i++) {
        SpeciesContextMapping scm = getSpeciesContextMapping(speciesContextSpecs[i].getSpeciesContext());
        SpeciesContextSpec.SpeciesContextSpecParameter diffParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_DiffusionRate);
        if (diffParm != null && (scm.isPDERequired())) {
            StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(diffParm, sm.getGeometryClass()), getIdentifierSubstitutions(diffParm.getExpression(), diffParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
        }
    }
    // 
    for (int i = 0; i < speciesContextSpecs.length; i++) {
        SpeciesContextSpec.SpeciesContextSpecParameter bc_xm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueXm);
        StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
        if (bc_xm != null && (bc_xm.getExpression() != null)) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_xm, sm.getGeometryClass()), getIdentifierSubstitutions(bc_xm.getExpression(), bc_xm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
        }
        SpeciesContextSpec.SpeciesContextSpecParameter bc_xp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueXp);
        if (bc_xp != null && (bc_xp.getExpression() != null)) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_xp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_xp.getExpression(), bc_xp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
        }
        SpeciesContextSpec.SpeciesContextSpecParameter bc_ym = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueYm);
        if (bc_ym != null && (bc_ym.getExpression() != null)) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_ym, sm.getGeometryClass()), getIdentifierSubstitutions(bc_ym.getExpression(), bc_ym.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
        }
        SpeciesContextSpec.SpeciesContextSpecParameter bc_yp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueYp);
        if (bc_yp != null && (bc_yp.getExpression() != null)) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_yp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_yp.getExpression(), bc_yp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
        }
        SpeciesContextSpec.SpeciesContextSpecParameter bc_zm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueZm);
        if (bc_zm != null && (bc_zm.getExpression() != null)) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_zm, sm.getGeometryClass()), getIdentifierSubstitutions(bc_zm.getExpression(), bc_zm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
        }
        SpeciesContextSpec.SpeciesContextSpecParameter bc_zp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueZp);
        if (bc_zp != null && (bc_zp.getExpression() != null)) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_zp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_zp.getExpression(), bc_zp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
        }
    }
    // 
    for (int i = 0; i < speciesContextSpecs.length; i++) {
        SpeciesContextSpec.SpeciesContextSpecParameter advection_velX = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityX);
        StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
        GeometryClass geometryClass = sm.getGeometryClass();
        if (advection_velX != null && (advection_velX.getExpression() != null)) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velX, geometryClass), getIdentifierSubstitutions(advection_velX.getExpression(), advection_velX.getUnitDefinition(), geometryClass), geometryClass));
        }
        SpeciesContextSpec.SpeciesContextSpecParameter advection_velY = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityY);
        if (advection_velY != null && (advection_velY.getExpression() != null)) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velY, geometryClass), getIdentifierSubstitutions(advection_velY.getExpression(), advection_velY.getUnitDefinition(), geometryClass), geometryClass));
        }
        SpeciesContextSpec.SpeciesContextSpecParameter advection_velZ = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityZ);
        if (advection_velZ != null && (advection_velZ.getExpression() != null)) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velZ, geometryClass), getIdentifierSubstitutions(advection_velZ.getExpression(), advection_velZ.getUnitDefinition(), geometryClass), geometryClass));
        }
    }
    // 
    // 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());
        }
    }
    // 
    // conversion factors
    // 
    varHash.addVariable(new Constant(model.getKMOLE().getName(), getIdentifierSubstitutions(model.getKMOLE().getExpression(), model.getKMOLE().getUnitDefinition(), null)));
    varHash.addVariable(new Constant(model.getN_PMOLE().getName(), getIdentifierSubstitutions(model.getN_PMOLE().getExpression(), model.getN_PMOLE().getUnitDefinition(), null)));
    varHash.addVariable(new Constant(model.getKMILLIVOLTS().getName(), getIdentifierSubstitutions(model.getKMILLIVOLTS().getExpression(), model.getKMILLIVOLTS().getUnitDefinition(), null)));
    varHash.addVariable(new Constant(model.getK_GHK().getName(), getIdentifierSubstitutions(model.getK_GHK().getExpression(), model.getK_GHK().getUnitDefinition(), null)));
    // 
    for (int i = 0; i < structureMappings.length; i++) {
        StructureMapping sm = structureMappings[i];
        if (simContext.getGeometry().getDimension() == 0) {
            StructureMappingParameter sizeParm = sm.getSizeParameter();
            if (sizeParm != null && sizeParm.getExpression() != null) {
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(sizeParm, sm.getGeometryClass()), getIdentifierSubstitutions(sizeParm.getExpression(), sizeParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
            } else {
                if (sm instanceof MembraneMapping) {
                    MembraneMapping mm = (MembraneMapping) sm;
                    StructureMappingParameter volFrac = mm.getVolumeFractionParameter();
                    if (volFrac != null && volFrac.getExpression() != null) {
                        varHash.addVariable(newFunctionOrConstant(getMathSymbol(volFrac, sm.getGeometryClass()), getIdentifierSubstitutions(volFrac.getExpression(), volFrac.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
                    }
                    StructureMappingParameter surfToVol = mm.getSurfaceToVolumeParameter();
                    if (surfToVol != null && surfToVol.getExpression() != null) {
                        varHash.addVariable(newFunctionOrConstant(getMathSymbol(surfToVol, sm.getGeometryClass()), getIdentifierSubstitutions(surfToVol.getExpression(), surfToVol.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
                    }
                }
            }
        } else {
            Parameter parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitArea);
            if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
            }
            parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitVolume);
            if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
            }
            parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitArea);
            if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
            }
            parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitVolume);
            if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
            }
        }
        StructureMappingParameter sizeParm = sm.getSizeParameter();
        if (sm.getGeometryClass() != null && sizeParm != null) {
            if (simContext.getGeometry().getDimension() == 0) {
                if (sizeParm.getExpression() != null) {
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(sizeParm, sm.getGeometryClass()), getIdentifierSubstitutions(sizeParm.getExpression(), sizeParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
                }
            } else {
                String compartmentName = sm.getGeometryClass().getName();
                VCUnitDefinition sizeUnit = sm.getSizeParameter().getUnitDefinition();
                String sizeFunctionName = null;
                if (sm instanceof MembraneMapping) {
                    MembraneMapping mm = (MembraneMapping) sm;
                    if (mm.getGeometryClass() instanceof SurfaceClass) {
                        sizeFunctionName = MathFunctionDefinitions.Function_regionArea_current.getFunctionName();
                    } else if (mm.getGeometryClass() instanceof SubVolume) {
                        sizeFunctionName = MathFunctionDefinitions.Function_regionVolume_current.getFunctionName();
                    }
                } else if (sm instanceof FeatureMapping) {
                    sizeFunctionName = MathFunctionDefinitions.Function_regionVolume_current.getFunctionName();
                } else {
                    throw new RuntimeException("structure mapping " + sm.getClass().getName() + " not yet supported");
                }
                Expression totalVolumeCorrection = sm.getStructureSizeCorrection(simContext, this);
                Expression sizeFunctionExpression = Expression.function(sizeFunctionName, new Expression[] { new Expression("'" + compartmentName + "'") });
                // sizeFunctionExpression.bindExpression(mathDesc);
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(sizeParm, sm.getGeometryClass()), getIdentifierSubstitutions(Expression.mult(totalVolumeCorrection, sizeFunctionExpression), sizeUnit, sm.getGeometryClass()), sm.getGeometryClass()));
            }
        }
    }
    // 
    for (SymbolTableEntry rrvSTE : rateRuleVarTargets) {
        if (rrvSTE instanceof SpeciesContext) {
            SpeciesContext rateRuleSpContext = (SpeciesContext) rrvSTE;
            // check if speciesContext has already been added as a vol/membrane var
            SpeciesContextMapping scm = getSpeciesContextMapping(rateRuleSpContext);
            if (scm.getVariable() instanceof VolVariable || scm.getVariable() instanceof MemVariable) {
                throw new RuntimeException("SpeciesContext '" + rrvSTE.getName() + "' has an equation and is also a rate rule variable, which is not allowed.");
            }
            // get the initial condition expression of the speciesContext.
            SpeciesContextSpec speciesContextSpec = simContext.getReactionContext().getSpeciesContextSpec(rateRuleSpContext);
            SpeciesContextSpecParameter scsInitParam = speciesContextSpec.getInitialConditionParameter();
            Expression scInitExpr = scsInitParam.getExpression();
            GeometryClass geometryClass = getDefaultGeometryClass(scInitExpr);
            VCUnitDefinition scsInitParamUnit = scsInitParam.getUnitDefinition();
            scInitExpr = getIdentifierSubstitutions(scInitExpr, scsInitParamUnit, geometryClass);
            EventAssignmentOrRateRuleInitParameter eap = null;
            try {
                // create an eventAssgnmentOrRateRuleInitParameter for the rate rule variable
                eap = addEventAssignmentOrRateRuleInitParameter(rateRuleSpContext, scInitExpr, PARAMETER_ROLE_EVENTASSIGN_OR_RATERULE_INITCONDN, scsInitParamUnit);
            } catch (PropertyVetoException e) {
                e.printStackTrace(System.out);
                throw new MappingException(e.getMessage());
            }
            // create a volVariable for this speciesContext (shouldn't have one already - since a speciesContext that has a rate rule should not participate in a reaction.
            VolVariable volVar = new VolVariable(rateRuleSpContext.getName(), null);
            varHash.addVariable(volVar);
            eventOrRateRuleVolVarHash.put(volVar, eap);
            // create the rate parameter
            RateRuleRateParameter rateParam = null;
            try {
                Expression origExp = simContext.getRateRule(rateRuleSpContext).getRateRuleExpression();
                VCUnitDefinition rateUnit = modelUnitSystem.getInstance_TBD();
                if (scsInitParamUnit != null && !scsInitParamUnit.equals(modelUnitSystem.getInstance_TBD())) {
                    rateUnit = scsInitParamUnit.divideBy(timeUnit);
                }
                Expression rateExpr = getIdentifierSubstitutions(origExp, rateUnit, geometryClass);
                rateParam = addRateRuleRateParameter(rateRuleSpContext, rateExpr, PARAMETER_ROLE_RATERULE_RATE, rateUnit);
            } catch (PropertyVetoException e) {
                e.printStackTrace(System.out);
                throw new MappingException(e.getMessage());
            }
            rateRuleRateParamHash.put(volVar, rateParam);
        }
    // end if (ste instanceof SC)
    }
    // 
    for (int i = 0; i < fieldMathMappingParameters.length; i++) {
        GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
        varHash.addVariable(newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass()));
    }
    // 
    // functions
    // 
    enum1 = getSpeciesContextMappings();
    RateRule[] rateRules = simContext.getRateRules();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
        if (scm.getVariable() == null && scm.getDependencyExpression() != null) {
            // check if speciesContext has a rateRule; then the speciesContext should not be added as a constant
            if (rateRules == null) {
                if (simContext.getRateRule(scm.getSpeciesContext()) == null) {
                    StructureMapping sm = simContext.getGeometryContext().getStructureMapping(scm.getSpeciesContext().getStructure());
                    if (sm.getGeometryClass() == null) {
                        Structure s = sm.getStructure();
                        if (s != null) {
                            throw new RuntimeException("unmapped structure " + s.getName());
                        }
                        throw new RuntimeException("structure mapping with no structure or mapping");
                    }
                    Variable dependentVariable = newFunctionOrConstant(getMathSymbol(scm.getSpeciesContext(), sm.getGeometryClass()), getIdentifierSubstitutions(scm.getDependencyExpression(), scm.getSpeciesContext().getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass());
                    dependentVariable.setDomain(new Domain(sm.getGeometryClass()));
                    varHash.addVariable(dependentVariable);
                }
            }
        }
    }
    BioEvent[] bioevents = simContext.getBioEvents();
    if (bioevents != null && bioevents.length > 0) {
        for (BioEvent be : bioevents) {
            // transform the bioEvent trigger/delay to math Event
            for (LocalParameter p : be.getEventParameters()) {
                if (p.getExpression() != null) {
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(p, null), getIdentifierSubstitutions(p.getExpression(), p.getUnitDefinition(), null), null));
                } else if (be.getParameter(BioEventParameterType.GeneralTriggerFunction) == p) {
                    // 
                    // use generated function here.
                    // 
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(p, null), getIdentifierSubstitutions(be.generateTriggerExpression(), p.getUnitDefinition(), null), null));
                }
            }
        }
    }
    // 
    // set Variables to MathDescription all at once with the order resolved by "VariableHash"
    // 
    mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
    // 
    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");
    }
    // 
    for (CompartmentSubdomainContext compartmentSubDomainContext : compartmentSubdomainContexts) {
        SubVolume subVolume = compartmentSubDomainContext.subvolume;
        CompartmentSubDomain subDomain = mathDesc.getCompartmentSubDomain(subVolume.getName());
        // 
        // assign boundary condition types
        // 
        StructureMapping[] mappedSMs = simContext.getGeometryContext().getStructureMappings(subVolume);
        FeatureMapping mappedFM = null;
        for (int i = 0; i < mappedSMs.length; i++) {
            if (mappedSMs[i] instanceof FeatureMapping) {
                if (mappedFM != null) {
                    lg.warn("WARNING:::: MathMapping.refreshMathDescription() ... assigning boundary condition types not unique");
                }
                mappedFM = (FeatureMapping) mappedSMs[i];
            }
        }
        if (mappedFM != null) {
            if (simContext.getGeometry().getDimension() > 0) {
                subDomain.setBoundaryConditionXm(mappedFM.getBoundaryConditionTypeXm());
                subDomain.setBoundaryConditionXp(mappedFM.getBoundaryConditionTypeXp());
            }
            if (simContext.getGeometry().getDimension() > 1) {
                subDomain.setBoundaryConditionYm(mappedFM.getBoundaryConditionTypeYm());
                subDomain.setBoundaryConditionYp(mappedFM.getBoundaryConditionTypeYp());
            }
            if (simContext.getGeometry().getDimension() > 2) {
                subDomain.setBoundaryConditionZm(mappedFM.getBoundaryConditionTypeZm());
                subDomain.setBoundaryConditionZp(mappedFM.getBoundaryConditionTypeZp());
            }
        }
        // 
        // create equations
        // 
        VolumeStructureAnalyzer structureAnalyzer = getVolumeStructureAnalyzer(subVolume);
        Enumeration<SpeciesContextMapping> enumSCM = getSpeciesContextMappings();
        while (enumSCM.hasMoreElements()) {
            SpeciesContextMapping scm = enumSCM.nextElement();
            SpeciesContext sc = scm.getSpeciesContext();
            StructureMapping sm = simContext.getGeometryContext().getStructureMapping(sc.getStructure());
            SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(sc);
            // 
            // if an independent volume variable, then create equation for it (if mapped to this subDomain)
            // 
            final GeometryClass gc = sm.getGeometryClass();
            if (gc == null || !gc.getName().equals(subDomain.getName())) {
                continue;
            }
            SpeciesContextSpecParameter initConcParameter = scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
            if ((scm.getVariable() instanceof VolumeRegionVariable) && scm.getDependencyExpression() == null) {
                VolumeRegionVariable volumeRegionVariable = (VolumeRegionVariable) scm.getVariable();
                Expression initial = getIdentifierSubstitutions(new Expression(initConcParameter, getNameScope()), initConcParameter.getUnitDefinition(), sm.getGeometryClass());
                Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass());
                VolumeRegionEquation volumeRegionEquation = new VolumeRegionEquation(volumeRegionVariable, initial);
                volumeRegionEquation.setVolumeRateExpression(rate);
                subDomain.addEquation(volumeRegionEquation);
            } else if (scm.getVariable() instanceof VolVariable && scm.getDependencyExpression() == null) {
                VolVariable variable = (VolVariable) scm.getVariable();
                Equation equation = null;
                if (sm.getGeometryClass() == subVolume) {
                    if (scm.isPDERequired()) {
                        // 
                        // species context belongs to this subDomain
                        // 
                        Expression initial = getIdentifierSubstitutions(new Expression(initConcParameter, getNameScope()), initConcParameter.getUnitDefinition(), sm.getGeometryClass());
                        Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass());
                        SpeciesContextSpecParameter diffusionParameter = scs.getDiffusionParameter();
                        Expression diffusion = getIdentifierSubstitutions(new Expression(diffusionParameter, getNameScope()), diffusionParameter.getUnitDefinition(), sm.getGeometryClass());
                        equation = new PdeEquation(variable, initial, rate, diffusion);
                        ((PdeEquation) equation).setBoundaryXm((scs.getBoundaryXmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXmParameter(), sm.getGeometryClass())));
                        ((PdeEquation) equation).setBoundaryXp((scs.getBoundaryXpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXpParameter(), sm.getGeometryClass())));
                        ((PdeEquation) equation).setBoundaryYm((scs.getBoundaryYmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYmParameter(), sm.getGeometryClass())));
                        ((PdeEquation) equation).setBoundaryYp((scs.getBoundaryYpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYpParameter(), sm.getGeometryClass())));
                        ((PdeEquation) equation).setBoundaryZm((scs.getBoundaryZmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZmParameter(), sm.getGeometryClass())));
                        ((PdeEquation) equation).setBoundaryZp((scs.getBoundaryZpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZpParameter(), sm.getGeometryClass())));
                        if (simContext.getGeometry().getDimension() >= 1) {
                            Expression velXExp = null;
                            if (scs.getVelocityXParameter().getExpression() != null) {
                                velXExp = new Expression(getMathSymbol(scs.getVelocityXParameter(), sm.getGeometryClass()));
                            } else {
                                SpatialQuantity[] velX_quantities = scs.getVelocityQuantities(QuantityComponent.X);
                                if (velX_quantities.length > 0) {
                                    int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(subVolume).length;
                                    if (velX_quantities.length == 1 && numRegions == 1) {
                                        velXExp = new Expression(getMathSymbol(velX_quantities[0], sm.getGeometryClass()));
                                    } else {
                                        throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
                                    }
                                }
                            }
                            ((PdeEquation) equation).setVelocityX(velXExp);
                        }
                        if (simContext.getGeometry().getDimension() >= 2) {
                            Expression velYExp = null;
                            if (scs.getVelocityYParameter().getExpression() != null) {
                                velYExp = new Expression(getMathSymbol(scs.getVelocityYParameter(), sm.getGeometryClass()));
                            } else {
                                SpatialQuantity[] velY_quantities = scs.getVelocityQuantities(QuantityComponent.Y);
                                if (velY_quantities.length > 0) {
                                    int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(subVolume).length;
                                    if (velY_quantities.length == 1 && numRegions == 1) {
                                        velYExp = new Expression(getMathSymbol(velY_quantities[0], sm.getGeometryClass()));
                                    } else {
                                        throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
                                    }
                                }
                            }
                            ((PdeEquation) equation).setVelocityY(velYExp);
                        }
                        if (simContext.getGeometry().getDimension() == 3) {
                            Expression velZExp = null;
                            if (scs.getVelocityZParameter().getExpression() != null) {
                                velZExp = new Expression(getMathSymbol(scs.getVelocityZParameter(), sm.getGeometryClass()));
                            } else {
                                SpatialQuantity[] velZ_quantities = scs.getVelocityQuantities(QuantityComponent.Z);
                                if (velZ_quantities.length > 0) {
                                    int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(subVolume).length;
                                    if (velZ_quantities.length == 1 && numRegions == 1) {
                                        velZExp = new Expression(getMathSymbol(velZ_quantities[0], sm.getGeometryClass()));
                                    } else {
                                        throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
                                    }
                                }
                            }
                            ((PdeEquation) equation).setVelocityZ(velZExp);
                        }
                        subDomain.replaceEquation(equation);
                    } else {
                        // 
                        // ODE - species context belongs to this subDomain
                        // 
                        Expression initial = new Expression(getMathSymbol(initConcParameter, null));
                        Expression rate = (scm.getRate() == null) ? new Expression(0.0) : getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass());
                        equation = new OdeEquation(variable, initial, rate);
                        subDomain.replaceEquation(equation);
                    }
                }
            }
        }
        // 
        // create fast system (if neccessary)
        // 
        SpeciesContextMapping[] fastSpeciesContextMappings = structureAnalyzer.getFastSpeciesContextMappings();
        if (fastSpeciesContextMappings != null) {
            FastSystem fastSystem = new FastSystem(mathDesc);
            for (int i = 0; i < fastSpeciesContextMappings.length; i++) {
                SpeciesContextMapping scm = fastSpeciesContextMappings[i];
                if (scm.getFastInvariant() == null) {
                    // 
                    // independant-fast variable, create a fastRate object
                    // 
                    Expression rate = getIdentifierSubstitutions(scm.getFastRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), subVolume);
                    FastRate fastRate = new FastRate(rate);
                    fastSystem.addFastRate(fastRate);
                } else {
                    // 
                    // dependant-fast variable, create a fastInvariant object
                    // 
                    Expression rate = getIdentifierSubstitutions(scm.getFastInvariant(), modelUnitSystem.getVolumeConcentrationUnit(), subVolume);
                    FastInvariant fastInvariant = new FastInvariant(rate);
                    fastSystem.addFastInvariant(fastInvariant);
                }
            }
            subDomain.setFastSystem(fastSystem);
            // constructor calls the 'refresh' method which constructs depemdency matrix, dependent/independent vars and pseudoconstants, etc.
            // FastSystemAnalyzer fs_analyzer =
            new FastSystemAnalyzer(fastSystem, mathDesc);
        }
        // 
        // create ode's for voltages to be calculated on unresolved membranes mapped to this subVolume
        // 
        Structure[] localStructures = simContext.getGeometryContext().getStructuresFromGeometryClass(subVolume);
        for (int sIndex = 0; sIndex < localStructures.length; sIndex++) {
            if (localStructures[sIndex] instanceof Membrane) {
                Membrane membrane = (Membrane) localStructures[sIndex];
                MembraneMapping membraneMapping = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(membrane);
                if ((membraneMapping.getGeometryClass() instanceof SubVolume) && membraneMapping.getCalculateVoltage()) {
                    MembraneElectricalDevice capacitiveDevice = potentialMapping.getCapacitiveDevice(membrane);
                    if (capacitiveDevice.getDependentVoltageExpression() == null) {
                        VolVariable vVar = (VolVariable) mathDesc.getVariable(getMathSymbol(capacitiveDevice.getVoltageSymbol(), membraneMapping.getGeometryClass()));
                        Expression initExp = new Expression(getMathSymbol(capacitiveDevice.getMembraneMapping().getInitialVoltageParameter(), membraneMapping.getGeometryClass()));
                        subDomain.addEquation(new OdeEquation(vVar, initExp, getIdentifierSubstitutions(potentialMapping.getOdeRHS(capacitiveDevice, this), membrane.getMembraneVoltage().getUnitDefinition().divideBy(timeUnit), membraneMapping.getGeometryClass())));
                    } else {
                    // 
                    // 
                    // 
                    }
                }
            }
        }
    }
    // 
    for (MembraneSubdomainContext memSubdomainContext : membraneSubdomainContexts) {
        MembraneSubDomain memSubDomain = memSubdomainContext.membraneSubdomain;
        SurfaceClass surfaceClass = memSubdomainContext.surfaceClass;
        for (SurfaceRegionObject surfaceRegionObject : memSubdomainContext.surfaceRegionObjects) {
            if (surfaceRegionObject.isQuantityCategoryEnabled(QuantityCategory.SurfaceVelocity)) {
                int dim = simContext.getGeometry().getDimension();
                if (dim != 2) {
                    throw new MappingException("Membrane Velocity only supported for 2D geometries");
                }
                if (simContext.getGeometry().getDimension() >= 1) {
                    SpatialQuantity velXQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.SurfaceVelocity, QuantityComponent.X);
                    Expression velXExp = new Expression(velXQuantity, simContext.getNameScope());
                    memSubDomain.setVelocityX(getIdentifierSubstitutions(velXExp, velXQuantity.getUnitDefinition(), surfaceClass));
                }
                if (simContext.getGeometry().getDimension() >= 2) {
                    SpatialQuantity velYQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.SurfaceVelocity, QuantityComponent.Y);
                    Expression velYExp = new Expression(velYQuantity, simContext.getNameScope());
                    memSubDomain.setVelocityY(getIdentifierSubstitutions(velYExp, velYQuantity.getUnitDefinition(), surfaceClass));
                }
                if (simContext.getGeometry().getDimension() == 3) {
                    SpatialQuantity velZQuantity = surfaceRegionObject.getSpatialQuantity(QuantityCategory.SurfaceVelocity, QuantityComponent.Z);
                    Expression velZExp = new Expression(velZQuantity, simContext.getNameScope());
                    // memSubDomain.setVelocityZ(getIdentifierSubstitutions(velZExp, velZQuantity.getUnitDefinition(), surfaceClass));
                    throw new MappingException("Membrane Velocity not supported for 2D problems");
                }
            }
        }
        // 
        // create equations for membrane-bound molecular species
        // 
        MembraneStructureAnalyzer membraneStructureAnalyzer = getMembraneStructureAnalyzer(surfaceClass);
        Enumeration<SpeciesContextMapping> enumSCM = getSpeciesContextMappings();
        while (enumSCM.hasMoreElements()) {
            SpeciesContextMapping scm = enumSCM.nextElement();
            SpeciesContext sc = scm.getSpeciesContext();
            SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(sc);
            StructureMapping sm = simContext.getGeometryContext().getStructureMapping(sc.getStructure());
            // 
            if ((scm.getVariable() instanceof MembraneRegionVariable) && scm.getDependencyExpression() == null) {
                MembraneRegionEquation equation = null;
                MembraneRegionVariable memRegionVar = (MembraneRegionVariable) scm.getVariable();
                if (sm.getGeometryClass() == surfaceClass) {
                    // 
                    // species context belongs to this subDomain
                    // 
                    Expression initial = new Expression(getMathSymbol(scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration), sm.getGeometryClass()));
                    Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass());
                    equation = new MembraneRegionEquation(memRegionVar, initial);
                    equation.setMembraneRateExpression(rate);
                    // equation.setUniformRateExpression(newUniformRateExpression);
                    memSubDomain.replaceEquation(equation);
                }
            } else if ((scm.getVariable() instanceof MemVariable) && scm.getDependencyExpression() == null) {
                // 
                if (sm.getGeometryClass() == surfaceClass) {
                    Equation equation = null;
                    MemVariable variable = (MemVariable) scm.getVariable();
                    if (scm.isPDERequired()) {
                        // 
                        // PDE
                        // 
                        // 
                        // species context belongs to this subDomain
                        // 
                        Expression initial = new Expression(getMathSymbol(scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration), sm.getGeometryClass()));
                        Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass());
                        Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), sm.getGeometryClass()));
                        equation = new PdeEquation(variable, initial, rate, diffusion);
                        ((PdeEquation) equation).setBoundaryXm((scs.getBoundaryXmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXmParameter(), sm.getGeometryClass())));
                        ((PdeEquation) equation).setBoundaryXp((scs.getBoundaryXpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXpParameter(), sm.getGeometryClass())));
                        ((PdeEquation) equation).setBoundaryYm((scs.getBoundaryYmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYmParameter(), sm.getGeometryClass())));
                        ((PdeEquation) equation).setBoundaryYp((scs.getBoundaryYpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYpParameter(), sm.getGeometryClass())));
                        ((PdeEquation) equation).setBoundaryZm((scs.getBoundaryZmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZmParameter(), sm.getGeometryClass())));
                        ((PdeEquation) equation).setBoundaryZp((scs.getBoundaryZpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZpParameter(), sm.getGeometryClass())));
                        memSubDomain.replaceEquation(equation);
                    } else {
                        // 
                        // ODE
                        // 
                        // 
                        // species context belongs to this subDomain
                        // 
                        Expression initial = new Expression(getMathSymbol(scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration), null));
                        Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass());
                        equation = new OdeEquation(variable, initial, rate);
                        memSubDomain.replaceEquation(equation);
                    }
                }
            }
        }
        Enumeration<SpeciesContextMapping> enum_scm = getSpeciesContextMappings();
        while (enum_scm.hasMoreElements()) {
            SpeciesContextMapping scm = enum_scm.nextElement();
            if (scm.isPDERequired() || scm.getVariable() instanceof VolumeRegionVariable) {
                // Species species = scm.getSpeciesContext().getSpecies();
                Variable var = scm.getVariable();
                final Domain dm = var.getDomain();
                if (dm != null) {
                    final String domainName = dm.getName();
                    if (sameName(domainName, memSubDomain.getInsideCompartment()) || sameName(domainName, memSubDomain.getOutsideCompartment())) {
                        JumpCondition jc = memSubDomain.getJumpCondition(var);
                        if (jc == null) {
                            // System.out.println("MathMapping.refreshMathDescription(), adding jump condition for diffusing variable "+var.getName()+" on membrane "+membraneStructureAnalyzer.getMembrane().getName());
                            if (var instanceof VolVariable) {
                                jc = new JumpCondition((VolVariable) var);
                            } else if (var instanceof VolumeRegionVariable) {
                                jc = new JumpCondition((VolumeRegionVariable) var);
                            } else {
                                throw new RuntimeException("unexpected Variable type " + var.getClass().getName());
                            }
                            memSubDomain.addJumpCondition(jc);
                        }
                    }
                }
            }
        }
        // 
        // set jump conditions for any volume variables or volume region variables that have explicitly defined fluxes
        // 
        ResolvedFlux[] resolvedFluxes = membraneStructureAnalyzer.getResolvedFluxes();
        if (resolvedFluxes != null) {
            for (int i = 0; i < resolvedFluxes.length; i++) {
                SpeciesContext sc = resolvedFluxes[i].getSpeciesContext();
                SpeciesContextMapping scm = getSpeciesContextMapping(sc);
                StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure());
                if (scm.getVariable() instanceof VolVariable && scm.isPDERequired()) {
                    VolVariable volVar = (VolVariable) scm.getVariable();
                    JumpCondition jc = memSubDomain.getJumpCondition(volVar);
                    if (jc == null) {
                        jc = new JumpCondition(volVar);
                        memSubDomain.addJumpCondition(jc);
                    }
                    Expression flux = getIdentifierSubstitutions(resolvedFluxes[i].getFluxExpression(), resolvedFluxes[i].getUnitDefinition(), membraneStructureAnalyzer.getSurfaceClass());
                    if (memSubDomain.getInsideCompartment().getName().equals(sm.getGeometryClass().getName())) {
                        jc.setInFlux(flux);
                    } else if (memSubDomain.getOutsideCompartment().getName().equals(sm.getGeometryClass().getName())) {
                        jc.setOutFlux(flux);
                    } else {
                        throw new RuntimeException("Application  " + simContext.getName() + " : " + scm.getSpeciesContext().getName() + " has spatially resolved flux at membrane " + scm.getSpeciesContext().getStructure().getName() + " with a non-local flux species " + scm.getSpeciesContext().getName());
                    }
                } else if (scm.getVariable() instanceof VolumeRegionVariable) {
                    VolumeRegionVariable volRegionVar = (VolumeRegionVariable) scm.getVariable();
                    JumpCondition jc = memSubDomain.getJumpCondition(volRegionVar);
                    if (jc == null) {
                        jc = new JumpCondition(volRegionVar);
                        memSubDomain.addJumpCondition(jc);
                    }
                    Expression flux = getIdentifierSubstitutions(resolvedFluxes[i].getFluxExpression(), resolvedFluxes[i].getUnitDefinition(), membraneStructureAnalyzer.getSurfaceClass());
                    if (memSubDomain.getInsideCompartment().getName().equals(sm.getGeometryClass().getName())) {
                        jc.setInFlux(flux);
                    } else if (memSubDomain.getOutsideCompartment().getName().equals(sm.getGeometryClass().getName())) {
                        jc.setOutFlux(flux);
                    } else {
                        throw new RuntimeException("Application  " + simContext.getName() + " : " + scm.getSpeciesContext().getName() + " has spatially resolved flux at membrane " + scm.getSpeciesContext().getStructure().getName() + " with a non-local flux species " + scm.getSpeciesContext().getName());
                    }
                } else {
                    throw new MappingException("Application  " + simContext.getName() + " : " + scm.getSpeciesContext().getName() + " has spatially resolved flux at membrane " + scm.getSpeciesContext().getStructure().getName() + ", but doesn't diffuse in compartment " + scm.getSpeciesContext().getStructure().getName());
                }
            }
        }
        // 
        // create fast system (if neccessary)
        // 
        SpeciesContextMapping[] fastSpeciesContextMappings = membraneStructureAnalyzer.getFastSpeciesContextMappings();
        if (fastSpeciesContextMappings != null) {
            FastSystem fastSystem = new FastSystem(mathDesc);
            for (int i = 0; i < fastSpeciesContextMappings.length; i++) {
                SpeciesContextMapping scm = fastSpeciesContextMappings[i];
                if (scm.getFastInvariant() == null) {
                    // 
                    // independant-fast variable, create a fastRate object
                    // 
                    VCUnitDefinition rateUnit = scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit);
                    FastRate fastRate = new FastRate(getIdentifierSubstitutions(scm.getFastRate(), rateUnit, surfaceClass));
                    fastSystem.addFastRate(fastRate);
                } else {
                    // 
                    // dependant-fast variable, create a fastInvariant object
                    // 
                    VCUnitDefinition invariantUnit = scm.getSpeciesContext().getUnitDefinition();
                    FastInvariant fastInvariant = new FastInvariant(getIdentifierSubstitutions(scm.getFastInvariant(), invariantUnit, surfaceClass));
                    fastSystem.addFastInvariant(fastInvariant);
                }
            }
            memSubDomain.setFastSystem(fastSystem);
            // constructor calls the 'refresh' method which constructs depemdency matrix, dependent/independent vars and pseudoconstants, etc.
            // FastSystemAnalyzer fs_analyzer =
            new FastSystemAnalyzer(fastSystem, mathDesc);
        }
        // 
        // create Membrane-region equations for potential of this resolved membrane
        // 
        Structure[] resolvedSurfaceStructures = membraneStructureAnalyzer.getStructures();
        for (int m = 0; m < resolvedSurfaceStructures.length; m++) {
            if (resolvedSurfaceStructures[m] instanceof Membrane) {
                Membrane membrane = (Membrane) resolvedSurfaceStructures[m];
                MembraneMapping membraneMapping = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(membrane);
                if (membraneMapping.getCalculateVoltage()) {
                    ElectricalDevice[] membraneDevices = potentialMapping.getElectricalDevices(membrane);
                    int numCapacitiveDevices = 0;
                    MembraneElectricalDevice capacitiveDevice = null;
                    for (int i = 0; i < membraneDevices.length; i++) {
                        if (membraneDevices[i] instanceof MembraneElectricalDevice) {
                            numCapacitiveDevices++;
                            capacitiveDevice = (MembraneElectricalDevice) membraneDevices[i];
                        }
                    }
                    if (numCapacitiveDevices != 1) {
                        throw new MappingException("expecting 1 capacitive electrical device on graph edge for membrane " + membrane.getName() + ", found '" + numCapacitiveDevices + "'");
                    }
                    if (mathDesc.getVariable(getMathSymbol(capacitiveDevice.getVoltageSymbol(), membraneMapping.getGeometryClass())) instanceof MembraneRegionVariable) {
                        MembraneRegionVariable vVar = (MembraneRegionVariable) mathDesc.getVariable(getMathSymbol(capacitiveDevice.getVoltageSymbol(), membraneMapping.getGeometryClass()));
                        Parameter initialVoltageParm = capacitiveDevice.getMembraneMapping().getInitialVoltageParameter();
                        Expression initExp = getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), capacitiveDevice.getMembraneMapping().getGeometryClass());
                        MembraneRegionEquation vEquation = new MembraneRegionEquation(vVar, initExp);
                        vEquation.setMembraneRateExpression(getIdentifierSubstitutions(potentialMapping.getOdeRHS(capacitiveDevice, this), membrane.getMembraneVoltage().getUnitDefinition().divideBy(timeUnit), capacitiveDevice.getMembraneMapping().getGeometryClass()));
                        memSubDomain.addEquation(vEquation);
                    }
                }
            }
        }
    }
    // create equations for event assignment or rate rule targets that are model params/species, etc.
    Set<VolVariable> hashKeySet = eventOrRateRuleVolVarHash.keySet();
    Iterator<VolVariable> volVarsIter = hashKeySet.iterator();
    // working under the assumption that we are dealing with non-spatial math, hence only one compartment domain!
    SubDomain subDomain = mathDesc.getSubDomains().nextElement();
    while (volVarsIter.hasNext()) {
        VolVariable volVar = volVarsIter.next();
        EventAssignmentOrRateRuleInitParameter initParam = eventOrRateRuleVolVarHash.get(volVar);
        // check event initial condition, it shouldn't contain vars, we have to do it here, coz we want to substitute functions...etc.
        Expression eapExp = MathUtilities.substituteFunctions(initParam.getExpression(), mathDesc);
        if (eapExp.getSymbols() != null) {
            for (String symbol : eapExp.getSymbols()) {
                SymbolTableEntry ste = eapExp.getSymbolBinding(symbol);
                if (ste instanceof VolVariable || ste instanceof MemVariable) {
                    throw new MathException("Variables are not allowed in Event assignment initial condition.\nEvent assignment target: " + volVar.getName() + " has variable (" + symbol + ") in its expression.");
                }
            }
        }
        Expression rateExpr = new Expression(0.0);
        RateRuleRateParameter rateParam = rateRuleRateParamHash.get(volVar);
        if (rateParam != null) {
            // this is a rate rule, get its expression.
            rateExpr = new Expression(getMathSymbol(rateParam, null));
        }
        Equation equation = new OdeEquation(volVar, new Expression(getMathSymbol(initParam, null)), rateExpr);
        subDomain.addEquation(equation);
    }
    // events - add events to math desc for event assignments that have parameters as target variables
    if (bioevents != null && bioevents.length > 0) {
        for (BioEvent be : bioevents) {
            // transform the bioEvent trigger/delay to math Event
            LocalParameter genTriggerParam = be.getParameter(BioEventParameterType.GeneralTriggerFunction);
            Expression mathTriggerExpr = getIdentifierSubstitutions(new Expression(genTriggerParam, be.getNameScope()), modelUnitSystem.getInstance_DIMENSIONLESS(), null);
            Delay mathDelay = null;
            LocalParameter delayParam = be.getParameter(BioEventParameterType.TriggerDelay);
            if (delayParam != null && delayParam.getExpression() != null && !delayParam.getExpression().compareEqual(new Expression(0.0))) {
                boolean bUseValsFromTriggerTime = be.getUseValuesFromTriggerTime();
                Expression mathDelayExpr = getIdentifierSubstitutions(new Expression(delayParam, be.getNameScope()), timeUnit, null);
                mathDelay = new Delay(bUseValsFromTriggerTime, mathDelayExpr);
            }
            // now deal with (bio)event Assignment translation to math EventAssignment
            ArrayList<EventAssignment> eventAssignments = be.getEventAssignments();
            ArrayList<Event.EventAssignment> mathEventAssignmentsList = new ArrayList<Event.EventAssignment>();
            if (eventAssignments != null) {
                for (EventAssignment ea : eventAssignments) {
                    SymbolTableEntry ste = simContext.getEntry(ea.getTarget().getName());
                    if (ste instanceof StructureSize) {
                        throw new RuntimeException("Event Assignment Variable for compartment size is not supported yet");
                    }
                    VCUnitDefinition eventAssignVarUnit = ste.getUnitDefinition();
                    Variable variable = varHash.getVariable(ste.getName());
                    Event.EventAssignment mathEA = new Event.EventAssignment(variable, getIdentifierSubstitutions(ea.getAssignmentExpression(), eventAssignVarUnit, null));
                    mathEventAssignmentsList.add(mathEA);
                }
            }
            // use the translated trigger, delay and event assignments to create (math) event
            Event mathEvent = new Event(be.getName(), mathTriggerExpr, mathDelay, mathEventAssignmentsList);
            mathDesc.addEvent(mathEvent);
        }
    }
    if (simContext.getMicroscopeMeasurement() != null && simContext.getMicroscopeMeasurement().getFluorescentSpecies().size() > 0) {
        MicroscopeMeasurement measurement = simContext.getMicroscopeMeasurement();
        Expression volumeConcExp = new Expression(0.0);
        Expression membraneDensityExp = new Expression(0.0);
        for (SpeciesContext speciesContext : measurement.getFluorescentSpecies()) {
            GeometryClass geometryClass = simContext.getGeometryContext().getStructureMapping(speciesContext.getStructure()).getGeometryClass();
            StructureMapping structureMapping = simContext.getGeometryContext().getStructureMapping(speciesContext.getStructure());
            StructureMappingParameter unitSizeParameter = structureMapping.getUnitSizeParameter();
            Expression mappedSpeciesContextExpression = Expression.mult(unitSizeParameter.getExpression(), new Expression(getMathSymbol(speciesContext, geometryClass)));
            VCUnitDefinition mappedSpeciesContextUnit = unitSizeParameter.getUnitDefinition().multiplyBy(speciesContext.getUnitDefinition());
            if (geometryClass instanceof SubVolume) {
                // volume function
                int dimension = 3;
                VCUnitDefinition desiredConcUnits = model.getUnitSystem().getInstance("molecules").divideBy(model.getUnitSystem().getLengthUnit().raiseTo(new ucar.units_vcell.RationalNumber(dimension)));
                Expression unitFactor = getUnitFactor(desiredConcUnits.divideBy(mappedSpeciesContextUnit));
                volumeConcExp = Expression.add(volumeConcExp, Expression.mult(unitFactor, mappedSpeciesContextExpression)).flatten();
            } else if (geometryClass instanceof SurfaceClass) {
                // membrane function
                int dimension = 2;
                VCUnitDefinition desiredSurfaceDensityUnits = model.getUnitSystem().getInstance("molecules").divideBy(model.getUnitSystem().getLengthUnit().raiseTo(new ucar.units_vcell.RationalNumber(dimension)));
                Expression unitFactor = getUnitFactor(desiredSurfaceDensityUnits.divideBy(mappedSpeciesContextUnit));
                membraneDensityExp = Expression.add(membraneDensityExp, Expression.mult(unitFactor, mappedSpeciesContextExpression)).flatten();
            } else {
                throw new MathException("unsupported geometry mapping for microscopy measurement");
            }
        }
        ConvolutionKernel kernel = measurement.getConvolutionKernel();
        if (kernel instanceof ExperimentalPSF) {
            if (!membraneDensityExp.isZero()) {
                throw new MappingException("membrane variables and functions not yet supported for Z projection in Microcopy Measurements");
            }
            ExperimentalPSF psf = (ExperimentalPSF) kernel;
            DataSymbol psfDataSymbol = psf.getPSFDataSymbol();
            if (psfDataSymbol instanceof FieldDataSymbol) {
                FieldDataSymbol fieldDataSymbol = (FieldDataSymbol) psfDataSymbol;
                String fieldDataName = ((FieldDataSymbol) psfDataSymbol).getExternalDataIdentifier().getName();
                Expression psfExp = Expression.function(FieldFunctionDefinition.FUNCTION_name, new Expression("'" + fieldDataName + "'"), new Expression("'" + fieldDataSymbol.getFieldDataVarName() + "'"), new Expression(fieldDataSymbol.getFieldDataVarTime()), new Expression("'" + fieldDataSymbol.getFieldDataVarType() + "'"));
                varHash.addVariable(new Function("__PSF__", psfExp, null));
            }
            Expression convExp = Expression.function(ConvFunctionDefinition.FUNCTION_name, volumeConcExp, new Expression("__PSF__"));
            varHash.addVariable(newFunctionOrConstant(measurement.getName(), convExp, null));
        } else if (kernel instanceof GaussianConvolutionKernel) {
            GaussianConvolutionKernel gaussianConvolutionKernel = (GaussianConvolutionKernel) kernel;
            GaussianConvolutionDataGeneratorKernel mathKernel = new GaussianConvolutionDataGeneratorKernel(gaussianConvolutionKernel.getSigmaXY_um(), gaussianConvolutionKernel.getSigmaZ_um());
            ConvolutionDataGenerator dataGenerator = new ConvolutionDataGenerator(measurement.getName(), mathKernel, volumeConcExp, membraneDensityExp);
            mathDesc.getPostProcessingBlock().addDataGenerator(dataGenerator);
        } else if (kernel instanceof ProjectionZKernel) {
            if (mathDesc.getGeometry().getDimension() == 3) {
                if (!membraneDensityExp.isZero()) {
                    throw new MappingException("membrane variables and functions not yet supported for Z projection in Microcopy Measurements");
                }
                ProjectionDataGenerator dataGenerator = new ProjectionDataGenerator(measurement.getName(), null, ProjectionDataGenerator.Axis.z, ProjectionDataGenerator.Operation.sum, volumeConcExp);
                mathDesc.getPostProcessingBlock().addDataGenerator(dataGenerator);
            } else {
                throw new MappingException("Z Projection is only supported in 3D spatial applications.");
            }
        }
    }
    // 
    for (int i = 0; i < fieldMathMappingParameters.length; i++) {
        if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
            GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
            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());
    }
// System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string begin ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
// System.out.println(mathDesc.getVCML());
// System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string end ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
}
Also used : MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) ArrayList(java.util.ArrayList) StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter) SpeciesContext(cbit.vcell.model.SpeciesContext) Feature(cbit.vcell.model.Feature) MemVariable(cbit.vcell.math.MemVariable) SubVolume(cbit.vcell.geometry.SubVolume) Vector(java.util.Vector) SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter) FastInvariant(cbit.vcell.math.FastInvariant) GaussianConvolutionKernel(cbit.vcell.mapping.MicroscopeMeasurement.GaussianConvolutionKernel) PropertyVetoException(java.beans.PropertyVetoException) FieldDataSymbol(cbit.vcell.data.FieldDataSymbol) DataSymbol(cbit.vcell.data.DataSymbol) VCUnitDefinition(cbit.vcell.units.VCUnitDefinition) FastSystem(cbit.vcell.math.FastSystem) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) ReactionStep(cbit.vcell.model.ReactionStep) MembraneRegionEquation(cbit.vcell.math.MembraneRegionEquation) SurfaceClass(cbit.vcell.geometry.SurfaceClass) VariableHash(cbit.vcell.math.VariableHash) ExperimentalPSF(cbit.vcell.mapping.MicroscopeMeasurement.ExperimentalPSF) ConvolutionDataGenerator(cbit.vcell.math.ConvolutionDataGenerator) GaussianConvolutionDataGeneratorKernel(cbit.vcell.math.ConvolutionDataGenerator.GaussianConvolutionDataGeneratorKernel) Structure(cbit.vcell.model.Structure) ModelUnitSystem(cbit.vcell.model.ModelUnitSystem) LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter) VoltageClampElectricalDevice(cbit.vcell.mapping.potential.VoltageClampElectricalDevice) StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter) Expression(cbit.vcell.parser.Expression) Model(cbit.vcell.model.Model) StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter) Parameter(cbit.vcell.model.Parameter) SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter) SimulationContextParameter(cbit.vcell.mapping.SimulationContext.SimulationContextParameter) ModelParameter(cbit.vcell.model.Model.ModelParameter) SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter) LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter) Event(cbit.vcell.math.Event) ProjectionDataGenerator(cbit.vcell.math.ProjectionDataGenerator) SurfaceRegionObject(cbit.vcell.mapping.spatial.SurfaceRegionObject) FieldDataSymbol(cbit.vcell.data.FieldDataSymbol) MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable) HashMap(java.util.HashMap) MathDescription(cbit.vcell.math.MathDescription) MembraneElectricalDevice(cbit.vcell.mapping.potential.MembraneElectricalDevice) Delay(cbit.vcell.math.Event.Delay) PropertyVetoException(java.beans.PropertyVetoException) PdeEquation(cbit.vcell.math.PdeEquation) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) PointSubDomain(cbit.vcell.math.PointSubDomain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) CurrentClampElectricalDevice(cbit.vcell.mapping.potential.CurrentClampElectricalDevice) SpatialQuantity(cbit.vcell.mapping.spatial.SpatialObject.SpatialQuantity) SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter) EventAssignment(cbit.vcell.mapping.BioEvent.EventAssignment) CurrentClampElectricalDevice(cbit.vcell.mapping.potential.CurrentClampElectricalDevice) MembraneElectricalDevice(cbit.vcell.mapping.potential.MembraneElectricalDevice) ElectricalDevice(cbit.vcell.mapping.potential.ElectricalDevice) VoltageClampElectricalDevice(cbit.vcell.mapping.potential.VoltageClampElectricalDevice) VolVariable(cbit.vcell.math.VolVariable) StructureSize(cbit.vcell.model.Structure.StructureSize) ProjectionZKernel(cbit.vcell.mapping.MicroscopeMeasurement.ProjectionZKernel) ModelParameter(cbit.vcell.model.Model.ModelParameter) OdeEquation(cbit.vcell.math.OdeEquation) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) PointSubDomain(cbit.vcell.math.PointSubDomain) Domain(cbit.vcell.math.Variable.Domain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) JumpCondition(cbit.vcell.math.JumpCondition) GeometryClass(cbit.vcell.geometry.GeometryClass) VolVariable(cbit.vcell.math.VolVariable) MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable) PointVariable(cbit.vcell.math.PointVariable) VolumeRegionVariable(cbit.vcell.math.VolumeRegionVariable) MemVariable(cbit.vcell.math.MemVariable) Variable(cbit.vcell.math.Variable) SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter) Constant(cbit.vcell.math.Constant) VolumeRegionVariable(cbit.vcell.math.VolumeRegionVariable) Function(cbit.vcell.math.Function) SymbolTableEntry(cbit.vcell.parser.SymbolTableEntry) Membrane(cbit.vcell.model.Membrane) VolumeRegionEquation(cbit.vcell.math.VolumeRegionEquation) PotentialMapping(cbit.vcell.mapping.potential.PotentialMapping) EventAssignment(cbit.vcell.mapping.BioEvent.EventAssignment) FieldFunctionArguments(cbit.vcell.field.FieldFunctionArguments) PdeEquation(cbit.vcell.math.PdeEquation) ComputeMembraneMetricEquation(cbit.vcell.math.ComputeMembraneMetricEquation) VolumeRegionEquation(cbit.vcell.math.VolumeRegionEquation) OdeEquation(cbit.vcell.math.OdeEquation) MembraneRegionEquation(cbit.vcell.math.MembraneRegionEquation) Equation(cbit.vcell.math.Equation) FastRate(cbit.vcell.math.FastRate) SimulationContextParameter(cbit.vcell.mapping.SimulationContext.SimulationContextParameter) ConvolutionKernel(cbit.vcell.mapping.MicroscopeMeasurement.ConvolutionKernel) GaussianConvolutionKernel(cbit.vcell.mapping.MicroscopeMeasurement.GaussianConvolutionKernel) MathException(cbit.vcell.math.MathException)

Example 53 with Structure

use of cbit.vcell.model.Structure in project vcell by virtualcell.

the class GeometryContext method refreshStructureMappings.

/**
 * This method was created by a SmartGuide.
 * @throws ExpressionBindingException
 */
public void refreshStructureMappings() throws MappingException, PropertyVetoException {
    // 
    // step through all structure mappings
    // 
    StructureMapping[] newStructureMappings = (StructureMapping[]) fieldStructureMappings.clone();
    // fill in geometryClass for those MembraneMappings that were not mapped explicitly (and topology is available)
    for (int j = 0; j < newStructureMappings.length; j++) {
        StructureMapping structureMapping = newStructureMappings[j];
        if (structureMapping instanceof MembraneMapping && structureMapping.getGeometryClass() == null) {
            MembraneMapping membraneMapping = (MembraneMapping) structureMapping;
            Feature insideFeature = getModel().getStructureTopology().getInsideFeature(membraneMapping.getMembrane());
            Feature outsideFeature = getModel().getStructureTopology().getOutsideFeature(membraneMapping.getMembrane());
            if (insideFeature != null && outsideFeature != null) {
                FeatureMapping insideFeatureMapping = (FeatureMapping) getStructureMapping(insideFeature);
                FeatureMapping outsideFeatureMapping = (FeatureMapping) getStructureMapping(outsideFeature);
                if (insideFeatureMapping.getGeometryClass() == outsideFeatureMapping.getGeometryClass()) {
                    membraneMapping.setGeometryClass(insideFeatureMapping.getGeometryClass());
                } else if (insideFeatureMapping.getGeometryClass() instanceof SubVolume && outsideFeatureMapping.getGeometryClass() instanceof SubVolume) {
                    SubVolume insideSubVolume = (SubVolume) insideFeatureMapping.getGeometryClass();
                    SubVolume outsideSubVolume = (SubVolume) outsideFeatureMapping.getGeometryClass();
                    SurfaceClass surfaceClass = getGeometry().getGeometrySurfaceDescription().getSurfaceClass(insideSubVolume, outsideSubVolume);
                    if (surfaceClass != null) {
                        membraneMapping.setGeometryClass(surfaceClass);
                    }
                }
            }
        }
    }
    for (int j = 0; j < newStructureMappings.length; j++) {
        StructureMapping structureMapping = newStructureMappings[j];
        Structure mappedStructure = structureMapping.getStructure();
        // SubVolume mappedSubvolume = structureMapping.getSubVolume();
        Structure newStructure = null;
        GeometryClass newGeometryClass = null;
        boolean structureFound = false;
        boolean geometryClassFound = false;
        // 
        // match up with structures defined within model
        // 
        Structure[] modelStructures = getModel().getStructures();
        for (int i = 0; i < modelStructures.length; i++) {
            Structure modelStructure = modelStructures[i];
            if (modelStructure.compareEqual(mappedStructure)) {
                structureFound = true;
                newStructure = modelStructure;
                break;
            }
        }
        // 
        // match up with geometryClasses defined within geometry
        // 
        GeometryClass[] geometryClasses = getGeometry().getGeometryClasses();
        for (int i = 0; i < geometryClasses.length; i++) {
            if (geometryClasses[i].compareEqual(structureMapping.getGeometryClass())) {
                geometryClassFound = true;
                newGeometryClass = geometryClasses[i];
                break;
            }
        }
        // 
        if (!(structureFound && geometryClassFound)) {
            newStructureMappings = (StructureMapping[]) BeanUtils.removeElement(newStructureMappings, structureMapping);
            j--;
        // //
        // // delete accompanied membrane mapping if exists
        // //
        // for (int i = 0; i < newStructureMappings.length; i++){
        // if (newStructureMappings[i] instanceof MembraneMapping){
        // MembraneMapping membraneMapping = (MembraneMapping)newStructureMappings[i];
        // if (membraneMapping.getMembrane()==null ||
        // membraneMapping.getMembrane().getInsideFeature() == structureMapping.getStructure() ||
        // membraneMapping.getMembrane().getOutsideFeature() == structureMapping.getStructure()){
        // newStructureMappings = (StructureMapping[])BeanUtils.removeElement(newStructureMappings,membraneMapping);
        // break;
        // }
        // }
        // }
        } else {
            // update references to Structure and SubVolume to correspond to those of Model and Geometry
            structureMapping.setGeometryClass(newGeometryClass);
        }
        if (structureFound) {
            structureMapping.setStructure(newStructure);
        }
    }
    // 
    // add default mappings for any new structures
    // 
    Structure[] structures = getModel().getStructures();
    for (int i = 0; i < structures.length; i++) {
        Structure structure = structures[i];
        StructureMapping sm = null;
        for (int j = 0; j < newStructureMappings.length; j++) {
            if (newStructureMappings[j].getStructure().compareEqual(structure)) {
                sm = newStructureMappings[j];
            }
        }
        if (sm == null) {
            if (structure instanceof Feature) {
                FeatureMapping fm = new FeatureMapping((Feature) structure, fieldSimulationContext, getModel().getUnitSystem());
                fm.setSimulationContext(this.fieldSimulationContext);
                newStructureMappings = (StructureMapping[]) BeanUtils.addElement(newStructureMappings, fm);
                if (getGeometry().getDimension() == 0) {
                    fm.setGeometryClass((CompartmentSubVolume) getGeometry().getGeometrySpec().getSubVolumes()[0]);
                }
            } else if (structure instanceof Membrane) {
                MembraneMapping mm = new MembraneMapping((Membrane) structure, fieldSimulationContext, getModel().getUnitSystem());
                mm.setSimulationContext(fieldSimulationContext);
                newStructureMappings = (StructureMapping[]) BeanUtils.addElement(newStructureMappings, mm);
                if (getGeometry().getDimension() == 0) {
                    mm.setGeometryClass((CompartmentSubVolume) getGeometry().getGeometrySpec().getSubVolumes()[0]);
                }
            } else {
                throw new MappingException("unsupported Structure Mapping for structure " + structure.getClass().toString());
            }
        }
    }
    if (newStructureMappings != fieldStructureMappings) {
        try {
            setStructureMappings(newStructureMappings);
        } catch (Exception e) {
            e.printStackTrace(System.out);
            throw new MappingException(e.getMessage());
        }
    }
    fixMembraneMappings();
}
Also used : GeometryClass(cbit.vcell.geometry.GeometryClass) SurfaceClass(cbit.vcell.geometry.SurfaceClass) Feature(cbit.vcell.model.Feature) PropertyVetoException(java.beans.PropertyVetoException) ExpressionException(cbit.vcell.parser.ExpressionException) ExpressionBindingException(cbit.vcell.parser.ExpressionBindingException) CompartmentSubVolume(cbit.vcell.geometry.CompartmentSubVolume) SubVolume(cbit.vcell.geometry.SubVolume) CompartmentSubVolume(cbit.vcell.geometry.CompartmentSubVolume) Membrane(cbit.vcell.model.Membrane) Structure(cbit.vcell.model.Structure)

Example 54 with Structure

use of cbit.vcell.model.Structure in project vcell by virtualcell.

the class GeometryContext method enforceHierarchicalBoundaryConditions.

public void enforceHierarchicalBoundaryConditions(StructureTopology structureTopology) {
    if (structureTopology != null) {
        for (StructureMapping sm : fieldStructureMappings) {
            // 
            // look for top level parent structure which is mapped to the same geometric domain,
            // 
            Structure topStructureForDomain = sm.getStructure();
            while (true) {
                Structure parent = structureTopology.getParentStructure(topStructureForDomain);
                if (parent == null || getStructureMapping(parent).getGeometryClass() != sm.getGeometryClass()) {
                    break;
                } else {
                    topStructureForDomain = parent;
                }
            }
            // 
            if (topStructureForDomain != sm.getStructure()) {
                StructureMapping parentSM = getStructureMapping(topStructureForDomain);
                sm.setBoundaryConditionTypeXm(parentSM.getBoundaryConditionTypeXm());
                sm.setBoundaryConditionTypeXp(parentSM.getBoundaryConditionTypeXp());
                sm.setBoundaryConditionTypeYm(parentSM.getBoundaryConditionTypeYm());
                sm.setBoundaryConditionTypeYp(parentSM.getBoundaryConditionTypeYp());
                sm.setBoundaryConditionTypeZm(parentSM.getBoundaryConditionTypeZm());
                sm.setBoundaryConditionTypeZp(parentSM.getBoundaryConditionTypeZp());
            }
        }
    }
}
Also used : Structure(cbit.vcell.model.Structure)

Example 55 with Structure

use of cbit.vcell.model.Structure in project vcell by virtualcell.

the class NetworkTransformer method transform.

private void transform(SimulationContext simContext, SimulationContext transformedSimulationContext, ArrayList<ModelEntityMapping> entityMappings, MathMappingCallback mathMappingCallback, NetworkGenerationRequirements networkGenerationRequirements) {
    String msg = "Generating network: flattening...";
    mathMappingCallback.setMessage(msg);
    TaskCallbackMessage tcm = new TaskCallbackMessage(TaskCallbackStatus.Clean, "");
    simContext.appendToConsole(tcm);
    tcm = new TaskCallbackMessage(TaskCallbackStatus.TaskStart, msg);
    simContext.appendToConsole(tcm);
    long startTime = System.currentTimeMillis();
    System.out.println("Convert to bngl, execute BNG, retrieve the results.");
    try {
        BNGOutputSpec outputSpec = generateNetwork(simContext, mathMappingCallback, networkGenerationRequirements);
        if (mathMappingCallback.isInterrupted()) {
            msg = "Canceled by user.";
            tcm = new TaskCallbackMessage(TaskCallbackStatus.Error, msg);
            simContext.appendToConsole(tcm);
            throw new UserCancelException(msg);
        }
        long endTime = System.currentTimeMillis();
        long elapsedTime = endTime - startTime;
        System.out.println("     " + elapsedTime + " milliseconds");
        Model model = transformedSimulationContext.getModel();
        ReactionContext reactionContext = transformedSimulationContext.getReactionContext();
        // ---- Parameters -----------------------------------------------------------------------------------------------
        startTime = System.currentTimeMillis();
        for (int i = 0; i < outputSpec.getBNGParams().length; i++) {
            BNGParameter p = outputSpec.getBNGParams()[i];
            // System.out.println(i+1 + ":\t\t"+ p.toString());
            if (model.getRbmModelContainer().getParameter(p.getName()) != null) {
                // if it's already there we don't try to add it again; this should be true for all of them!
                continue;
            }
            String s = p.getName();
            FakeSeedSpeciesInitialConditionsParameter fakeICParam = FakeSeedSpeciesInitialConditionsParameter.fromString(s);
            if (speciesEquivalenceMap.containsKey(fakeICParam)) {
                // we get rid of the fake parameters we use as keys
                continue;
            }
            FakeReactionRuleRateParameter fakeKineticParam = FakeReactionRuleRateParameter.fromString(s);
            if (fakeKineticParam != null) {
                System.out.println("found fakeKineticParam " + fakeKineticParam.fakeParameterName);
                // we get rid of the fake parameters we use as keys
                continue;
            }
            throw new RuntimeException("unexpected parameter " + p.getName() + " in internal BNG processing");
        // Expression exp = new Expression(p.getValue());
        // exp.bindExpression(model.getRbmModelContainer().getSymbolTable());
        // model.getRbmModelContainer().addParameter(p.getName(), exp, model.getUnitSystem().getInstance_TBD());
        }
        endTime = System.currentTimeMillis();
        elapsedTime = endTime - startTime;
        msg = "Adding " + outputSpec.getBNGParams().length + " parameters to model, " + elapsedTime + " ms";
        System.out.println(msg);
        // ---- Species ------------------------------------------------------------------------------------------------------------
        mathMappingCallback.setMessage("generating network: adding species...");
        mathMappingCallback.setProgressFraction(progressFractionQuota / 4.0f);
        startTime = System.currentTimeMillis();
        System.out.println("\nSpecies :");
        // the reactions will need this map to recover the names of species knowing only the networkFileIndex
        HashMap<Integer, String> speciesMap = new HashMap<Integer, String>();
        LinkedHashMap<String, Species> sMap = new LinkedHashMap<String, Species>();
        LinkedHashMap<String, SpeciesContext> scMap = new LinkedHashMap<String, SpeciesContext>();
        LinkedHashMap<String, BNGSpecies> crossMap = new LinkedHashMap<String, BNGSpecies>();
        List<SpeciesContext> noMapForThese = new ArrayList<SpeciesContext>();
        // final int decimalTickCount = Math.max(outputSpec.getBNGSpecies().length/10, 1);
        for (int i = 0; i < outputSpec.getBNGSpecies().length; i++) {
            BNGSpecies s = outputSpec.getBNGSpecies()[i];
            // System.out.println(i+1 + ":\t\t"+ s.toString());
            String key = s.getConcentration().infix();
            FakeSeedSpeciesInitialConditionsParameter fakeParam = FakeSeedSpeciesInitialConditionsParameter.fromString(key);
            if (fakeParam != null) {
                Pair<SpeciesContext, Expression> value = speciesEquivalenceMap.get(fakeParam);
                // the species context of the original model
                SpeciesContext originalsc = value.one;
                Expression initial = value.two;
                // replace the fake initial condition with the real one
                s.setConcentration(initial);
                // we'll have to find the species context from the cloned model which correspond to the original species
                SpeciesContext sc = model.getSpeciesContext(originalsc.getName());
                // System.out.println(sc.getName() + ", " + sc.getSpecies().getCommonName() + "   ...is one of the original seed species.");
                // existing name
                speciesMap.put(s.getNetworkFileIndex(), sc.getName());
                sMap.put(sc.getName(), sc.getSpecies());
                scMap.put(sc.getName(), sc);
                crossMap.put(sc.getName(), s);
                noMapForThese.add(sc);
                continue;
            }
            // all these species are new!
            // generate unique name for the species
            int count = 0;
            String speciesName = null;
            String nameRoot = "s";
            String speciesPatternNameString = s.extractName();
            while (true) {
                speciesName = nameRoot + count;
                if (Model.isNameUnused(speciesName, model) && !sMap.containsKey(speciesName) && !scMap.containsKey(speciesName)) {
                    break;
                }
                count++;
            }
            // newly created name
            speciesMap.put(s.getNetworkFileIndex(), speciesName);
            SpeciesContext speciesContext;
            if (s.hasCompartment()) {
                String speciesPatternCompartmentString = s.extractCompartment();
                speciesContext = new SpeciesContext(new Species(speciesName, s.getName()), model.getStructure(speciesPatternCompartmentString), null);
            } else {
                speciesContext = new SpeciesContext(new Species(speciesName, s.getName()), model.getStructure(0), null);
            }
            speciesContext.setName(speciesName);
            try {
                if (speciesPatternNameString != null) {
                    SpeciesPattern sp = RbmUtils.parseSpeciesPattern(speciesPatternNameString, model);
                    speciesContext.setSpeciesPattern(sp);
                }
            } catch (ParseException e) {
                e.printStackTrace();
                throw new RuntimeException("Bad format for species pattern string: " + e.getMessage());
            }
            // speciesContext.setSpeciesPatternString(speciesPatternString);
            // model.addSpecies(speciesContext.getSpecies());
            // model.addSpeciesContext(speciesContext);
            sMap.put(speciesName, speciesContext.getSpecies());
            scMap.put(speciesName, speciesContext);
            crossMap.put(speciesName, s);
            // }
            if (mathMappingCallback.isInterrupted()) {
                msg = "Canceled by user.";
                tcm = new TaskCallbackMessage(TaskCallbackStatus.Error, msg);
                simContext.appendToConsole(tcm);
                throw new UserCancelException(msg);
            }
        // if(i%50 == 0) {
        // System.out.println(i+"");
        // }
        // if(i%decimalTickCount == 0) {
        // int multiplier = i/decimalTickCount;
        // float progress = progressFractionQuota/4.0f + progressFractionQuotaSpecies*multiplier;
        // mathMappingCallback.setProgressFraction(progress);
        // }
        }
        for (SpeciesContext sc1 : model.getSpeciesContexts()) {
            boolean found = false;
            for (Map.Entry<String, SpeciesContext> entry : scMap.entrySet()) {
                SpeciesContext sc2 = entry.getValue();
                if (sc1.getName().equals(sc2.getName())) {
                    found = true;
                    // System.out.println("found species context " + sc1.getName() + " of species " + sc1.getSpecies().getCommonName() + " // " + sc2.getSpecies().getCommonName());
                    break;
                }
            }
            if (found == false) {
                // we add to the map the species context and the species which exist in the model but which are not in the map yet
                // the only ones in this situation should be plain species which were not given to bngl for flattening (they are flat already)
                // System.out.println("species context " + sc1.getName() + " not found in the map. Adding it.");
                scMap.put(sc1.getName(), sc1);
                sMap.put(sc1.getName(), sc1.getSpecies());
                noMapForThese.add(sc1);
            }
        }
        for (Species s1 : model.getSpecies()) {
            boolean found = false;
            for (Map.Entry<String, Species> entry : sMap.entrySet()) {
                Species s2 = entry.getValue();
                if (s1.getCommonName().equals(s2.getCommonName())) {
                    found = true;
                    // System.out.println("found species " + s1.getCommonName());
                    break;
                }
            }
            if (found == false) {
                System.err.println("species " + s1.getCommonName() + " not found in the map!");
            }
        }
        SpeciesContext[] sca = new SpeciesContext[scMap.size()];
        scMap.values().toArray(sca);
        Species[] sa = new HashSet<Species>(sMap.values()).toArray(new Species[0]);
        model.setSpecies(sa);
        model.setSpeciesContexts(sca);
        boolean isSpatial = transformedSimulationContext.getGeometry().getDimension() > 0;
        for (SpeciesContext sc : sca) {
            if (noMapForThese.contains(sc)) {
                continue;
            }
            SpeciesContextSpec scs = reactionContext.getSpeciesContextSpec(sc);
            Parameter param = scs.getParameter(SpeciesContextSpec.ROLE_InitialConcentration);
            BNGSpecies s = crossMap.get(sc.getName());
            param.setExpression(s.getConcentration());
            SpeciesContext origSpeciesContext = simContext.getModel().getSpeciesContext(s.getName());
            if (origSpeciesContext != null) {
                ModelEntityMapping em = new ModelEntityMapping(origSpeciesContext, sc);
                entityMappings.add(em);
            } else {
                ModelEntityMapping em = new ModelEntityMapping(new GeneratedSpeciesSymbolTableEntry(sc), sc);
                if (isSpatial) {
                    scs.initializeForSpatial();
                }
                entityMappings.add(em);
            }
        }
        // for(SpeciesContext sc : sca) {		// clean all the species patterns from the flattened species, we have no sp now
        // sc.setSpeciesPattern(null);
        // }
        endTime = System.currentTimeMillis();
        elapsedTime = endTime - startTime;
        msg = "Adding " + outputSpec.getBNGSpecies().length + " species to model, " + elapsedTime + " ms";
        System.out.println(msg);
        // ---- Reactions -----------------------------------------------------------------------------------------------------
        mathMappingCallback.setMessage("generating network: adding reactions...");
        mathMappingCallback.setProgressFraction(progressFractionQuota / 4.0f * 3.0f);
        startTime = System.currentTimeMillis();
        System.out.println("\nReactions :");
        Map<String, HashSet<String>> ruleKeyMap = new HashMap<String, HashSet<String>>();
        Map<String, BNGReaction> directBNGReactionsMap = new HashMap<String, BNGReaction>();
        Map<String, BNGReaction> reverseBNGReactionsMap = new HashMap<String, BNGReaction>();
        for (int i = 0; i < outputSpec.getBNGReactions().length; i++) {
            BNGReaction r = outputSpec.getBNGReactions()[i];
            if (!r.isRuleReversed()) {
                // direct
                directBNGReactionsMap.put(r.getKey(), r);
            } else {
                reverseBNGReactionsMap.put(r.getKey(), r);
            }
            // 
            // for each rule name, store set of keySets (number of unique keysets are number of generated reactions from this ruleName).
            // 
            HashSet<String> keySet = ruleKeyMap.get(r.getRuleName());
            if (keySet == null) {
                keySet = new HashSet<String>();
                ruleKeyMap.put(r.getRuleName(), keySet);
            }
            keySet.add(r.getKey());
        }
        Map<String, ReactionStep> reactionStepMap = new HashMap<String, ReactionStep>();
        for (int i = 0; i < outputSpec.getBNGReactions().length; i++) {
            BNGReaction bngReaction = outputSpec.getBNGReactions()[i];
            // System.out.println(i+1 + ":\t\t"+ r.writeReaction());
            String baseName = bngReaction.getRuleName();
            String reactionName = null;
            HashSet<String> keySetsForThisRule = ruleKeyMap.get(bngReaction.getRuleName());
            if (keySetsForThisRule.size() == 1 && model.getReactionStep(bngReaction.getRuleName()) == null && !reactionStepMap.containsKey(bngReaction.getRuleName())) {
                // we can reuse the reaction rule labels
                reactionName = bngReaction.getRuleName();
            } else {
                reactionName = bngReaction.getRuleName() + "_0";
                while (true) {
                    if (model.getReactionStep(reactionName) == null && !reactionStepMap.containsKey(reactionName)) {
                        // we can reuse the reaction rule labels
                        break;
                    }
                    reactionName = TokenMangler.getNextEnumeratedToken(reactionName);
                }
            }
            // 
            if (directBNGReactionsMap.containsValue(bngReaction)) {
                BNGReaction forwardBNGReaction = bngReaction;
                BNGReaction reverseBNGReaction = reverseBNGReactionsMap.get(bngReaction.getKey());
                String name = forwardBNGReaction.getRuleName();
                if (name.endsWith(ReactionRule.DirectHalf)) {
                    name = name.substring(0, name.indexOf(ReactionRule.DirectHalf));
                }
                if (name.endsWith(ReactionRule.InverseHalf)) {
                    name = name.substring(0, name.indexOf(ReactionRule.InverseHalf));
                }
                ReactionRule rr = model.getRbmModelContainer().getReactionRule(name);
                Structure structure = rr.getStructure();
                boolean bReversible = reverseBNGReaction != null;
                SimpleReaction sr = new SimpleReaction(model, structure, reactionName, bReversible);
                for (int j = 0; j < forwardBNGReaction.getReactants().length; j++) {
                    BNGSpecies s = forwardBNGReaction.getReactants()[j];
                    String scName = speciesMap.get(s.getNetworkFileIndex());
                    SpeciesContext sc = model.getSpeciesContext(scName);
                    Reactant reactant = sr.getReactant(scName);
                    if (reactant == null) {
                        int stoichiometry = 1;
                        sr.addReactant(sc, stoichiometry);
                    } else {
                        int stoichiometry = reactant.getStoichiometry();
                        stoichiometry += 1;
                        reactant.setStoichiometry(stoichiometry);
                    }
                }
                for (int j = 0; j < forwardBNGReaction.getProducts().length; j++) {
                    BNGSpecies s = forwardBNGReaction.getProducts()[j];
                    String scName = speciesMap.get(s.getNetworkFileIndex());
                    SpeciesContext sc = model.getSpeciesContext(scName);
                    Product product = sr.getProduct(scName);
                    if (product == null) {
                        int stoichiometry = 1;
                        sr.addProduct(sc, stoichiometry);
                    } else {
                        int stoichiometry = product.getStoichiometry();
                        stoichiometry += 1;
                        product.setStoichiometry(stoichiometry);
                    }
                }
                MassActionKinetics targetKinetics = new MassActionKinetics(sr);
                sr.setKinetics(targetKinetics);
                KineticsParameter kforward = targetKinetics.getForwardRateParameter();
                KineticsParameter kreverse = targetKinetics.getReverseRateParameter();
                String kforwardNewName = rr.getKineticLaw().getLocalParameter(RbmKineticLawParameterType.MassActionForwardRate).getName();
                if (!kforward.getName().equals(kforwardNewName)) {
                    targetKinetics.renameParameter(kforward.getName(), kforwardNewName);
                    kforward = targetKinetics.getForwardRateParameter();
                }
                final String kreverseNewName = rr.getKineticLaw().getLocalParameter(RbmKineticLawParameterType.MassActionReverseRate).getName();
                if (!kreverse.getName().equals(kreverseNewName)) {
                    targetKinetics.renameParameter(kreverse.getName(), kreverseNewName);
                    kreverse = targetKinetics.getReverseRateParameter();
                }
                applyKineticsExpressions(forwardBNGReaction, kforward, targetKinetics);
                if (reverseBNGReaction != null) {
                    applyKineticsExpressions(reverseBNGReaction, kreverse, targetKinetics);
                }
                // String fieldParameterName = kforward.getName();
                // fieldParameterName += "_" + r.getRuleName();
                // kforward.setName(fieldParameterName);
                reactionStepMap.put(reactionName, sr);
            } else if (reverseBNGReactionsMap.containsValue(bngReaction) && !directBNGReactionsMap.containsKey(bngReaction.getKey())) {
                // reverse only (must be irreversible)
                BNGReaction reverseBNGReaction = reverseBNGReactionsMap.get(bngReaction.getKey());
                ReactionRule rr = model.getRbmModelContainer().getReactionRule(reverseBNGReaction.extractRuleName());
                Structure structure = rr.getStructure();
                boolean bReversible = false;
                SimpleReaction sr = new SimpleReaction(model, structure, reactionName, bReversible);
                for (int j = 0; j < reverseBNGReaction.getReactants().length; j++) {
                    BNGSpecies s = reverseBNGReaction.getReactants()[j];
                    String scName = speciesMap.get(s.getNetworkFileIndex());
                    SpeciesContext sc = model.getSpeciesContext(scName);
                    Reactant reactant = sr.getReactant(scName);
                    if (reactant == null) {
                        int stoichiometry = 1;
                        sr.addReactant(sc, stoichiometry);
                    } else {
                        int stoichiometry = reactant.getStoichiometry();
                        stoichiometry += 1;
                        reactant.setStoichiometry(stoichiometry);
                    }
                }
                for (int j = 0; j < reverseBNGReaction.getProducts().length; j++) {
                    BNGSpecies s = reverseBNGReaction.getProducts()[j];
                    String scName = speciesMap.get(s.getNetworkFileIndex());
                    SpeciesContext sc = model.getSpeciesContext(scName);
                    Product product = sr.getProduct(scName);
                    if (product == null) {
                        int stoichiometry = 1;
                        sr.addProduct(sc, stoichiometry);
                    } else {
                        int stoichiometry = product.getStoichiometry();
                        stoichiometry += 1;
                        product.setStoichiometry(stoichiometry);
                    }
                }
                MassActionKinetics k = new MassActionKinetics(sr);
                sr.setKinetics(k);
                KineticsParameter kforward = k.getForwardRateParameter();
                KineticsParameter kreverse = k.getReverseRateParameter();
                String kforwardNewName = rr.getKineticLaw().getLocalParameter(RbmKineticLawParameterType.MassActionForwardRate).getName();
                if (!kforward.getName().equals(kforwardNewName)) {
                    k.renameParameter(kforward.getName(), kforwardNewName);
                    kforward = k.getForwardRateParameter();
                }
                final String kreverseNewName = rr.getKineticLaw().getLocalParameter(RbmKineticLawParameterType.MassActionReverseRate).getName();
                if (!kreverse.getName().equals(kreverseNewName)) {
                    k.renameParameter(kreverse.getName(), kreverseNewName);
                    kreverse = k.getReverseRateParameter();
                }
                applyKineticsExpressions(reverseBNGReaction, kforward, k);
                // String fieldParameterName = kforward.getName();
                // fieldParameterName += "_" + r.getRuleName();
                // kforward.setName(fieldParameterName);
                reactionStepMap.put(reactionName, sr);
            }
        }
        for (ReactionStep rs : model.getReactionSteps()) {
            reactionStepMap.put(rs.getName(), rs);
        }
        ReactionStep[] reactionSteps = new ReactionStep[reactionStepMap.size()];
        reactionStepMap.values().toArray(reactionSteps);
        model.setReactionSteps(reactionSteps);
        if (mathMappingCallback.isInterrupted()) {
            msg = "Canceled by user.";
            tcm = new TaskCallbackMessage(TaskCallbackStatus.Error, msg);
            simContext.appendToConsole(tcm);
            throw new UserCancelException(msg);
        }
        endTime = System.currentTimeMillis();
        elapsedTime = endTime - startTime;
        msg = "Adding " + outputSpec.getBNGReactions().length + " reactions to model, " + elapsedTime + " ms";
        System.out.println(msg);
        // clean all the reaction rules
        model.getRbmModelContainer().getReactionRuleList().clear();
        // ---- Observables -------------------------------------------------------------------------------------------------
        mathMappingCallback.setMessage("generating network: adding observables...");
        mathMappingCallback.setProgressFraction(progressFractionQuota / 8.0f * 7.0f);
        startTime = System.currentTimeMillis();
        System.out.println("\nObservables :");
        RbmModelContainer rbmmc = model.getRbmModelContainer();
        for (int i = 0; i < outputSpec.getObservableGroups().length; i++) {
            ObservableGroup o = outputSpec.getObservableGroups()[i];
            if (rbmmc.getParameter(o.getObservableGroupName()) != null) {
                System.out.println("   ...already exists.");
                // if it's already there we don't try to add it again; this should be true for all of them!
                continue;
            }
            ArrayList<Expression> terms = new ArrayList<Expression>();
            for (int j = 0; j < o.getListofSpecies().length; j++) {
                Expression term = Expression.mult(new Expression(o.getSpeciesMultiplicity()[j]), new Expression(speciesMap.get(o.getListofSpecies()[j].getNetworkFileIndex())));
                terms.add(term);
            }
            Expression exp = Expression.add(terms.toArray(new Expression[terms.size()])).flatten();
            exp.bindExpression(rbmmc.getSymbolTable());
            RbmObservable originalObservable = rbmmc.getObservable(o.getObservableGroupName());
            VCUnitDefinition observableUnitDefinition = originalObservable.getUnitDefinition();
            rbmmc.removeObservable(originalObservable);
            Parameter newParameter = rbmmc.addParameter(o.getObservableGroupName(), exp, observableUnitDefinition);
            RbmObservable origObservable = simContext.getModel().getRbmModelContainer().getObservable(o.getObservableGroupName());
            ModelEntityMapping em = new ModelEntityMapping(origObservable, newParameter);
            entityMappings.add(em);
        }
        if (mathMappingCallback.isInterrupted()) {
            msg = "Canceled by user.";
            tcm = new TaskCallbackMessage(TaskCallbackStatus.Error, msg);
            simContext.appendToConsole(tcm);
            throw new UserCancelException(msg);
        }
        endTime = System.currentTimeMillis();
        elapsedTime = endTime - startTime;
        msg = "Adding " + outputSpec.getObservableGroups().length + " observables to model, " + elapsedTime + " ms";
        System.out.println(msg);
    } catch (PropertyVetoException ex) {
        ex.printStackTrace(System.out);
        throw new RuntimeException(ex.getMessage());
    } catch (ExpressionBindingException ex) {
        ex.printStackTrace(System.out);
        throw new RuntimeException(ex.getMessage());
    } catch (ModelException ex) {
        ex.printStackTrace(System.out);
        throw new RuntimeException(ex.getMessage());
    } catch (ExpressionException ex) {
        ex.printStackTrace(System.out);
        throw new RuntimeException(ex.getMessage());
    } catch (ClassNotFoundException ex) {
        throw new RuntimeException(ex.getMessage());
    } catch (IOException ex) {
        throw new RuntimeException(ex.getMessage());
    }
    System.out.println("Done transforming");
    msg = "Generating math...";
    System.out.println(msg);
    mathMappingCallback.setMessage(msg);
    mathMappingCallback.setProgressFraction(progressFractionQuota);
}
Also used : HashMap(java.util.HashMap) LinkedHashMap(java.util.LinkedHashMap) UserCancelException(org.vcell.util.UserCancelException) ArrayList(java.util.ArrayList) Product(cbit.vcell.model.Product) SpeciesContext(cbit.vcell.model.SpeciesContext) FakeSeedSpeciesInitialConditionsParameter(org.vcell.model.rbm.FakeSeedSpeciesInitialConditionsParameter) Reactant(cbit.vcell.model.Reactant) BNGOutputSpec(cbit.vcell.bionetgen.BNGOutputSpec) ExpressionException(cbit.vcell.parser.ExpressionException) LinkedHashMap(java.util.LinkedHashMap) FakeReactionRuleRateParameter(org.vcell.model.rbm.FakeReactionRuleRateParameter) KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter) RbmModelContainer(cbit.vcell.model.Model.RbmModelContainer) Species(cbit.vcell.model.Species) BNGSpecies(cbit.vcell.bionetgen.BNGSpecies) HashSet(java.util.HashSet) BNGParameter(cbit.vcell.bionetgen.BNGParameter) ModelException(cbit.vcell.model.ModelException) ObservableGroup(cbit.vcell.bionetgen.ObservableGroup) RbmObservable(cbit.vcell.model.RbmObservable) PropertyVetoException(java.beans.PropertyVetoException) BNGReaction(cbit.vcell.bionetgen.BNGReaction) VCUnitDefinition(cbit.vcell.units.VCUnitDefinition) ReactionStep(cbit.vcell.model.ReactionStep) Map(java.util.Map) HashMap(java.util.HashMap) LinkedHashMap(java.util.LinkedHashMap) SpeciesPattern(org.vcell.model.rbm.SpeciesPattern) Structure(cbit.vcell.model.Structure) SimpleReaction(cbit.vcell.model.SimpleReaction) ReactionRule(cbit.vcell.model.ReactionRule) IOException(java.io.IOException) ExpressionBindingException(cbit.vcell.parser.ExpressionBindingException) Expression(cbit.vcell.parser.Expression) Model(cbit.vcell.model.Model) FakeSeedSpeciesInitialConditionsParameter(org.vcell.model.rbm.FakeSeedSpeciesInitialConditionsParameter) Parameter(cbit.vcell.model.Parameter) KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter) LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter) BNGParameter(cbit.vcell.bionetgen.BNGParameter) FakeReactionRuleRateParameter(org.vcell.model.rbm.FakeReactionRuleRateParameter) MassActionKinetics(cbit.vcell.model.MassActionKinetics) ParseException(org.vcell.model.bngl.ParseException) BNGSpecies(cbit.vcell.bionetgen.BNGSpecies)

Aggregations

Structure (cbit.vcell.model.Structure)159 SpeciesContext (cbit.vcell.model.SpeciesContext)57 Membrane (cbit.vcell.model.Membrane)47 PropertyVetoException (java.beans.PropertyVetoException)42 Feature (cbit.vcell.model.Feature)36 Model (cbit.vcell.model.Model)35 ArrayList (java.util.ArrayList)35 ReactionStep (cbit.vcell.model.ReactionStep)33 Expression (cbit.vcell.parser.Expression)33 ReactionRule (cbit.vcell.model.ReactionRule)27 ExpressionException (cbit.vcell.parser.ExpressionException)27 BioModel (cbit.vcell.biomodel.BioModel)23 StructureMapping (cbit.vcell.mapping.StructureMapping)22 SpeciesPattern (org.vcell.model.rbm.SpeciesPattern)22 Species (cbit.vcell.model.Species)21 MolecularType (org.vcell.model.rbm.MolecularType)20 ReactionParticipant (cbit.vcell.model.ReactionParticipant)19 SimpleReaction (cbit.vcell.model.SimpleReaction)19 SimulationContext (cbit.vcell.mapping.SimulationContext)18 ModelException (cbit.vcell.model.ModelException)18