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

Example 16 with MappingException

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

the class StochtestCompareService method compareOne.

public void compareOne() throws IllegalArgumentException, SQLException, DataAccessException, XmlParseException, PropertyVetoException, ExpressionException, MappingException, GeometryException, ImageException, IOException {
    StochtestCompare stochtestCompare = StochtestDbUtils.acceptNextWaitingStochtestCompare(conFactory);
    String biomodelXML = null;
    if (stochtestCompare != null) {
        try {
            StochtestRun stochtestRun1 = StochtestDbUtils.getStochtestRun(conFactory, stochtestCompare.stochtestRun1ref);
            StochtestRun stochtestRun2 = StochtestDbUtils.getStochtestRun(conFactory, stochtestCompare.stochtestRun2ref);
            if (stochtestRun1.status != StochtestRunStatus.complete) {
                throw new RuntimeException("incomplete run status found: " + stochtestRun1.status.name());
            }
            if (stochtestRun2.status != StochtestRunStatus.complete) {
                throw new RuntimeException("incomplete run status found: " + stochtestRun2.status.name());
            }
            TimeSeriesMultitrialData data1 = StochtestFileUtils.readData(StochtestFileUtils.getStochtestRunDataFile(baseDir, stochtestRun1));
            TimeSeriesMultitrialData data2 = StochtestFileUtils.readData(StochtestFileUtils.getStochtestRunDataFile(baseDir, stochtestRun2));
            SummaryStatistics results = TimeSeriesMultitrialData.statisticsSummary(data1, data2);
            XmlUtil.writeXMLStringToFile(results.results(), new File(baseDir, "stochtestcompare_" + stochtestCompare.key + "_summary.txt").getPath(), false);
            StochtestFileUtils.writeVarDiffData(new File(baseDir, "stochtestcompare_" + stochtestCompare.key + "_vardiff.csv"), data1, data2);
            StochtestFileUtils.writeKolmogorovSmirnovTest(new File(baseDir, "stochtestcompare_" + stochtestCompare.key + "_kolmogorovSmirnov.csv"), data1, data2);
            StochtestFileUtils.writeChiSquareTest(new File(baseDir, "stochtestcompare_" + stochtestCompare.key + "_chiSquared.csv"), data1, data2);
            StochtestCompareStatus status = (results.pass()) ? StochtestCompareStatus.not_verydifferent : StochtestCompareStatus.verydifferent;
            StochtestDbUtils.finalizeAcceptedStochtestCompare(conFactory, stochtestCompare, status, null, results);
        } catch (Exception e) {
            StochtestDbUtils.finalizeAcceptedStochtestCompare(conFactory, stochtestCompare, StochtestCompare.StochtestCompareStatus.failed, e.getMessage(), null);
            // 
            // write exception trace to .txt file
            // 
            StringWriter stringWriter = new StringWriter();
            PrintWriter printWriter = new PrintWriter(stringWriter);
            e.printStackTrace(printWriter);
            printWriter.flush();
            System.out.println(stringWriter.getBuffer().toString());
            XmlUtil.writeXMLStringToFile(stringWriter.getBuffer().toString(), new File(baseDir, "stochtestrun_" + stochtestCompare.key + "_error.txt").getPath(), false);
        }
    } else {
        System.out.println("no compare jobs waiting");
        try {
            Thread.sleep(5000);
        } catch (InterruptedException e) {
            // TODO Auto-generated catch block
            e.printStackTrace();
        }
    }
}

Example 17 with MappingException

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

the class StructureMappingCartoonTool method menuAction.

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

Example 18 with MappingException

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

the class StochMathMapping_4_8 method refreshMathDescription.

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

Example 19 with MappingException

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

the class StochMathMapping_4_8 method getProbabilityRate.

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

Example 20 with MappingException

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

the class MathMapping_4_8 method getIdentifierSubstitutions.

/**
 * Substitutes appropriate variables for speciesContext bindings
 *
 * @return cbit.vcell.parser.Expression
 * @param origExp cbit.vcell.parser.Expression
 * @param structureMapping cbit.vcell.mapping.StructureMapping
 */
protected Expression getIdentifierSubstitutions(Expression origExp, VCUnitDefinition desiredExpUnitDef, StructureMapping structureMapping) throws ExpressionException, MappingException {
    String[] symbols = origExp.getSymbols();
    if (symbols == null) {
        return origExp;
    }
    VCUnitDefinition expUnitDef = null;
    try {
        VCUnitEvaluator unitEvaluator = new VCUnitEvaluator(simContext.getModel().getUnitSystem());
        expUnitDef = unitEvaluator.getUnitDefinition(origExp);
        if (desiredExpUnitDef == null) {
            String expStr = origExp.renameBoundSymbols(getNameScope()).infix();
            System.out.println("...........exp='" + expStr + "', desiredUnits are null");
            localIssueList.add(new Issue(origExp, issueContext, IssueCategory.Units, "expected=[null], observed=[" + expUnitDef.getSymbol() + "]", Issue.SEVERITY_WARNING));
        } else if (expUnitDef == null) {
            String expStr = origExp.renameBoundSymbols(getNameScope()).infix();
            System.out.println("...........exp='" + expStr + "', evaluated Units are null");
            localIssueList.add(new Issue(origExp, issueContext, IssueCategory.Units, "expected=[" + desiredExpUnitDef.getSymbol() + "], observed=[null]", Issue.SEVERITY_WARNING));
        } else if (desiredExpUnitDef.isTBD()) {
            String expStr = origExp.renameBoundSymbols(getNameScope()).infix();
            System.out.println("...........exp='" + expStr + "', desiredUnits are [" + desiredExpUnitDef.getSymbol() + "] and expression units are [" + expUnitDef.getSymbol() + "]");
            localIssueList.add(new Issue(origExp, issueContext, IssueCategory.Units, "expected=[" + desiredExpUnitDef.getSymbol() + "], observed=[" + expUnitDef.getSymbol() + "] for exp = " + expStr, Issue.SEVERITY_WARNING));
        } else if (!desiredExpUnitDef.isEquivalent(expUnitDef) && !expUnitDef.isTBD()) {
            String expStr = origExp.renameBoundSymbols(getNameScope()).infix();
            System.out.println("...........exp='" + expStr + "', desiredUnits are [" + desiredExpUnitDef.getSymbol() + "] and expression units are [" + expUnitDef.getSymbol() + "]");
            localIssueList.add(new Issue(origExp, issueContext, IssueCategory.Units, "expected=[" + desiredExpUnitDef.getSymbol() + "], observed=[" + expUnitDef.getSymbol() + "] for exp = " + expStr, Issue.SEVERITY_WARNING));
        }
    } catch (VCUnitException e) {
        String expStr = origExp.renameBoundSymbols(getNameScope()).infix();
        System.out.println(".........exp='" + expStr + "' exception='" + e.getMessage() + "'");
        localIssueList.add(new Issue(origExp, issueContext, IssueCategory.Units, "expected=[" + ((desiredExpUnitDef != null) ? (desiredExpUnitDef.getSymbol()) : ("null")) + "], exception=" + e.getMessage(), Issue.SEVERITY_WARNING));
    } catch (ExpressionException e) {
        String expStr = origExp.renameBoundSymbols(getNameScope()).infix();
        System.out.println(".........exp='" + expStr + "' exception='" + e.getMessage() + "'");
        localIssueList.add(new Issue(origExp, issueContext, IssueCategory.Units, "expected=[" + ((desiredExpUnitDef != null) ? (desiredExpUnitDef.getSymbol()) : ("null")) + "], exception=" + e.getMessage(), Issue.SEVERITY_WARNING));
    } catch (Exception e) {
        e.printStackTrace(System.out);
        localIssueList.add(new Issue(origExp, issueContext, IssueCategory.Units, "expected=[" + ((desiredExpUnitDef != null) ? (desiredExpUnitDef.getSymbol()) : ("null")) + "], exception=" + e.getMessage(), Issue.SEVERITY_WARNING));
    }
    Expression newExp = new Expression(origExp);
    for (int i = 0; i < symbols.length; i++) {
        SymbolTableEntry ste = origExp.getSymbolBinding(symbols[i]);
        if (ste == null) {
            throw new ExpressionBindingException("symbol '" + symbols[i] + "' not bound");
        // ste = simContext.getGeometryContext().getModel().getSpeciesContext(symbols[i]);
        }
        if (ste != null) {
            String newName = getMathSymbol(ste, structureMapping);
            if (!newName.equals(symbols[i])) {
                newExp.substituteInPlace(new Expression(symbols[i]), new Expression(newName));
            }
        }
    }
    return newExp;
}