Search in sources :

Example 6 with FastInvariant

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

the class MathMapping_4_8 method refreshMathDescription.

/**
 * This method was created in VisualAge.
 */
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();
    Model model = simContext.getModel();
    StructureTopology structTopology = model.getStructureTopology();
    // 
    // verify that all structures are mapped to subvolumes and all subvolumes are mapped to a structure
    // 
    Structure[] structures = simContext.getGeometryContext().getModel().getStructures();
    for (int i = 0; i < structures.length; i++) {
        StructureMapping sm = simContext.getGeometryContext().getStructureMapping(structures[i]);
        if (sm == null || (sm instanceof FeatureMapping && getSubVolume((FeatureMapping) sm) == null)) {
            throw new MappingException("model structure '" + structures[i].getName() + "' not mapped to a geometry subdomain");
        }
        if (sm != null && (sm instanceof MembraneMapping) && ((MembraneMapping) sm).getVolumeFractionParameter() != null) {
            Expression volFractExp = ((MembraneMapping) sm).getVolumeFractionParameter().getExpression();
            if (volFractExp != null) {
                try {
                    double volFract = volFractExp.evaluateConstant();
                    if (volFract >= 1.0) {
                        throw new MappingException("model structure '" + structTopology.getInsideFeature(((MembraneMapping) sm).getMembrane()).getName() + "' has volume fraction >= 1.0");
                    }
                } catch (ExpressionException e) {
                }
            }
        }
    }
    SubVolume[] subVolumes = simContext.getGeometryContext().getGeometry().getGeometrySpec().getSubVolumes();
    for (int i = 0; i < subVolumes.length; i++) {
        if (getStructures(subVolumes[i]) == null || getStructures(subVolumes[i]).length == 0) {
            throw new MappingException("geometry subdomain '" + subVolumes[i].getName() + "' not mapped from a model structure");
        }
    }
    // deals with model parameters
    Hashtable<VolVariable, EventAssignmentInitParameter> eventVolVarHash = new Hashtable<VolVariable, EventAssignmentInitParameter>();
    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) {
                for (EventAssignment ea : be.getEventAssignments()) {
                    if (!eventAssignTargets.contains(ea.getTarget())) {
                        eventAssignTargets.add(ea.getTarget());
                    }
                }
            }
        }
        for (int j = 0; j < modelParameters.length; j++) {
            Expression modelParamExpr = getIdentifierSubstitutions(modelParameters[j].getExpression(), modelParameters[j].getUnitDefinition(), null);
            if (eventAssignTargets.contains(modelParameters[j])) {
                EventAssignmentInitParameter eap = null;
                try {
                    eap = addEventAssignmentInitParameter(modelParameters[j].getName(), modelParameters[j].getExpression(), PARAMETER_ROLE_EVENTASSIGN_INITCONDN, modelParameters[j].getUnitDefinition());
                } catch (PropertyVetoException e) {
                    e.printStackTrace(System.out);
                    throw new MappingException(e.getMessage());
                }
                // varHash.addVariable(newFunctionOrConstant(getMathSymbol(eap, null), modelParamExpr));
                VolVariable volVar = new VolVariable(modelParameters[j].getName(), nullDomain);
                varHash.addVariable(volVar);
                eventVolVarHash.put(volVar, eap);
            } else {
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], null), modelParamExpr));
            }
        }
    } else {
        // 
        for (int pass = 0; pass < 2; pass++) {
            for (int j = 0; j < modelParameters.length; j++) {
                Hashtable<String, Expression> structMappingVariantsHash = new Hashtable<String, Expression>();
                for (int k = 0; k < structureMappings.length; k++) {
                    String paramVariantName = null;
                    Expression paramVariantExpr = null;
                    if (modelParameters[j].getExpression().getSymbols() == null) {
                        paramVariantName = modelParameters[j].getName();
                        paramVariantExpr = getIdentifierSubstitutions(modelParameters[j].getExpression(), modelParameters[j].getUnitDefinition(), null);
                    } else {
                        paramVariantName = modelParameters[j].getName() + "_" + TokenMangler.fixTokenStrict(structureMappings[k].getStructure().getName());
                        // if the expression has symbols that do not belong in that structureMapping, do not create the variant.
                        Expression exp1 = modelParameters[j].getExpression();
                        Expression flattenedModelParamExpr = substituteGlobalParameters(exp1);
                        String[] symbols = flattenedModelParamExpr.getSymbols();
                        boolean bValid = true;
                        Structure sm_struct = structureMappings[k].getStructure();
                        if (symbols != null) {
                            for (int ii = 0; ii < symbols.length; ii++) {
                                SpeciesContext sc = model.getSpeciesContext(symbols[ii]);
                                if (sc != null) {
                                    // symbol[ii] is a speciesContext, check its structure with structureMapping[k].structure. If they are the same or
                                    // if it is the adjacent membrane(s), allow variant expression to be created. Else, continue.
                                    Structure sp_struct = sc.getStructure();
                                    if (sp_struct.compareEqual(sm_struct)) {
                                        bValid = bValid && true;
                                    } else {
                                        // if the 2 structures are not the same, are they adjacent? then 'bValid' is true, else false.
                                        if ((sm_struct instanceof Feature) && (sp_struct instanceof Membrane)) {
                                            Feature sm_feature = (Feature) sm_struct;
                                            Membrane sp_mem = (Membrane) sp_struct;
                                            if (sp_mem.compareEqual(structTopology.getParentStructure(sm_feature)) || (structTopology.getInsideFeature(sp_mem).compareEqual(sm_feature) || structTopology.getOutsideFeature(sp_mem).compareEqual(sm_feature))) {
                                                bValid = bValid && true;
                                            } else {
                                                bValid = bValid && false;
                                                break;
                                            }
                                        } else if ((sm_struct instanceof Membrane) && (sp_struct instanceof Feature)) {
                                            Feature sp_feature = (Feature) sp_struct;
                                            Membrane sm_mem = (Membrane) sm_struct;
                                            if (sm_mem.compareEqual(structTopology.getParentStructure(sp_feature)) || (structTopology.getInsideFeature(sm_mem).compareEqual(sp_feature) || structTopology.getOutsideFeature(sm_mem).compareEqual(sp_feature))) {
                                                bValid = bValid && true;
                                            } else {
                                                bValid = bValid && false;
                                                break;
                                            }
                                        } else {
                                            bValid = bValid && false;
                                            break;
                                        }
                                    }
                                }
                            }
                        }
                        if (bValid) {
                            if (pass == 0) {
                                paramVariantExpr = new Expression("VCELL_TEMPORARY_EXPRESSION_PLACEHOLDER");
                            } else {
                                paramVariantExpr = getIdentifierSubstitutions(modelParameters[j].getExpression(), modelParameters[j].getUnitDefinition(), structureMappings[k]);
                            }
                        }
                    }
                    if (paramVariantExpr != null) {
                        structMappingVariantsHash.put(paramVariantName, paramVariantExpr);
                    }
                }
                globalParamVariantsHash.put(modelParameters[j], structMappingVariantsHash);
            }
        }
        // 
        for (int j = 0; j < modelParameters.length; j++) {
            if (modelParameters[j].getExpression().getSymbols() == null) {
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], null), getIdentifierSubstitutions(modelParameters[j].getExpression(), modelParameters[j].getUnitDefinition(), null)));
            } else {
                Hashtable<String, Expression> smVariantsHash = globalParamVariantsHash.get(modelParameters[j]);
                for (int k = 0; k < structureMappings.length; k++) {
                    String variantName = modelParameters[j].getName() + "_" + TokenMangler.fixTokenStrict(structureMappings[k].getStructure().getName());
                    Expression variantExpr = smVariantsHash.get(variantName);
                    if (variantExpr != null) {
                        varHash.addVariable(newFunctionOrConstant(variantName, variantExpr));
                    }
                }
            }
        }
    }
    // 
    // gather only those reactionSteps that are not "excluded"
    // 
    ReactionSpec[] reactionSpecs = simContext.getReactionContext().getReactionSpecs();
    Vector<ReactionStep> rsList = new Vector<ReactionStep>();
    for (int i = 0; i < reactionSpecs.length; i++) {
        if (reactionSpecs[i].isExcluded() == false) {
            rsList.add(reactionSpecs[i].getReactionStep());
        }
    }
    ReactionStep[] reactionSteps = new ReactionStep[rsList.size()];
    rsList.copyInto(reactionSteps);
    // 
    for (int i = 0; i < reactionSteps.length; i++) {
        Kinetics.UnresolvedParameter[] unresolvedParameters = reactionSteps[i].getKinetics().getUnresolvedParameters();
        if (unresolvedParameters != null && unresolvedParameters.length > 0) {
            StringBuffer buffer = new StringBuffer();
            for (int j = 0; j < unresolvedParameters.length; j++) {
                if (j > 0) {
                    buffer.append(", ");
                }
                buffer.append(unresolvedParameters[j].getName());
            }
            throw new MappingException(reactionSteps[i].getDisplayType() + " '" + reactionSteps[i].getName() + "' contains unresolved identifier(s): " + buffer);
        }
    }
    // 
    // create new MathDescription (based on simContext's previous MathDescription if possible)
    // 
    MathDescription oldMathDesc = simContext.getMathDescription();
    mathDesc = null;
    if (oldMathDesc != null) {
        if (oldMathDesc.getVersion() != null) {
            mathDesc = new MathDescription(oldMathDesc.getVersion());
        } else {
            mathDesc = new MathDescription(oldMathDesc.getName());
        }
    } else {
        mathDesc = new MathDescription(simContext.getName() + "_generated");
    }
    // 
    // 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());
        }
    }
    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;
            }
        }
    }
    // (simContext.getGeometry().getDimension() == 0);
    potentialMapping = new PotentialMapping(simContext, this);
    potentialMapping.computeMath();
    if (bCalculatePotential) {
        // 
        // 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), getIdentifierSubstitutions(specificCapacitanceParm.getExpression(), specificCapacitanceParm.getUnitDefinition(), memMapping)));
                ElectricalDevice.ElectricalDeviceParameter transmembraneCurrentParm = membraneElectricalDevice.getParameterFromRole(ElectricalDevice.ROLE_TransmembraneCurrent);
                ElectricalDevice.ElectricalDeviceParameter totalCurrentParm = membraneElectricalDevice.getParameterFromRole(ElectricalDevice.ROLE_TotalCurrent);
                ElectricalDevice.ElectricalDeviceParameter capacitanceParm = membraneElectricalDevice.getParameterFromRole(ElectricalDevice.ROLE_Capacitance);
                if (totalCurrentParm != null && /* totalCurrentDensityParm.getExpression()!=null && */
                memMapping.getCalculateVoltage()) {
                    Expression totalCurrentDensityExp = (totalCurrentParm.getExpression() != null) ? (totalCurrentParm.getExpression()) : (new Expression(0.0));
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrentParm, membraneElectricalDevice.getMembraneMapping()), getIdentifierSubstitutions(totalCurrentDensityExp, totalCurrentParm.getUnitDefinition(), membraneElectricalDevice.getMembraneMapping())));
                }
                if (transmembraneCurrentParm != null && transmembraneCurrentParm.getExpression() != null && memMapping.getCalculateVoltage()) {
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(transmembraneCurrentParm, membraneElectricalDevice.getMembraneMapping()), getIdentifierSubstitutions(transmembraneCurrentParm.getExpression(), transmembraneCurrentParm.getUnitDefinition(), membraneElectricalDevice.getMembraneMapping())));
                }
                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, membraneElectricalDevice.getMembraneMapping()), getIdentifierSubstitutions(Expression.mult(memMapping.getNullSizeParameterValue(), specificCapacitanceParm.getExpression()), capacitanceParm.getUnitDefinition(), membraneElectricalDevice.getMembraneMapping())));
                    } else {
                        varHash.addVariable(newFunctionOrConstant(getMathSymbol(capacitanceParm, membraneElectricalDevice.getMembraneMapping()), getIdentifierSubstitutions(capacitanceParm.getExpression(), capacitanceParm.getUnitDefinition(), membraneElectricalDevice.getMembraneMapping())));
                    }
                }
                // 
                if (membraneElectricalDevice.getDependentVoltageExpression() == null) {
                    // is Voltage Independent?
                    StructureMapping.StructureMappingParameter initialVoltageParm = memMapping.getInitialVoltageParameter();
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(initialVoltageParm, memMapping), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), memMapping)));
                } else // 
                // membrane forced potential
                // 
                {
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(memMapping.getMembrane().getMembraneVoltage(), memMapping), getIdentifierSubstitutions(membraneElectricalDevice.getDependentVoltageExpression(), memMapping.getMembrane().getMembraneVoltage().getUnitDefinition(), memMapping)));
                }
            } 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();
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrentParm, null), getIdentifierSubstitutions(totalCurrentParm.getExpression(), totalCurrentParm.getUnitDefinition(), null)));
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(currentParm, null), getIdentifierSubstitutions(currentParm.getExpression(), currentParm.getUnitDefinition(), null)));
                // 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(), null)));
                    }
                }
            } else if (devices[j] instanceof VoltageClampElectricalDevice) {
                VoltageClampElectricalDevice voltageClampDevice = (VoltageClampElectricalDevice) devices[j];
                // 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, null), getIdentifierSubstitutions(totalCurrent.getExpression(), totalCurrent.getUnitDefinition(), null)));
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(totalCurrentParm, null), getIdentifierSubstitutions(totalCurrentParm.getExpression(), totalCurrentParm.getUnitDefinition(), null)));
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(voltageParm, null), getIdentifierSubstitutions(voltageParm.getExpression(), voltageParm.getUnitDefinition(), null)));
                // 
                // 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], null), getIdentifierSubstitutions(parameters[k].getExpression(), parameters[k].getUnitDefinition(), null)));
                    }
                }
            }
        }
    } 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), getIdentifierSubstitutions(memMapping.getInitialVoltageParameter().getExpression(), memMapping.getInitialVoltageParameter().getUnitDefinition(), memMapping)));
            }
        }
    }
    // 
    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) {
                    if (membraneMapping.getCalculateVoltage() && bCalculatePotential) {
                        if (getResolved(membraneMapping)) {
                            // 
                            if (mathDesc.getVariable(Membrane.MEMBRANE_VOLTAGE_REGION_NAME) == null) {
                                // varHash.addVariable(new MembraneRegionVariable(MembraneVoltage.MEMBRANE_VOLTAGE_REGION_NAME));
                                varHash.addVariable(new MembraneRegionVariable(getMathSymbol(membraneVoltage, membraneMapping), nullDomain));
                            }
                        } else {
                            // 
                            // spatially unresolved membrane, and must solve for potential ... make VolVariable for this compartment
                            // 
                            varHash.addVariable(new VolVariable(getMathSymbol(membraneVoltage, membraneMapping), nullDomain));
                        }
                        Parameter initialVoltageParm = membraneMapping.getInitialVoltageParameter();
                        Variable initVoltageFunction = newFunctionOrConstant(getMathSymbol(initialVoltageParm, membraneMapping), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), membraneMapping));
                        varHash.addVariable(initVoltageFunction);
                    } else {
                        // 
                        // don't calculate voltage, still may need it though
                        // 
                        Parameter initialVoltageParm = membraneMapping.getInitialVoltageParameter();
                        Variable voltageFunction = newFunctionOrConstant(getMathSymbol(membraneMapping.getMembrane().getMembraneVoltage(), membraneMapping), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), membraneMapping));
                        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();
        StructureMapping sm = simContext.getGeometryContext().getStructureMapping(rs.getStructure());
        if (parameters != null) {
            for (int i = 0; i < parameters.length; i++) {
                if (((parameters[i].getRole() == Kinetics.ROLE_CurrentDensity) || (parameters[i].getRole() == Kinetics.ROLE_LumpedCurrent)) && (parameters[i].getExpression() == null || parameters[i].getExpression().isZero())) {
                    continue;
                }
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[i], sm), getIdentifierSubstitutions(parameters[i].getExpression(), parameters[i].getUnitDefinition(), sm)));
            }
        }
    }
    // 
    // initial constants (either function or constant)
    // 
    SpeciesContextSpec[] speciesContextSpecs = simContext.getReactionContext().getSpeciesContextSpecs();
    for (int i = 0; i < speciesContextSpecs.length; i++) {
        SpeciesContextSpecParameter initParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
        if (initParm != null) {
            Expression initExpr = new Expression(initParm.getExpression());
            StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
            String[] symbols = initExpr.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 = initExpr.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);
                        initExpr.substituteInPlace(new Expression(spC.getName()), scsInitExpr);
                    }
                }
            }
            // now create the appropriate function for the current speciesContextSpec.
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(initParm, sm), getIdentifierSubstitutions(initExpr, initParm.getUnitDefinition(), sm)));
        }
    }
    // 
    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), getIdentifierSubstitutions(diffParm.getExpression(), diffParm.getUnitDefinition(), sm)));
        }
    }
    // 
    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), getIdentifierSubstitutions(bc_xm.getExpression(), bc_xm.getUnitDefinition(), sm)));
        }
        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), getIdentifierSubstitutions(bc_xp.getExpression(), bc_xp.getUnitDefinition(), sm)));
        }
        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), getIdentifierSubstitutions(bc_ym.getExpression(), bc_ym.getUnitDefinition(), sm)));
        }
        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), getIdentifierSubstitutions(bc_yp.getExpression(), bc_yp.getUnitDefinition(), sm)));
        }
        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), getIdentifierSubstitutions(bc_zm.getExpression(), bc_zm.getUnitDefinition(), sm)));
        }
        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), getIdentifierSubstitutions(bc_zp.getExpression(), bc_zp.getUnitDefinition(), sm)));
        }
    }
    // 
    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());
        if (advection_velX != null && (advection_velX.getExpression() != null)) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velX, sm), getIdentifierSubstitutions(advection_velX.getExpression(), advection_velX.getUnitDefinition(), sm)));
        }
        SpeciesContextSpec.SpeciesContextSpecParameter advection_velY = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityY);
        if (advection_velY != null && (advection_velY.getExpression() != null)) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velY, sm), getIdentifierSubstitutions(advection_velY.getExpression(), advection_velY.getUnitDefinition(), sm)));
        }
        SpeciesContextSpec.SpeciesContextSpecParameter advection_velZ = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityZ);
        if (advection_velZ != null && (advection_velZ.getExpression() != null)) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velZ, sm), getIdentifierSubstitutions(advection_velZ.getExpression(), advection_velZ.getUnitDefinition(), sm)));
        }
    }
    // 
    // constant species (either function or constant)
    // 
    enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
        if (scm.getVariable() instanceof Constant) {
            varHash.addVariable(scm.getVariable());
        }
    }
    // 
    // conversion factors
    // 
    varHash.addVariable(new Constant(getMathSymbol(model.getKMOLE(), null), 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)));
    // 
    // geometric functions
    // 
    ModelUnitSystem modelUnitSystem = simContext.getModel().getUnitSystem();
    VCUnitDefinition lengthInverseUnit = modelUnitSystem.getLengthUnit().getInverse();
    for (int i = 0; i < structureMappings.length; i++) {
        StructureMapping sm = structureMappings[i];
        Parameter parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumeFraction);
        if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm)));
        }
        parm = sm.getParameterFromRole(StructureMapping.ROLE_SurfaceToVolumeRatio);
        if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm)));
        }
        if (sm instanceof MembraneMapping && !getResolved(sm)) {
            MembraneMapping mm = (MembraneMapping) sm;
            parm = ((MembraneMapping) sm).getVolumeFractionParameter();
            if (parm.getExpression() == null) {
                throw new MappingException("volume fraction not specified for feature '" + structTopology.getInsideFeature(mm.getMembrane()).getName() + "', please refer to Structure Mapping in Application '" + simContext.getName() + "'");
            }
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), modelUnitSystem.getInstance_DIMENSIONLESS(), sm)));
            parm = mm.getSurfaceToVolumeParameter();
            if (parm.getExpression() == null) {
                throw new MappingException("surface to volume ratio not specified for membrane '" + mm.getMembrane().getName() + "', please refer to Structure Mapping in Application '" + simContext.getName() + "'");
            }
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), lengthInverseUnit, sm)));
        }
        StructureMappingParameter sizeParm = sm.getSizeParameter();
        if (sizeParm != null) {
            if (simContext.getGeometry().getDimension() == 0) {
                if (sizeParm.getExpression() != null) {
                    try {
                        double value = sizeParm.getExpression().evaluateConstant();
                        varHash.addVariable(new Constant(getMathSymbol(sizeParm, sm), new Expression(value)));
                    } catch (ExpressionException e) {
                        // varHash.addVariable(new Function(getMathSymbol(parm,sm),getIdentifierSubstitutions(parm.getExpression(),parm.getUnitDefinition(),sm)));
                        e.printStackTrace(System.out);
                        throw new MappingException("Size of structure:" + sm.getNameScope().getName() + " cannot be evaluated as constant.");
                    }
                }
            } else {
                String compartmentName = null;
                VCUnitDefinition sizeUnit = sm.getSizeParameter().getUnitDefinition();
                String sizeFunctionName = null;
                if (sm instanceof MembraneMapping) {
                    MembraneMapping mm = (MembraneMapping) sm;
                    if (getResolved(mm)) {
                        FeatureMapping fm_inside = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(structTopology.getInsideFeature(mm.getMembrane()));
                        FeatureMapping fm_outside = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(structTopology.getOutsideFeature(mm.getMembrane()));
                        compartmentName = getSubVolume(fm_inside).getName() + "_" + getSubVolume(fm_outside).getName();
                        sizeFunctionName = MathFunctionDefinitions.Function_regionArea_current.getFunctionName();
                    } else {
                        FeatureMapping fm_inside = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(structTopology.getInsideFeature(mm.getMembrane()));
                        FeatureMapping fm_outside = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(structTopology.getOutsideFeature(mm.getMembrane()));
                        if (getSubVolume(fm_inside) == getSubVolume(fm_outside)) {
                            compartmentName = getSubVolume(fm_inside).getName();
                            sizeFunctionName = MathFunctionDefinitions.Function_regionVolume_current.getFunctionName();
                        } else {
                            throw new RuntimeException("unexpected structure mapping for membrane '" + mm.getMembrane().getName() + "'");
                        }
                    }
                } else if (sm instanceof FeatureMapping) {
                    FeatureMapping fm = (FeatureMapping) sm;
                    compartmentName = getSubVolume(fm).getName();
                    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), getIdentifierSubstitutions(Expression.mult(totalVolumeCorrection, sizeFunctionExpression), sizeUnit, sm)));
                parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitArea);
                if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm)));
                }
                parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitVolume);
                if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm)));
                }
                parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitArea);
                if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm)));
                }
                parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitVolume);
                if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm)));
                }
            }
        }
    }
    // 
    for (int i = 0; i < fieldMathMappingParameters.length; i++) {
        varHash.addVariable(newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], null), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), null)));
    }
    // 
    // functions
    // 
    enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
        if (scm.getVariable() == null && scm.getDependencyExpression() != null) {
            StructureMapping sm = simContext.getGeometryContext().getStructureMapping(scm.getSpeciesContext().getStructure());
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(scm.getSpeciesContext(), sm), getIdentifierSubstitutions(scm.getDependencyExpression(), scm.getSpeciesContext().getUnitDefinition(), sm)));
        }
    }
    // 
    // 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");
    }
    // 
    // volume subdomains
    // 
    subVolumes = simContext.getGeometryContext().getGeometry().getGeometrySpec().getSubVolumes();
    VCUnitDefinition timeUnit = modelUnitSystem.getTimeUnit();
    for (int j = 0; j < subVolumes.length; j++) {
        SubVolume subVolume = (SubVolume) subVolumes[j];
        // 
        // get priority of subDomain
        // 
        int priority;
        Feature spatialFeature = getResolvedFeature(subVolume);
        if (spatialFeature == null) {
            if (simContext.getGeometryContext().getGeometry().getDimension() > 0) {
                throw new MappingException("no compartment (in Physiology) is mapped to subdomain '" + subVolume.getName() + "' (in Geometry)");
            } else {
                priority = CompartmentSubDomain.NON_SPATIAL_PRIORITY;
            }
        } else {
            // now does not have to match spatial feature, *BUT* needs to be unique
            priority = j;
        }
        // 
        // create subDomain
        // 
        CompartmentSubDomain subDomain = new CompartmentSubDomain(subVolume.getName(), priority);
        mathDesc.addSubDomain(subDomain);
        // 
        if (spatialFeature != null) {
            FeatureMapping fm = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(spatialFeature);
            subDomain.setBoundaryConditionXm(fm.getBoundaryConditionTypeXm());
            subDomain.setBoundaryConditionXp(fm.getBoundaryConditionTypeXp());
            if (simContext.getGeometry().getDimension() > 1) {
                subDomain.setBoundaryConditionYm(fm.getBoundaryConditionTypeYm());
                subDomain.setBoundaryConditionYp(fm.getBoundaryConditionTypeYp());
            }
            if (simContext.getGeometry().getDimension() > 2) {
                subDomain.setBoundaryConditionZm(fm.getBoundaryConditionTypeZm());
                subDomain.setBoundaryConditionZp(fm.getBoundaryConditionTypeZp());
            }
        }
        // 
        // create equations
        // 
        VolumeStructureAnalyzer structureAnalyzer = getVolumeStructureAnalyzer(subVolume);
        Enumeration<SpeciesContextMapping> enumSCM = getSpeciesContextMappings();
        while (enumSCM.hasMoreElements()) {
            SpeciesContextMapping scm = enumSCM.nextElement();
            // 
            if (scm.getVariable() instanceof VolVariable && scm.getDependencyExpression() == null) {
                SpeciesContext sc = scm.getSpeciesContext();
                StructureMapping sm = simContext.getGeometryContext().getStructureMapping(sc.getStructure());
                SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(sc);
                VolVariable variable = (VolVariable) scm.getVariable();
                Equation equation = null;
                if ((scm.isPDERequired()) && sm instanceof FeatureMapping) {
                    // 
                    if (getSubVolume((FeatureMapping) sm) == subVolume) {
                        // 
                        // species context belongs to this subDomain
                        // 
                        Expression initial = new Expression(getMathSymbol(scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration), sm));
                        Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()));
                        Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), sm));
                        equation = new PdeEquation(variable, initial, rate, diffusion);
                        ((PdeEquation) equation).setBoundaryXm((scs.getBoundaryXmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXmParameter(), sm)));
                        ((PdeEquation) equation).setBoundaryXp((scs.getBoundaryXpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXpParameter(), sm)));
                        ((PdeEquation) equation).setBoundaryYm((scs.getBoundaryYmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYmParameter(), sm)));
                        ((PdeEquation) equation).setBoundaryYp((scs.getBoundaryYpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYpParameter(), sm)));
                        ((PdeEquation) equation).setBoundaryZm((scs.getBoundaryZmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZmParameter(), sm)));
                        ((PdeEquation) equation).setBoundaryZp((scs.getBoundaryZpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZpParameter(), sm)));
                        ((PdeEquation) equation).setVelocityX((scs.getVelocityXParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getVelocityXParameter(), sm)));
                        ((PdeEquation) equation).setVelocityY((scs.getVelocityYParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getVelocityYParameter(), sm)));
                        ((PdeEquation) equation).setVelocityZ((scs.getVelocityZParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getVelocityZParameter(), sm)));
                        subDomain.replaceEquation(equation);
                    } else {
                        Expression initial = new Expression(0.0);
                        Expression rate = new Expression(0.0);
                        Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), sm));
                        equation = new PdeEquation(variable, initial, rate, diffusion);
                        if (subDomain.getEquation(variable) == null) {
                            subDomain.addEquation(equation);
                        }
                    }
                } else {
                    // 
                    // ODE
                    // 
                    SubVolume mappedSubVolume = null;
                    if (sm instanceof FeatureMapping) {
                        mappedSubVolume = getSubVolume((FeatureMapping) sm);
                    } else if (sm instanceof MembraneMapping) {
                        // membrane is mapped to that of the inside feature
                        FeatureMapping featureMapping = (FeatureMapping) simContext.getGeometryContext().getStructureMapping(structTopology.getInsideFeature((Membrane) sm.getStructure()));
                        mappedSubVolume = getSubVolume(featureMapping);
                    }
                    if (mappedSubVolume == subVolume) {
                        // 
                        // species context belongs to this subDomain
                        // 
                        Expression initial = new Expression(getMathSymbol(scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration), null));
                        Expression rate = (scm.getRate() == null) ? new Expression(0.0) : getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()));
                        equation = new OdeEquation(variable, initial, rate);
                        subDomain.replaceEquation(equation);
                    } else {
                        Expression initial = new Expression(0.0);
                        Expression rate = new Expression(0.0);
                        equation = new OdeEquation(variable, initial, rate);
                        if (subDomain.getEquation(variable) == null) {
                            subDomain.addEquation(equation);
                        }
                    }
                }
            }
        }
        // 
        // create fast system (if neccessary)
        // 
        SpeciesContextMapping[] fastSpeciesContextMappings = structureAnalyzer.getFastSpeciesContextMappings();
        VCUnitDefinition subDomainUnit = modelUnitSystem.getVolumeConcentrationUnit();
        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), simContext.getGeometryContext().getStructureMapping(getResolvedFeature(subVolume)));
                    FastRate fastRate = new FastRate(rate);
                    fastSystem.addFastRate(fastRate);
                } else {
                    // 
                    // dependant-fast variable, create a fastInvariant object
                    // 
                    Expression rate = getIdentifierSubstitutions(scm.getFastInvariant(), subDomainUnit, simContext.getGeometryContext().getStructureMapping(getResolvedFeature(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 = getStructures(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 (!getResolved(membraneMapping) && membraneMapping.getCalculateVoltage()) {
                    MembraneElectricalDevice capacitiveDevice = potentialMapping.getCapacitiveDevice(membrane);
                    if (capacitiveDevice.getDependentVoltageExpression() == null) {
                        VolVariable vVar = (VolVariable) mathDesc.getVariable(getMathSymbol(capacitiveDevice.getVoltageSymbol(), membraneMapping));
                        Expression initExp = new Expression(getMathSymbol(capacitiveDevice.getMembraneMapping().getInitialVoltageParameter(), membraneMapping));
                        subDomain.addEquation(new OdeEquation(vVar, initExp, getIdentifierSubstitutions(potentialMapping.getOdeRHS(capacitiveDevice, this), membrane.getMembraneVoltage().getUnitDefinition().divideBy(timeUnit), membraneMapping)));
                    } else {
                    // 
                    // 
                    // 
                    }
                }
            }
        }
    }
    // 
    for (int k = 0; k < subVolumes.length; k++) {
        SubVolume subVolume = (SubVolume) subVolumes[k];
        // 
        // if there is a spatially resolved membrane surrounding this subVolume, then create a membraneSubDomain
        // 
        structures = getStructures(subVolume);
        Membrane membrane = null;
        if (structures != null) {
            for (int j = 0; j < structures.length; j++) {
                if (structures[j] instanceof Membrane && getResolved(simContext.getGeometryContext().getStructureMapping(structures[j]))) {
                    membrane = (Membrane) structures[j];
                }
            }
        }
        if (membrane == null) {
            continue;
        }
        SubVolume outerSubVolume = getSubVolume(((FeatureMapping) simContext.getGeometryContext().getStructureMapping(structTopology.getOutsideFeature(membrane))));
        SubVolume innerSubVolume = getSubVolume(((FeatureMapping) simContext.getGeometryContext().getStructureMapping(structTopology.getInsideFeature(membrane))));
        if (innerSubVolume != subVolume) {
            throw new MappingException("membrane " + membrane.getName() + " improperly mapped to inner subVolume " + innerSubVolume.getName());
        }
        // 
        // get priority of subDomain
        // 
        // Feature spatialFeature = simContext.getGeometryContext().getResolvedFeature(subVolume);
        // int priority = spatialFeature.getPriority();
        // 
        // create subDomain
        // 
        CompartmentSubDomain outerCompartment = mathDesc.getCompartmentSubDomain(outerSubVolume.getName());
        CompartmentSubDomain innerCompartment = mathDesc.getCompartmentSubDomain(innerSubVolume.getName());
        SurfaceClass surfaceClass = simContext.getGeometry().getGeometrySurfaceDescription().getSurfaceClass(innerSubVolume, outerSubVolume);
        MembraneSubDomain memSubDomain = new MembraneSubDomain(innerCompartment, outerCompartment, surfaceClass.getName());
        mathDesc.addSubDomain(memSubDomain);
        // 
        // create equations for membrane-bound molecular species
        // 
        MembraneStructureAnalyzer membraneStructureAnalyzer = getMembraneStructureAnalyzer(membrane);
        Enumeration<SpeciesContextMapping> enumSCM = getSpeciesContextMappings();
        while (enumSCM.hasMoreElements()) {
            SpeciesContextMapping scm = enumSCM.nextElement();
            SpeciesContext sc = scm.getSpeciesContext();
            SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(sc);
            // 
            if ((scm.getVariable() instanceof MemVariable) && scm.getDependencyExpression() == null) {
                // 
                // independant variable, create an equation object
                // 
                Equation equation = null;
                MemVariable variable = (MemVariable) scm.getVariable();
                MembraneMapping mm = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(sc.getStructure());
                if (scm.isPDERequired()) {
                    // 
                    if (mm.getMembrane() == membrane) {
                        // 
                        // species context belongs to this subDomain
                        // 
                        Expression initial = new Expression(getMathSymbol(scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration), mm));
                        Expression rate = getIdentifierSubstitutions(scm.getRate(), scm.getSpeciesContext().getUnitDefinition().divideBy(timeUnit), simContext.getGeometryContext().getStructureMapping(sc.getStructure()));
                        Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), mm));
                        equation = new PdeEquation(variable, initial, rate, diffusion);
                        ((PdeEquation) equation).setBoundaryXm((scs.getBoundaryXmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXmParameter(), mm)));
                        ((PdeEquation) equation).setBoundaryXp((scs.getBoundaryXpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryXpParameter(), mm)));
                        ((PdeEquation) equation).setBoundaryYm((scs.getBoundaryYmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYmParameter(), mm)));
                        ((PdeEquation) equation).setBoundaryYp((scs.getBoundaryYpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryYpParameter(), mm)));
                        ((PdeEquation) equation).setBoundaryZm((scs.getBoundaryZmParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZmParameter(), mm)));
                        ((PdeEquation) equation).setBoundaryZp((scs.getBoundaryZpParameter().getExpression() == null) ? (null) : new Expression(getMathSymbol(scs.getBoundaryZpParameter(), mm)));
                        memSubDomain.replaceEquation(equation);
                    } else {
                        Expression initial = new Expression(0.0);
                        Expression rate = new Expression(0.0);
                        Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), mm));
                        equation = new PdeEquation(variable, initial, rate, diffusion);
                        if (memSubDomain.getEquation(variable) == null) {
                            memSubDomain.addEquation(equation);
                        }
                    }
                } else {
                    // 
                    if (mm.getMembrane() == membrane) {
                        // 
                        // 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()));
                        equation = new OdeEquation(variable, initial, rate);
                        memSubDomain.replaceEquation(equation);
                    } else {
                        Expression initial = new Expression(0.0);
                        Expression rate = new Expression(0.0);
                        equation = new OdeEquation(variable, initial, rate);
                        if (memSubDomain.getEquation(variable) == null) {
                            memSubDomain.addEquation(equation);
                        }
                    }
                }
            }
        }
        // 
        // create dummy jump conditions for all volume variables that diffuse and/or advect
        // 
        Enumeration<SpeciesContextMapping> enum_scm = getSpeciesContextMappings();
        while (enum_scm.hasMoreElements()) {
            SpeciesContextMapping scm = enum_scm.nextElement();
            if (scm.isPDERequired()) {
                // Species species = scm.getSpeciesContext().getSpecies();
                Variable var = scm.getVariable();
                if (var instanceof VolVariable && (scm.isPDERequired())) {
                    JumpCondition jc = memSubDomain.getJumpCondition((VolVariable) var);
                    if (jc == null) {
                        // System.out.println("MathMapping.refreshMathDescription(), adding jump condition for diffusing variable "+var.getName()+" on membrane "+membraneStructureAnalyzer.getMembrane().getName());
                        jc = new JumpCondition((VolVariable) var);
                        memSubDomain.addJumpCondition(jc);
                    }
                }
            }
        }
        // 
        // create jump conditions for any volume variables that bind to membrane or have explicitly defined fluxes
        // 
        ResolvedFlux[] resolvedFluxes = membraneStructureAnalyzer.getResolvedFluxes();
        if (resolvedFluxes != null) {
            for (int i = 0; i < resolvedFluxes.length; i++) {
                Species species = resolvedFluxes[i].getSpecies();
                SpeciesContext sc = simContext.getReactionContext().getModel().getSpeciesContext(species, structTopology.getInsideFeature(membraneStructureAnalyzer.getMembrane()));
                if (sc == null) {
                    sc = simContext.getReactionContext().getModel().getSpeciesContext(species, structTopology.getOutsideFeature(membraneStructureAnalyzer.getMembrane()));
                }
                SpeciesContextMapping scm = getSpeciesContextMapping(sc);
                // if (scm.getVariable() instanceof VolVariable && scm.isDiffusing()){
                if (scm.getVariable() instanceof VolVariable && ((MembraneStructureAnalyzer.bNoFluxIfFixed || (scm.isPDERequired())))) {
                    if (MembraneStructureAnalyzer.bNoFluxIfFixed && !scm.isPDERequired()) {
                        MembraneStructureAnalyzer.bNoFluxIfFixedExercised = true;
                    }
                    JumpCondition jc = memSubDomain.getJumpCondition((VolVariable) scm.getVariable());
                    if (jc == null) {
                        jc = new JumpCondition((VolVariable) scm.getVariable());
                        memSubDomain.addJumpCondition(jc);
                    }
                    Expression inFlux = getIdentifierSubstitutions(resolvedFluxes[i].inFluxExpression, resolvedFluxes[i].getUnitDefinition(), simContext.getGeometryContext().getStructureMapping(membraneStructureAnalyzer.getMembrane()));
                    jc.setInFlux(inFlux);
                    Expression outFlux = getIdentifierSubstitutions(resolvedFluxes[i].outFluxExpression, resolvedFluxes[i].getUnitDefinition(), simContext.getGeometryContext().getStructureMapping(membraneStructureAnalyzer.getMembrane()));
                    jc.setOutFlux(outFlux);
                } else {
                    throw new MappingException("APPLICATION  " + simContext.getName() + " : " + scm.getSpeciesContext().getName() + " has spatially resolved flux at membrane " + membrane.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);
                    MembraneMapping membraneMapping = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(membraneStructureAnalyzer.getMembrane());
                    FastRate fastRate = new FastRate(getIdentifierSubstitutions(scm.getFastRate(), rateUnit, membraneMapping));
                    fastSystem.addFastRate(fastRate);
                } else {
                    // 
                    // dependant-fast variable, create a fastInvariant object
                    // 
                    VCUnitDefinition invariantUnit = scm.getSpeciesContext().getUnitDefinition();
                    MembraneMapping membraneMapping = (MembraneMapping) simContext.getGeometryContext().getStructureMapping(membraneStructureAnalyzer.getMembrane());
                    FastInvariant fastInvariant = new FastInvariant(getIdentifierSubstitutions(scm.getFastInvariant(), invariantUnit, membraneMapping));
                    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
        // 
        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)) instanceof MembraneRegionVariable) {
                MembraneRegionVariable vVar = (MembraneRegionVariable) mathDesc.getVariable(getMathSymbol(capacitiveDevice.getVoltageSymbol(), membraneMapping));
                Parameter initialVoltageParm = capacitiveDevice.getMembraneMapping().getInitialVoltageParameter();
                Expression initExp = getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), capacitiveDevice.getMembraneMapping());
                MembraneRegionEquation vEquation = new MembraneRegionEquation(vVar, initExp);
                vEquation.setMembraneRateExpression(getIdentifierSubstitutions(potentialMapping.getOdeRHS(capacitiveDevice, this), membrane.getMembraneVoltage().getUnitDefinition().divideBy(timeUnit), capacitiveDevice.getMembraneMapping()));
                memSubDomain.addEquation(vEquation);
            }
        }
    }
    // create equations for event assign targets that are model params/strutureSize, etc.
    Set<VolVariable> hashKeySet = eventVolVarHash.keySet();
    Iterator<VolVariable> volVarsIter = hashKeySet.iterator();
    // working under teh 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();
        EventAssignmentInitParameter eap = eventVolVarHash.get(volVar);
        Expression rateExpr = new Expression(0.0);
        Equation equation = new OdeEquation(volVar, new Expression(getMathSymbol(eap, null)), rateExpr);
        subDomain.addEquation(equation);
    }
    // events - add events to math desc and odes for event assignments that have parameters as target variables
    BioEvent[] bioevents = simContext.getBioEvents();
    if (bioevents != null && bioevents.length > 0) {
        for (BioEvent be : bioevents) {
            // transform the bioEvent trigger/delay to math Event
            Expression mathTriggerExpr = getIdentifierSubstitutions(be.generateTriggerExpression(), modelUnitSystem.getInstance_DIMENSIONLESS(), null);
            Delay mathDelay = null;
            if (be.getParameter(BioEventParameterType.TriggerDelay) != null) {
                boolean bUseValsFromTriggerTime = be.getUseValuesFromTriggerTime();
                Expression mathDelayExpr = getIdentifierSubstitutions(be.getParameter(BioEventParameterType.TriggerDelay).getExpression(), 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>();
            for (EventAssignment ea : eventAssignments) {
                SymbolTableEntry ste = simContext.getEntry(ea.getTarget().getName());
                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 (!mathDesc.isValid()) {
        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 : MembraneMapping(cbit.vcell.mapping.MembraneMapping) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) ArrayList(java.util.ArrayList) SpeciesContext(cbit.vcell.model.SpeciesContext) StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter) Feature(cbit.vcell.model.Feature) MemVariable(cbit.vcell.math.MemVariable) KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter) SubVolume(cbit.vcell.geometry.SubVolume) Vector(java.util.Vector) SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter) StructureTopology(cbit.vcell.model.Model.StructureTopology) ReactionSpec(cbit.vcell.mapping.ReactionSpec) FastInvariant(cbit.vcell.math.FastInvariant) PropertyVetoException(java.beans.PropertyVetoException) 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) StructureMapping(cbit.vcell.mapping.StructureMapping) FeatureMapping(cbit.vcell.mapping.FeatureMapping) Structure(cbit.vcell.model.Structure) ModelUnitSystem(cbit.vcell.model.ModelUnitSystem) Hashtable(java.util.Hashtable) StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter) Expression(cbit.vcell.parser.Expression) Model(cbit.vcell.model.Model) ProxyParameter(cbit.vcell.model.ProxyParameter) StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter) Parameter(cbit.vcell.model.Parameter) KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter) SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter) LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter) ModelParameter(cbit.vcell.model.Model.ModelParameter) SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter) BioEvent(cbit.vcell.mapping.BioEvent) Event(cbit.vcell.math.Event) BioEvent(cbit.vcell.mapping.BioEvent) MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable) MathDescription(cbit.vcell.math.MathDescription) SpeciesContextMapping(cbit.vcell.mapping.SpeciesContextMapping) SpeciesContextSpec(cbit.vcell.mapping.SpeciesContextSpec) ExpressionException(cbit.vcell.parser.ExpressionException) Delay(cbit.vcell.math.Event.Delay) MappingException(cbit.vcell.mapping.MappingException) PropertyVetoException(java.beans.PropertyVetoException) PdeEquation(cbit.vcell.math.PdeEquation) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) Species(cbit.vcell.model.Species) SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter) EventAssignment(cbit.vcell.mapping.BioEvent.EventAssignment) VolVariable(cbit.vcell.math.VolVariable) ModelParameter(cbit.vcell.model.Model.ModelParameter) OdeEquation(cbit.vcell.math.OdeEquation) JumpCondition(cbit.vcell.math.JumpCondition) MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable) VolVariable(cbit.vcell.math.VolVariable) MemVariable(cbit.vcell.math.MemVariable) Variable(cbit.vcell.math.Variable) SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter) Constant(cbit.vcell.math.Constant) SymbolTableEntry(cbit.vcell.parser.SymbolTableEntry) Membrane(cbit.vcell.model.Membrane) EventAssignment(cbit.vcell.mapping.BioEvent.EventAssignment) OdeEquation(cbit.vcell.math.OdeEquation) MembraneRegionEquation(cbit.vcell.math.MembraneRegionEquation) PdeEquation(cbit.vcell.math.PdeEquation) Equation(cbit.vcell.math.Equation) FastRate(cbit.vcell.math.FastRate)

Example 7 with FastInvariant

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

the class FastSystemAnalyzer method checkLinearity.

/**
 * @exception java.lang.Exception The exception description.
 */
private void checkLinearity() throws MathException, ExpressionException {
    Enumeration<Variable> enum1 = fastVarList.elements();
    while (enum1.hasMoreElements()) {
        Variable var = enum1.nextElement();
        // 
        // d invariant
        // for each variable, make sure ------------- = constant;
        // d Var
        // 
        Enumeration<FastInvariant> enum_fi = fastSystem.getFastInvariants();
        while (enum_fi.hasMoreElements()) {
            FastInvariant fi = enum_fi.nextElement();
            Expression exp = fi.getFunction().differentiate(var.getName());
            exp = MathUtilities.substituteFunctions(exp, this).flatten();
            if (!exp.isNumeric()) {
                // If expression is in terms of 'x','y','z' - then its ok - relax the constant requirement.
                String[] symbols = exp.getSymbols();
                for (int i = 0; i < symbols.length; i++) {
                    if (!symbols[i].equals(ReservedVariable.X.getName()) && !symbols[i].equals(ReservedVariable.Y.getName()) && !symbols[i].equals(ReservedVariable.Z.getName()) && !symbols[i].equals(ReservedVariable.TIME.getName())) {
                        throw new MathException("FastInvariant " + fi.getFunction().toString() + " isn't linear, d/d(" + var.getName() + ") = " + exp.toString());
                    }
                }
            }
        }
    }
}
Also used : ReservedVariable(cbit.vcell.math.ReservedVariable) Variable(cbit.vcell.math.Variable) VolVariable(cbit.vcell.math.VolVariable) MemVariable(cbit.vcell.math.MemVariable) Expression(cbit.vcell.parser.Expression) MathException(cbit.vcell.math.MathException) FastInvariant(cbit.vcell.math.FastInvariant)

Example 8 with FastInvariant

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

the class MathTestingUtilities method constructOdesForSensitivity.

// 
// This method is used to solve for sensitivity of variables to a given parameter.
// The mathDescription and the sensitivity parameter are passed as arguments.
// New variables and ODEs are constructed according to the rule listed below and are added to the mathDescription.
// The method returns the modified mathDescription.
// 
public static MathDescription constructOdesForSensitivity(MathDescription mathDesc, Constant sensParam) throws ExpressionException, MathException, MappingException {
    // 
    // For each ODE :
    // 
    // dX/dt = F(X, P)
    // 
    // (where P is the sensitivity parameter)
    // 
    // we create two other ODEs :
    // 
    // dX_1/dt = F(X_1, P_1)    and
    // 
    // dX_2/dt = F(X_2, P_2)
    // 
    // where P_1 = P + epsilon, and
    // P_2 = P - epsilon.
    // 
    // We keep the initial conditions for both the new ODEs to be the same,
    // i.e., X_1_init = X_2_init.
    // 
    // Then, solving for X_1 & X_2, sensitivity of X wrt P can be computed as :
    // 
    // dX = (X_1 - X_2)
    // --	 -----------   .
    // dP	 (P_1 - P_2)
    // 
    // 
    // REMOVE PRINTS AFTER CHECKING !!!
    System.out.println(" \n\n------------  Old Math Description -----------------");
    System.out.println(mathDesc.getVCML_database());
    if (mathDesc.getGeometry().getDimension() > 0) {
        throw new RuntimeException("Suppport for Spatial systems not yet implemented.");
    }
    VariableHash varHash = new VariableHash();
    Enumeration<Variable> enumVar = mathDesc.getVariables();
    while (enumVar.hasMoreElements()) {
        varHash.addVariable(enumVar.nextElement());
    }
    // 
    // Get 2 values of senstivity parameter (P + epsilon) & (P - epsilon)
    // 
    Constant epsilon = new Constant("epsilon", new Expression(sensParam.getConstantValue() * 1e-3));
    Constant sensParam1 = new Constant(sensParam.getName() + "_1", new Expression(sensParam.getConstantValue() + epsilon.getConstantValue()));
    Constant sensParam2 = new Constant(sensParam.getName() + "_2", new Expression(sensParam.getConstantValue() - epsilon.getConstantValue()));
    // 
    // Iterate through each subdomain (only 1 in compartmental case), and each equation in the subdomain
    // 
    Enumeration<SubDomain> subDomainEnum = mathDesc.getSubDomains();
    // 
    // Create a vector of equations to store the 2 equations for each ODE variable in the subdomain.
    // Later, add it to the equations list in the subdomain.
    // 
    Vector<Equation> equnsVector = new Vector<Equation>();
    Vector<Variable> varsVector = new Vector<Variable>();
    Vector<Variable> var1s = new Vector<Variable>();
    Vector<Variable> var2s = new Vector<Variable>();
    while (subDomainEnum.hasMoreElements()) {
        SubDomain subDomain = subDomainEnum.nextElement();
        Enumeration<Equation> equationEnum = subDomain.getEquations();
        Domain domain = new Domain(subDomain);
        while (equationEnum.hasMoreElements()) {
            Equation equation = equationEnum.nextElement();
            if (equation instanceof OdeEquation) {
                OdeEquation odeEquation = (OdeEquation) equation;
                // Similar to substituteWithExactSolutions, to bind and substitute functions in the ODE
                Expression substitutedRateExp = substituteFunctions(odeEquation.getRateExpression(), mathDesc);
                String varName = odeEquation.getVariable().getName();
                VolVariable var = new VolVariable(varName, domain);
                varsVector.addElement(var);
                // 
                // Create the variable var1, and get the initExpr and rateExpr from the original ODE.
                // Substitute the new vars (var1 and param1) in the old initExpr and rateExpr and create a new ODE
                // 
                String varName1 = new String("__" + varName + "_1");
                Expression initExpr1 = odeEquation.getInitialExpression();
                Expression rateExpr1 = new Expression(substitutedRateExp);
                rateExpr1.substituteInPlace(new Expression(varName), new Expression(varName1));
                rateExpr1.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam1.getName()));
                VolVariable var1 = new VolVariable(varName1, domain);
                var1s.addElement(var1);
                OdeEquation odeEqun1 = new OdeEquation(var1, initExpr1, rateExpr1);
                equnsVector.addElement(odeEqun1);
                // 
                // Create the variable var2, and get the initExpr and rateExpr from the original ODE.
                // Substitute the new vars (var2 and param2) in the old initExpr and rateExpr and create a new ODE
                // 
                String varName2 = new String("__" + varName + "_2");
                Expression initExpr2 = odeEquation.getInitialExpression();
                Expression rateExpr2 = new Expression(substitutedRateExp);
                rateExpr2.substituteInPlace(new Expression(varName), new Expression(varName2));
                rateExpr2.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam2.getName()));
                VolVariable var2 = new VolVariable(varName2, domain);
                var2s.addElement(var2);
                OdeEquation odeEqun2 = new OdeEquation(var2, initExpr2, rateExpr2);
                equnsVector.addElement(odeEqun2);
                // 
                // Create a function for the sensitivity function expression (X1-X2)/(P1-P2), and save in varHash
                // 
                Expression diffVar = Expression.add(new Expression(var1.getName()), Expression.negate(new Expression(var2.getName())));
                Expression diffParam = Expression.add(new Expression(sensParam1.getName()), Expression.negate(new Expression(sensParam2.getName())));
                Expression sensitivityExpr = Expression.mult(diffVar, Expression.invert(diffParam));
                Function sens_Func = new Function("__sens" + varName + "_wrt_" + sensParam.getName(), sensitivityExpr, domain);
                varHash.addVariable(epsilon);
                varHash.addVariable(sensParam1);
                varHash.addVariable(sensParam2);
                varHash.addVariable(var1);
                varHash.addVariable(var2);
                varHash.addVariable(sens_Func);
            } else {
                // sensitivity not implemented for PDEs or other equation types.
                throw new RuntimeException("SolverTest.constructedExactMath(): equation type " + equation.getClass().getName() + " not yet implemented");
            }
        }
        // 
        // Need to substitute the new variables in the new ODEs.
        // i.e., if Rate Expr for ODE_1 for variable X_1 contains variable Y, variable Z, etc.
        // Rate Expr is already substituted with X_1, but it also needs substitute Y with Y_1, Z with Z_1, etc.
        // So get the volume variables, from the vectors for vars, var_1s and var_2s
        // Substitute the rate expressions for the newly added ODEs in equnsVector.
        // 
        Variable[] vars = (Variable[]) BeanUtils.getArray(varsVector, Variable.class);
        Variable[] var_1s = (Variable[]) BeanUtils.getArray(var1s, Variable.class);
        Variable[] var_2s = (Variable[]) BeanUtils.getArray(var2s, Variable.class);
        Vector<Equation> newEqunsVector = new Vector<Equation>();
        for (int i = 0; i < equnsVector.size(); i++) {
            Equation equn = equnsVector.elementAt(i);
            Expression initEx = equn.getInitialExpression();
            Expression rateEx = equn.getRateExpression();
            for (int j = 0; j < vars.length; j++) {
                if (equn.getVariable().getName().endsWith("_1")) {
                    rateEx.substituteInPlace(new Expression(vars[j].getName()), new Expression(var_1s[j].getName()));
                } else if (equn.getVariable().getName().endsWith("_2")) {
                    rateEx.substituteInPlace(new Expression(vars[j].getName()), new Expression(var_2s[j].getName()));
                }
            }
            OdeEquation odeEqun = new OdeEquation(equn.getVariable(), initEx, rateEx);
            newEqunsVector.addElement(odeEqun);
        }
        // 
        for (int i = 0; i < newEqunsVector.size(); i++) {
            mathDesc.getSubDomain(subDomain.getName()).addEquation((Equation) newEqunsVector.elementAt(i));
        }
        // 
        // FAST SYSTEM
        // If the subdomain has a fast system, create a new fast system by substituting the high-low variables/parameters
        // in the expressions for the fastInvariants and fastRates and adding them to the fast system.
        // 
        Vector<FastInvariant> invarsVector = new Vector<FastInvariant>();
        Vector<FastRate> ratesVector = new Vector<FastRate>();
        Enumeration<FastInvariant> fastInvarsEnum = null;
        Enumeration<FastRate> fastRatesEnum = null;
        // Get the fast invariants and fast rates in the system.
        FastSystem fastSystem = subDomain.getFastSystem();
        if (fastSystem != null) {
            fastInvarsEnum = fastSystem.getFastInvariants();
            fastRatesEnum = fastSystem.getFastRates();
            // 
            while (fastInvarsEnum.hasMoreElements()) {
                FastInvariant fastInvar = fastInvarsEnum.nextElement();
                Expression fastInvarExpr = fastInvar.getFunction();
                fastInvarExpr = MathUtilities.substituteFunctions(fastInvarExpr, mathDesc);
                Expression fastInvarExpr1 = new Expression(fastInvarExpr);
                Expression fastInvarExpr2 = new Expression(fastInvarExpr);
                for (int i = 0; i < vars.length; i++) {
                    fastInvarExpr1.substituteInPlace(new Expression(vars[i].getName()), new Expression(var_1s[i].getName()));
                    fastInvarExpr2.substituteInPlace(new Expression(vars[i].getName()), new Expression(var_2s[i].getName()));
                }
                fastInvarExpr1.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam1.getName()));
                FastInvariant fastInvar1 = new FastInvariant(fastInvarExpr1);
                invarsVector.addElement(fastInvar1);
                fastInvarExpr2.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam2.getName()));
                FastInvariant fastInvar2 = new FastInvariant(fastInvarExpr2);
                invarsVector.addElement(fastInvar2);
            }
            // Add the newly created fast invariants to the existing list of fast invariants in the fast system.
            for (int i = 0; i < invarsVector.size(); i++) {
                FastInvariant inVar = (FastInvariant) invarsVector.elementAt(i);
                fastSystem.addFastInvariant(inVar);
            }
            // 
            while (fastRatesEnum.hasMoreElements()) {
                FastRate fastRate = fastRatesEnum.nextElement();
                Expression fastRateExpr = fastRate.getFunction();
                fastRateExpr = MathUtilities.substituteFunctions(fastRateExpr, mathDesc);
                Expression fastRateExpr1 = new Expression(fastRateExpr);
                Expression fastRateExpr2 = new Expression(fastRateExpr);
                for (int i = 0; i < vars.length; i++) {
                    fastRateExpr1.substituteInPlace(new Expression(vars[i].getName()), new Expression(var_1s[i].getName()));
                    fastRateExpr2.substituteInPlace(new Expression(vars[i].getName()), new Expression(var_2s[i].getName()));
                }
                fastRateExpr1.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam1.getName()));
                FastRate fastRate1 = new FastRate(fastRateExpr1);
                ratesVector.addElement(fastRate1);
                fastRateExpr2.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam2.getName()));
                FastRate fastRate2 = new FastRate(fastRateExpr2);
                ratesVector.addElement(fastRate2);
            }
            // Add the newly created fast rates to the existing list of fast rates in the fast system.
            for (int i = 0; i < ratesVector.size(); i++) {
                FastRate rate = (FastRate) ratesVector.elementAt(i);
                fastSystem.addFastRate(rate);
            }
        }
    }
    // Reset all variables in mathDesc.
    mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
    // REMOVE PRINTS AFTER CHECKING
    System.out.println(" \n\n------------  New Math Description -----------------");
    System.out.println(mathDesc.getVCML_database());
    return mathDesc;
}
Also used : MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable) InsideVariable(cbit.vcell.math.InsideVariable) SensVariable(cbit.vcell.solver.ode.SensVariable) FilamentVariable(cbit.vcell.math.FilamentVariable) VolVariable(cbit.vcell.math.VolVariable) VolumeRegionVariable(cbit.vcell.math.VolumeRegionVariable) ReservedVariable(cbit.vcell.math.ReservedVariable) MemVariable(cbit.vcell.math.MemVariable) OutsideVariable(cbit.vcell.math.OutsideVariable) FilamentRegionVariable(cbit.vcell.math.FilamentRegionVariable) Variable(cbit.vcell.math.Variable) VariableHash(cbit.vcell.math.VariableHash) Constant(cbit.vcell.math.Constant) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) Function(cbit.vcell.math.Function) Vector(java.util.Vector) VolVariable(cbit.vcell.math.VolVariable) OdeEquation(cbit.vcell.math.OdeEquation) PdeEquation(cbit.vcell.math.PdeEquation) Equation(cbit.vcell.math.Equation) FastRate(cbit.vcell.math.FastRate) FastInvariant(cbit.vcell.math.FastInvariant) OdeEquation(cbit.vcell.math.OdeEquation) Expression(cbit.vcell.parser.Expression) FastSystem(cbit.vcell.math.FastSystem) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) Domain(cbit.vcell.math.Variable.Domain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)

Example 9 with FastInvariant

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

the class IDAFileWriter method writeEquations.

/**
 * Insert the method's description here.
 * Creation date: (3/8/00 10:31:52 PM)
 */
protected String writeEquations(HashMap<Discontinuity, String> discontinuityNameMap) throws MathException, ExpressionException {
    Simulation simulation = simTask.getSimulation();
    StringBuffer sb = new StringBuffer();
    MathDescription mathDescription = simulation.getMathDescription();
    if (mathDescription.hasFastSystems()) {
        // 
        // define vector of original variables
        // 
        SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
        CompartmentSubDomain subDomain = (CompartmentSubDomain) mathDescription.getSubDomains().nextElement();
        FastSystem fastSystem = subDomain.getFastSystem();
        FastSystemAnalyzer fs_Analyzer = new FastSystemAnalyzer(fastSystem, simSymbolTable);
        int numIndependent = fs_Analyzer.getNumIndependentVariables();
        int systemDim = mathDescription.getStateVariableNames().size();
        int numDependent = systemDim - numIndependent;
        // 
        // get all variables from fast system (dependent and independent)
        // 
        HashSet<String> fastSystemVarHash = new HashSet<String>();
        Enumeration<Variable> dependentfastSystemVarEnum = fs_Analyzer.getDependentVariables();
        while (dependentfastSystemVarEnum.hasMoreElements()) {
            fastSystemVarHash.add(dependentfastSystemVarEnum.nextElement().getName());
        }
        Enumeration<Variable> independentfastSystemVarEnum = fs_Analyzer.getIndependentVariables();
        while (independentfastSystemVarEnum.hasMoreElements()) {
            fastSystemVarHash.add(independentfastSystemVarEnum.nextElement().getName());
        }
        // 
        // get all equations including for variables that are not in the fastSystem (ode equations for "slow system")
        // 
        RationalExpMatrix origInitVector = new RationalExpMatrix(systemDim, 1);
        RationalExpMatrix origSlowRateVector = new RationalExpMatrix(systemDim, 1);
        RationalExpMatrix origVarColumnVector = new RationalExpMatrix(systemDim, 1);
        Enumeration<Equation> enumEquations = subDomain.getEquations();
        int varIndex = 0;
        while (enumEquations.hasMoreElements()) {
            Equation equation = enumEquations.nextElement();
            origVarColumnVector.set_elem(varIndex, 0, new RationalExp(equation.getVariable().getName()));
            Expression rateExpr = equation.getRateExpression();
            rateExpr.bindExpression(varsSymbolTable);
            rateExpr = MathUtilities.substituteFunctions(rateExpr, varsSymbolTable);
            origSlowRateVector.set_elem(varIndex, 0, new RationalExp("(" + rateExpr.flatten().infix() + ")"));
            Expression initExpr = new Expression(equation.getInitialExpression());
            initExpr.substituteInPlace(new Expression("t"), new Expression(0.0));
            initExpr = MathUtilities.substituteFunctions(initExpr, varsSymbolTable).flatten();
            origInitVector.set_elem(varIndex, 0, new RationalExp("(" + initExpr.flatten().infix() + ")"));
            varIndex++;
        }
        // 
        // make symbolic matrix for fast invariants (from FastSystem's fastInvariants as well as a new fast invariant for each variable not included in the fast system.
        // 
        RationalExpMatrix fastInvarianceMatrix = new RationalExpMatrix(numDependent, systemDim);
        int row = 0;
        for (int i = 0; i < origVarColumnVector.getNumRows(); i++) {
            // 
            if (!fastSystemVarHash.contains(origVarColumnVector.get(i, 0).infixString())) {
                fastInvarianceMatrix.set_elem(row, i, RationalExp.ONE);
                row++;
            }
        }
        Enumeration<FastInvariant> enumFastInvariants = fastSystem.getFastInvariants();
        while (enumFastInvariants.hasMoreElements()) {
            FastInvariant fastInvariant = enumFastInvariants.nextElement();
            Expression fastInvariantExpression = fastInvariant.getFunction();
            for (int col = 0; col < systemDim; col++) {
                Expression coeff = fastInvariantExpression.differentiate(origVarColumnVector.get(col, 0).infixString()).flatten();
                coeff = simSymbolTable.substituteFunctions(coeff);
                fastInvarianceMatrix.set_elem(row, col, RationalExpUtils.getRationalExp(coeff));
            }
            row++;
        }
        for (int i = 0; i < systemDim; i++) {
            sb.append("VAR " + origVarColumnVector.get(i, 0).infixString() + " INIT " + origInitVector.get(i, 0).infixString() + ";\n");
        }
        RationalExpMatrix fullMatrix = null;
        RationalExpMatrix inverseFullMatrix = null;
        RationalExpMatrix newSlowRateVector = null;
        try {
            RationalExpMatrix fastMat = ((RationalExpMatrix) fastInvarianceMatrix.transpose().findNullSpace());
            fullMatrix = new RationalExpMatrix(systemDim, systemDim);
            row = 0;
            for (int i = 0; i < fastInvarianceMatrix.getNumRows(); i++) {
                for (int col = 0; col < systemDim; col++) {
                    fullMatrix.set_elem(row, col, fastInvarianceMatrix.get(i, col));
                }
                row++;
            }
            for (int i = 0; i < fastMat.getNumRows(); i++) {
                for (int col = 0; col < systemDim; col++) {
                    fullMatrix.set_elem(row, col, fastMat.get(i, col));
                }
                row++;
            }
            inverseFullMatrix = new RationalExpMatrix(systemDim, systemDim);
            inverseFullMatrix.identity();
            RationalExpMatrix copyOfFullMatrix = new RationalExpMatrix(fullMatrix);
            copyOfFullMatrix.gaussianElimination(inverseFullMatrix);
            newSlowRateVector = new RationalExpMatrix(numDependent, 1);
            newSlowRateVector.matmul(fastInvarianceMatrix, origSlowRateVector);
        } catch (MatrixException ex) {
            ex.printStackTrace();
            throw new MathException(ex.getMessage());
        }
        sb.append("TRANSFORM\n");
        for (row = 0; row < systemDim; row++) {
            for (int col = 0; col < systemDim; col++) {
                sb.append(fullMatrix.get(row, col).getConstant().doubleValue() + " ");
            }
            sb.append("\n");
        }
        sb.append("INVERSETRANSFORM\n");
        for (row = 0; row < systemDim; row++) {
            for (int col = 0; col < systemDim; col++) {
                sb.append(inverseFullMatrix.get(row, col).getConstant().doubleValue() + " ");
            }
            sb.append("\n");
        }
        int numDifferential = numDependent;
        int numAlgebraic = numIndependent;
        sb.append("RHS DIFFERENTIAL " + numDifferential + " ALGEBRAIC " + numAlgebraic + "\n");
        int equationIndex = 0;
        while (equationIndex < numDependent) {
            // print row of mass matrix followed by slow rate corresponding to fast invariant
            Expression slowRateExp = new Expression(newSlowRateVector.get(equationIndex, 0).infixString()).flatten();
            slowRateExp.bindExpression(simSymbolTable);
            slowRateExp = MathUtilities.substituteFunctions(slowRateExp, varsSymbolTable).flatten();
            Vector<Discontinuity> v = slowRateExp.getDiscontinuities();
            for (Discontinuity od : v) {
                od = getSubsitutedAndFlattened(od, varsSymbolTable);
                String dname = discontinuityNameMap.get(od);
                if (dname == null) {
                    dname = ROOT_VARIABLE_PREFIX + discontinuityNameMap.size();
                    discontinuityNameMap.put(od, dname);
                }
                slowRateExp.substituteInPlace(od.getDiscontinuityExp(), new Expression("(" + dname + "==1)"));
            }
            sb.append(slowRateExp.infix() + ";\n");
            equationIndex++;
        }
        Enumeration<FastRate> enumFastRates = fastSystem.getFastRates();
        while (enumFastRates.hasMoreElements()) {
            // print the fastRate for this row
            Expression fastRateExp = new Expression(enumFastRates.nextElement().getFunction());
            fastRateExp = MathUtilities.substituteFunctions(fastRateExp, varsSymbolTable).flatten();
            Vector<Discontinuity> v = fastRateExp.getDiscontinuities();
            for (Discontinuity od : v) {
                od = getSubsitutedAndFlattened(od, varsSymbolTable);
                String dname = discontinuityNameMap.get(od);
                if (dname == null) {
                    dname = ROOT_VARIABLE_PREFIX + discontinuityNameMap.size();
                    discontinuityNameMap.put(od, dname);
                }
                fastRateExp.substituteInPlace(od.getDiscontinuityExp(), new Expression("(" + dname + "==1)"));
            }
            sb.append(fastRateExp.flatten().infix() + ";\n");
            equationIndex++;
        }
    } else {
        for (int i = 0; i < getStateVariableCount(); i++) {
            StateVariable stateVar = getStateVariable(i);
            Expression initExpr = new Expression(stateVar.getInitialRateExpression());
            initExpr = MathUtilities.substituteFunctions(initExpr, varsSymbolTable);
            initExpr.substituteInPlace(new Expression("t"), new Expression(0.0));
            sb.append("VAR " + stateVar.getVariable().getName() + " INIT " + initExpr.flatten().infix() + ";\n");
        }
        sb.append("TRANSFORM\n");
        for (int row = 0; row < getStateVariableCount(); row++) {
            for (int col = 0; col < getStateVariableCount(); col++) {
                sb.append((row == col ? 1 : 0) + " ");
            }
            sb.append("\n");
        }
        sb.append("INVERSETRANSFORM\n");
        for (int row = 0; row < getStateVariableCount(); row++) {
            for (int col = 0; col < getStateVariableCount(); col++) {
                sb.append((row == col ? 1 : 0) + " ");
            }
            sb.append("\n");
        }
        sb.append("RHS DIFFERENTIAL " + getStateVariableCount() + " ALGEBRAIC 0\n");
        for (int i = 0; i < getStateVariableCount(); i++) {
            StateVariable stateVar = getStateVariable(i);
            Expression rateExpr = new Expression(stateVar.getRateExpression());
            rateExpr = MathUtilities.substituteFunctions(rateExpr, varsSymbolTable).flatten();
            Vector<Discontinuity> v = rateExpr.getDiscontinuities();
            for (Discontinuity od : v) {
                od = getSubsitutedAndFlattened(od, varsSymbolTable);
                String dname = discontinuityNameMap.get(od);
                if (dname == null) {
                    dname = ROOT_VARIABLE_PREFIX + discontinuityNameMap.size();
                    discontinuityNameMap.put(od, dname);
                }
                rateExpr.substituteInPlace(od.getDiscontinuityExp(), new Expression("(" + dname + "==1)"));
            }
            sb.append(rateExpr.infix() + ";\n");
        }
    }
    return sb.toString();
}
Also used : Variable(cbit.vcell.math.Variable) MathDescription(cbit.vcell.math.MathDescription) MatrixException(cbit.vcell.matrix.MatrixException) RationalExpMatrix(cbit.vcell.matrix.RationalExpMatrix) HashSet(java.util.HashSet) Discontinuity(cbit.vcell.parser.Discontinuity) SimulationSymbolTable(cbit.vcell.solver.SimulationSymbolTable) Equation(cbit.vcell.math.Equation) FastRate(cbit.vcell.math.FastRate) RationalExp(cbit.vcell.matrix.RationalExp) FastInvariant(cbit.vcell.math.FastInvariant) FastSystemAnalyzer(cbit.vcell.mapping.FastSystemAnalyzer) Simulation(cbit.vcell.solver.Simulation) FastSystem(cbit.vcell.math.FastSystem) Expression(cbit.vcell.parser.Expression) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) MathException(cbit.vcell.math.MathException)

Example 10 with FastInvariant

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

the class XmlReader method getFastSystem.

/**
 * This method returns a FastSystemImplicit from a XML Element.
 * Creation date: (5/18/2001 2:38:56 PM)
 * @return cbit.vcell.math.FastSystemImplicit
 * @param param org.jdom.Element
 * @exception cbit.vcell.xml.XmlParseException The exception description.
 */
private FastSystem getFastSystem(Element param, MathDescription mathDesc) throws XmlParseException {
    // Create a new FastSystem
    FastSystem fastSystem = new FastSystem(mathDesc);
    // Process the FastInvariants
    Iterator<Element> iterator = param.getChildren(XMLTags.FastInvariantTag, vcNamespace).iterator();
    FastInvariant fastInvariant = null;
    while (iterator.hasNext()) {
        Element tempElement = (Element) iterator.next();
        String temp = tempElement.getText();
        try {
            Expression newExp = unMangleExpression(temp);
            fastInvariant = new FastInvariant(newExp);
            fastSystem.addFastInvariant(fastInvariant);
        } catch (MathException e) {
            e.printStackTrace();
            throw new XmlParseException("A MathException was fired when adding the FastInvariant " + fastInvariant + ", to a FastSystem!" + " : ", e);
        }
    }
    // Process the FastRate
    iterator = param.getChildren(XMLTags.FastRateTag, vcNamespace).iterator();
    FastRate fastRate = null;
    while (iterator.hasNext()) {
        Element tempElement = (Element) iterator.next();
        String temp = tempElement.getText();
        try {
            Expression newExp = unMangleExpression(temp);
            fastRate = new FastRate(newExp);
            fastSystem.addFastRate(fastRate);
        } catch (MathException e) {
            e.printStackTrace();
            throw new XmlParseException("A MathException was fired when adding the FastRate " + fastRate + ", to a FastSystem!", e);
        }
    }
    return fastSystem;
}
Also used : FastSystem(cbit.vcell.math.FastSystem) Expression(cbit.vcell.parser.Expression) MathException(cbit.vcell.math.MathException) Element(org.jdom.Element) FastRate(cbit.vcell.math.FastRate) FastInvariant(cbit.vcell.math.FastInvariant)

Aggregations

FastInvariant (cbit.vcell.math.FastInvariant)17 FastRate (cbit.vcell.math.FastRate)13 Variable (cbit.vcell.math.Variable)13 Expression (cbit.vcell.parser.Expression)13 MemVariable (cbit.vcell.math.MemVariable)11 VolVariable (cbit.vcell.math.VolVariable)11 MathException (cbit.vcell.math.MathException)9 ReservedVariable (cbit.vcell.math.ReservedVariable)9 FastSystem (cbit.vcell.math.FastSystem)7 CompartmentSubDomain (cbit.vcell.math.CompartmentSubDomain)5 Constant (cbit.vcell.math.Constant)5 Equation (cbit.vcell.math.Equation)5 SubDomain (cbit.vcell.math.SubDomain)5 RationalExp (cbit.vcell.matrix.RationalExp)5 Function (cbit.vcell.math.Function)4 JumpCondition (cbit.vcell.math.JumpCondition)4 MembraneSubDomain (cbit.vcell.math.MembraneSubDomain)4 OdeEquation (cbit.vcell.math.OdeEquation)4 PdeEquation (cbit.vcell.math.PdeEquation)4 MathDescription (cbit.vcell.math.MathDescription)3