Search in sources :

Example 56 with VCUnitDefinition

use of cbit.vcell.units.VCUnitDefinition in project vcell by virtualcell.

the class NetworkTransformer method transform.

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

Example 57 with VCUnitDefinition

use of cbit.vcell.units.VCUnitDefinition in project vcell by virtualcell.

the class ParameterContext method gatherIssues.

public void gatherIssues(IssueContext issueContext, List<Issue> issueList, ParameterRoleEnum userDefinedRole) {
    // 
    for (int i = 0; fieldUnresolvedParameters != null && i < fieldUnresolvedParameters.length; i++) {
        issueList.add(new Issue(fieldUnresolvedParameters[i], issueContext, IssueCategory.UnresolvedParameter, "Unresolved parameter '" + fieldUnresolvedParameters[i].getName(), Issue.SEVERITY_ERROR));
    }
    // 
    for (int i = 0; fieldParameters != null && i < fieldParameters.length; i++) {
        if (fieldParameters[i].getRole() == userDefinedRole) {
            try {
                if (!isReferenced(fieldParameters[i], 0)) {
                    issueList.add(new Issue(fieldParameters[i], issueContext, IssueCategory.KineticsUnreferencedParameter, "Unreferenced Kinetic Parameter '" + fieldParameters[i].getName(), Issue.SEVERITY_WARNING));
                }
            } catch (ExpressionException e) {
                issueList.add(new Issue(fieldParameters[i], issueContext, IssueCategory.KineticsExpressionError, "error resolving expression " + e.getMessage(), Issue.SEVERITY_WARNING));
            }
        }
    }
    // 
    if (fieldParameters != null) {
        for (LocalParameter parameter : this.fieldParameters) {
            if (parameter.getExpression() == null) {
                issueList.add(new Issue(parameter, issueContext, IssueCategory.KineticsExpressionMissing, "expression is missing", Issue.SEVERITY_INFO));
            } else {
                Expression exp = parameter.getExpression();
                String[] symbols = exp.getSymbols();
                String issueMessagePrefix = "parameter '" + parameter.getName() + "' ";
                if (symbols != null) {
                    for (int j = 0; j < symbols.length; j++) {
                        SymbolTableEntry ste = exp.getSymbolBinding(symbols[j]);
                        if (ste instanceof LocalProxyParameter) {
                            ste = ((LocalProxyParameter) ste).getTarget();
                        }
                        if (ste == null) {
                            issueList.add(new Issue(parameter, issueContext, IssueCategory.KineticsExpressionUndefinedSymbol, issueMessagePrefix + "references undefined symbol '" + symbols[j] + "'", Issue.SEVERITY_ERROR));
                        // } else if (ste instanceof SpeciesContext) {
                        // if (!getReactionStep().getModel().contains((SpeciesContext)ste)) {
                        // issueList.add(new Issue(parameter,issueContext,IssueCategory.KineticsExpressionUndefinedSymbol, issueMessagePrefix + "references undefined species '"+symbols[j]+"'",Issue.SEVERITY_ERROR));
                        // }
                        // if (reactionStep.countNumReactionParticipants((SpeciesContext)ste) == 0){
                        // issueList.add(new Issue(parameter,issueContext,IssueCategory.KineticsExpressionNonParticipantSymbol, issueMessagePrefix + "references species context '"+symbols[j]+"', but it is not a reactant/product/catalyst of this reaction",Issue.SEVERITY_WARNING));
                        // }
                        // } else if (ste instanceof ModelParameter) {
                        // if (!getReactionStep().getModel().contains((ModelParameter)ste)) {
                        // issueList.add(new Issue(parameter,issueContext,IssueCategory.KineticsExpressionUndefinedSymbol, issueMessagePrefix + "references undefined global parameter '"+symbols[j]+"'",Issue.SEVERITY_ERROR));
                        // }
                        }
                    }
                }
            }
        }
        // looking for local param which masks a global and issueing a warning
        for (LocalParameter parameter : fieldParameters) {
            String name = parameter.getName();
            SymbolTableEntry ste = nameScope.getExternalEntry(name, this);
            String steName;
            if (ste != null) {
                if (ste instanceof Displayable) {
                    steName = ((Displayable) ste).getDisplayType() + " " + ste.getName();
                } else {
                    steName = ste.getClass().getSimpleName() + " " + ste.getName();
                }
                String msg = steName + " is overriden by a local parameter " + name;
                issueList.add(new Issue(parameter, issueContext, IssueCategory.Identifiers, msg, Issue.SEVERITY_WARNING));
            }
        }
    }
    try {
        // 
        // determine unit consistency for each expression
        // 
        VCUnitSystem unitSystem = unitSystemProvider.getUnitSystem();
        VCUnitEvaluator unitEvaluator = new VCUnitEvaluator(unitSystem);
        for (int i = 0; i < fieldParameters.length; i++) {
            if (fieldParameters[i].getExpression() == null) {
                continue;
            }
            try {
                VCUnitDefinition paramUnitDef = fieldParameters[i].getUnitDefinition();
                VCUnitDefinition expUnitDef = unitEvaluator.getUnitDefinition(fieldParameters[i].getExpression());
                if (paramUnitDef == null) {
                    issueList.add(new Issue(fieldParameters[i], issueContext, IssueCategory.Units, "defined unit is null", Issue.SEVERITY_WARNING));
                } else if (paramUnitDef.isTBD()) {
                    issueList.add(new Issue(fieldParameters[i], issueContext, IssueCategory.Units, "undefined unit " + unitSystem.getInstance_TBD().getSymbol(), Issue.SEVERITY_WARNING));
                } else if (expUnitDef == null) {
                    issueList.add(new Issue(fieldParameters[i], issueContext, IssueCategory.Units, "computed unit is null", Issue.SEVERITY_WARNING));
                } else if (paramUnitDef.isTBD() || (!paramUnitDef.isEquivalent(expUnitDef) && !expUnitDef.isTBD())) {
                    issueList.add(new Issue(fieldParameters[i], issueContext, IssueCategory.Units, "inconsistent units, defined=[" + fieldParameters[i].getUnitDefinition().getSymbol() + "], computed=[" + expUnitDef.getSymbol() + "]", Issue.SEVERITY_WARNING));
                }
            } catch (VCUnitException e) {
                issueList.add(new Issue(fieldParameters[i], issueContext, IssueCategory.Units, e.getMessage(), Issue.SEVERITY_WARNING));
            } catch (ExpressionException e) {
                issueList.add(new Issue(fieldParameters[i], issueContext, IssueCategory.Units, e.getMessage(), Issue.SEVERITY_WARNING));
            }
        }
    } catch (Throwable e) {
        issueList.add(new Issue(parameterPolicy.getIssueSource(), issueContext, IssueCategory.Units, "unexpected exception: " + e.getMessage(), Issue.SEVERITY_INFO));
    }
    // 
    for (int i = 0; i < fieldParameters.length; i++) {
        RealInterval simpleBounds = parameterPolicy.getConstraintBounds(fieldParameters[i].getRole());
        if (simpleBounds != null) {
            String parmName = fieldParameters[i].getName();
            issueList.add(new SimpleBoundsIssue(fieldParameters[i], issueContext, simpleBounds, "parameter " + parmName + ": must be within " + simpleBounds.toString()));
        }
    }
}
Also used : Displayable(org.vcell.util.Displayable) VCUnitSystem(cbit.vcell.units.VCUnitSystem) SimpleBoundsIssue(cbit.vcell.model.SimpleBoundsIssue) Issue(org.vcell.util.Issue) SimpleBoundsIssue(cbit.vcell.model.SimpleBoundsIssue) RealInterval(net.sourceforge.interval.ia_math.RealInterval) ExpressionException(cbit.vcell.parser.ExpressionException) VCUnitException(cbit.vcell.units.VCUnitException) VCUnitEvaluator(cbit.vcell.parser.VCUnitEvaluator) SymbolTableEntry(cbit.vcell.parser.SymbolTableEntry) VCUnitDefinition(cbit.vcell.units.VCUnitDefinition) Expression(cbit.vcell.parser.Expression)

Example 58 with VCUnitDefinition

use of cbit.vcell.units.VCUnitDefinition in project vcell by virtualcell.

the class ParticleMathMapping method refreshMathDescription.

/**
 * This method was created in VisualAge.
 */
private void refreshMathDescription() throws MappingException, MatrixException, MathException, ExpressionException, ModelException {
    getSimulationContext().checkValidity();
    if (getSimulationContext().getGeometry().getDimension() == 0) {
        throw new MappingException("particle math mapping requires spatial geometry - dimension >= 1");
    }
    StructureMapping[] structureMappings = getSimulationContext().getGeometryContext().getStructureMappings();
    for (int i = 0; i < structureMappings.length; i++) {
        if (structureMappings[i] instanceof MembraneMapping) {
            if (((MembraneMapping) structureMappings[i]).getCalculateVoltage()) {
                throw new MappingException("electric potential not yet supported for particle models");
            }
        }
    }
    // 
    // fail if any events
    // 
    BioEvent[] bioEvents = getSimulationContext().getBioEvents();
    if (bioEvents != null && bioEvents.length > 0) {
        throw new MappingException("events not yet supported for particle-based models");
    }
    // 
    // gather only those reactionSteps that are not "excluded"
    // 
    ReactionSpec[] reactionSpecs = getSimulationContext().getReactionContext().getReactionSpecs();
    Vector<ReactionStep> rsList = new Vector<ReactionStep>();
    for (int i = 0; i < reactionSpecs.length; i++) {
        if (reactionSpecs[i].isExcluded() == false) {
            if (reactionSpecs[i].isFast()) {
                throw new MappingException("fast reactions not supported for particle models");
            }
            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);
        }
    }
    // 
    // temporarily place all variables in a hashtable (before binding) and discarding duplicates (check for equality)
    // 
    VariableHash varHash = new VariableHash();
    // //
    // // verify that all structures are mapped to geometry classes and all geometry classes are mapped to a structure
    // //
    // Structure structures[] = getSimulationContext().getGeometryContext().getModel().getStructures();
    // for (int i = 0; i < structures.length; i++){
    // StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(structures[i]);
    // if (sm==null || (sm.getGeometryClass() == null)){
    // throw new MappingException("model structure '"+structures[i].getName()+"' not mapped to a geometry subdomain");
    // }
    // if (sm.getUnitSizeParameter()!=null){
    // Expression unitSizeExp = sm.getUnitSizeParameter().getExpression();
    // if(unitSizeExp != null)
    // {
    // try {
    // double unitSize = unitSizeExp.evaluateConstant();
    // if (unitSize != 1.0){
    // throw new MappingException("model structure '"+sm.getStructure().getName()+"' unit size = "+unitSize+" != 1.0 ... partial volume or surface mapping not yet supported for particles");
    // }
    // }catch (ExpressionException e){
    // e.printStackTrace(System.out);
    // throw new MappingException("couldn't evaluate unit size for model structure '"+sm.getStructure().getName()+"' : "+e.getMessage());
    // }
    // }
    // }
    // }
    // {
    // GeometryClass[] geometryClass = getSimulationContext().getGeometryContext().getGeometry().getGeometryClasses();
    // for (int i = 0; i < geometryClass.length; i++){
    // Structure[] mappedStructures = getSimulationContext().getGeometryContext().getStructuresFromGeometryClass(geometryClass[i]);
    // if (mappedStructures==null || mappedStructures.length==0){
    // throw new MappingException("geometryClass '"+geometryClass[i].getName()+"' not mapped from a model structure");
    // }
    // }
    // }
    // deals with model parameters
    Model model = getSimulationContext().getModel();
    ModelUnitSystem modelUnitSystem = model.getUnitSystem();
    ModelParameter[] modelParameters = model.getModelParameters();
    // populate in globalParameterVariants hashtable
    for (int j = 0; j < modelParameters.length; j++) {
        Expression modelParamExpr = modelParameters[j].getExpression();
        GeometryClass geometryClass = getDefaultGeometryClass(modelParamExpr);
        modelParamExpr = getIdentifierSubstitutions(modelParamExpr, modelParameters[j].getUnitDefinition(), geometryClass);
        varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], geometryClass), modelParamExpr, geometryClass));
    }
    // 
    // create new MathDescription (based on simContext's previous MathDescription if possible)
    // 
    MathDescription oldMathDesc = getSimulationContext().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(getSimulationContext().getName() + "_generated");
    }
    // 
    // volume particle variables
    // 
    Enumeration<SpeciesContextMapping> enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = enum1.nextElement();
        if (scm.getVariable() instanceof ParticleVariable) {
            if (!(mathDesc.getVariable(scm.getVariable().getName()) instanceof ParticleVariable)) {
                varHash.addVariable(scm.getVariable());
            }
        }
    }
    varHash.addVariable(new Constant(getMathSymbol(model.getPI_CONSTANT(), null), getIdentifierSubstitutions(model.getPI_CONSTANT().getExpression(), model.getPI_CONSTANT().getUnitDefinition(), null)));
    varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT().getExpression(), model.getFARADAY_CONSTANT().getUnitDefinition(), null)));
    varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT_NMOLE(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT_NMOLE().getExpression(), model.getFARADAY_CONSTANT_NMOLE().getUnitDefinition(), null)));
    varHash.addVariable(new Constant(getMathSymbol(model.getGAS_CONSTANT(), null), getIdentifierSubstitutions(model.getGAS_CONSTANT().getExpression(), model.getGAS_CONSTANT().getUnitDefinition(), null)));
    varHash.addVariable(new Constant(getMathSymbol(model.getTEMPERATURE(), null), getIdentifierSubstitutions(new Expression(getSimulationContext().getTemperatureKelvin()), model.getTEMPERATURE().getUnitDefinition(), null)));
    // 
    for (int j = 0; j < structureMappings.length; j++) {
        if (structureMappings[j] instanceof MembraneMapping) {
            MembraneMapping membraneMapping = (MembraneMapping) structureMappings[j];
            GeometryClass geometryClass = membraneMapping.getGeometryClass();
            // 
            // don't calculate voltage, still may need it though
            // 
            Parameter initialVoltageParm = membraneMapping.getInitialVoltageParameter();
            Variable voltageFunction = newFunctionOrConstant(getMathSymbol(membraneMapping.getMembrane().getMembraneVoltage(), geometryClass), getIdentifierSubstitutions(initialVoltageParm.getExpression(), initialVoltageParm.getUnitDefinition(), geometryClass), geometryClass);
            varHash.addVariable(voltageFunction);
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(membraneMapping.getMembrane().getMembraneVoltage(), membraneMapping.getGeometryClass()), getIdentifierSubstitutions(membraneMapping.getInitialVoltageParameter().getExpression(), membraneMapping.getInitialVoltageParameter().getUnitDefinition(), membraneMapping.getGeometryClass()), membraneMapping.getGeometryClass()));
        }
    }
    // 
    for (int j = 0; j < reactionSteps.length; j++) {
        ReactionStep rs = reactionSteps[j];
        if (getSimulationContext().getReactionContext().getReactionSpec(rs).isExcluded()) {
            continue;
        }
        Kinetics.KineticsParameter[] parameters = rs.getKinetics().getKineticsParameters();
        GeometryClass geometryClass = null;
        if (rs.getStructure() != null) {
            geometryClass = getSimulationContext().getGeometryContext().getStructureMapping(rs.getStructure()).getGeometryClass();
        }
        if (parameters != null) {
            for (int i = 0; i < parameters.length; i++) {
                // Reaction rate, currentDensity, LumpedCurrent and null parameters are not going to displayed in the particle math description.
                if (((parameters[i].getRole() == Kinetics.ROLE_CurrentDensity) || (parameters[i].getRole() == Kinetics.ROLE_LumpedCurrent) || (parameters[i].getRole() == Kinetics.ROLE_ReactionRate)) || (parameters[i].getExpression() == null)) {
                    continue;
                }
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[i], geometryClass), getIdentifierSubstitutions(parameters[i].getExpression(), parameters[i].getUnitDefinition(), geometryClass), geometryClass));
            }
        }
    }
    // 
    // initial constants (either function or constant)
    // 
    SpeciesContextSpec[] speciesContextSpecs = getSimulationContext().getReactionContext().getSpeciesContextSpecs();
    for (int i = 0; i < speciesContextSpecs.length; i++) {
        SpeciesContextSpecParameter initParm = null;
        Expression initExpr = null;
        if (getSimulationContext().isUsingConcentration()) {
            initParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
            initExpr = new Expression(initParm.getExpression());
        // if (speciesContextSpecs[i].getSpeciesContext().getStructure() instanceof Feature) {
        // initExpr = Expression.div(initExpr, new Expression(model.getKMOLE, getNameScope())).flatten();
        // }
        } else {
            initParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_InitialCount);
            initExpr = new Expression(initParm.getExpression());
        }
        if (initExpr != null) {
            StructureMapping sm = getSimulationContext().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 = getSimulationContext().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.getGeometryClass()), getIdentifierSubstitutions(initExpr, initParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
        }
    }
    // 
    for (int i = 0; i < speciesContextSpecs.length; i++) {
        SpeciesContextSpec.SpeciesContextSpecParameter diffParm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_DiffusionRate);
        if (diffParm != null) {
            StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(diffParm, sm.getGeometryClass()), getIdentifierSubstitutions(diffParm.getExpression(), diffParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
        }
    }
    // 
    for (int i = 0; i < speciesContextSpecs.length; i++) {
        SpeciesContextSpec.SpeciesContextSpecParameter bc_xm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueXm);
        StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
        if (bc_xm != null && (bc_xm.getExpression() != null)) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_xm, sm.getGeometryClass()), getIdentifierSubstitutions(bc_xm.getExpression(), bc_xm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
        }
        SpeciesContextSpec.SpeciesContextSpecParameter bc_xp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueXp);
        if (bc_xp != null && (bc_xp.getExpression() != null)) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_xp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_xp.getExpression(), bc_xp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
        }
        SpeciesContextSpec.SpeciesContextSpecParameter bc_ym = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueYm);
        if (bc_ym != null && (bc_ym.getExpression() != null)) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_ym, sm.getGeometryClass()), getIdentifierSubstitutions(bc_ym.getExpression(), bc_ym.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
        }
        SpeciesContextSpec.SpeciesContextSpecParameter bc_yp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueYp);
        if (bc_yp != null && (bc_yp.getExpression() != null)) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_yp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_yp.getExpression(), bc_yp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
        }
        SpeciesContextSpec.SpeciesContextSpecParameter bc_zm = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueZm);
        if (bc_zm != null && (bc_zm.getExpression() != null)) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_zm, sm.getGeometryClass()), getIdentifierSubstitutions(bc_zm.getExpression(), bc_zm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
        }
        SpeciesContextSpec.SpeciesContextSpecParameter bc_zp = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_BoundaryValueZp);
        if (bc_zp != null && (bc_zp.getExpression() != null)) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(bc_zp, sm.getGeometryClass()), getIdentifierSubstitutions(bc_zp.getExpression(), bc_zp.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
        }
    }
    // 
    for (int i = 0; i < speciesContextSpecs.length; i++) {
        SpeciesContextSpec.SpeciesContextSpecParameter advection_velX = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityX);
        StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
        GeometryClass geometryClass = sm.getGeometryClass();
        if (advection_velX != null && (advection_velX.getExpression() != null)) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velX, geometryClass), getIdentifierSubstitutions(advection_velX.getExpression(), advection_velX.getUnitDefinition(), geometryClass), geometryClass));
        }
        SpeciesContextSpec.SpeciesContextSpecParameter advection_velY = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityY);
        if (advection_velY != null && (advection_velY.getExpression() != null)) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velY, geometryClass), getIdentifierSubstitutions(advection_velY.getExpression(), advection_velY.getUnitDefinition(), geometryClass), geometryClass));
        }
        SpeciesContextSpec.SpeciesContextSpecParameter advection_velZ = speciesContextSpecs[i].getParameterFromRole(SpeciesContextSpec.ROLE_VelocityZ);
        if (advection_velZ != null && (advection_velZ.getExpression() != null)) {
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(advection_velZ, geometryClass), getIdentifierSubstitutions(advection_velZ.getExpression(), advection_velZ.getUnitDefinition(), geometryClass), geometryClass));
        }
    }
    // 
    // constant species (either function or constant)
    // 
    enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
        if (scm.getVariable() instanceof Constant) {
            varHash.addVariable(scm.getVariable());
        }
    }
    // 
    // conversion factors
    // 
    varHash.addVariable(new Constant(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.getKMILLIVOLTS(), null), getIdentifierSubstitutions(model.getKMILLIVOLTS().getExpression(), model.getKMILLIVOLTS().getUnitDefinition(), null)));
    varHash.addVariable(new Constant(getMathSymbol(model.getK_GHK(), null), getIdentifierSubstitutions(model.getK_GHK().getExpression(), model.getK_GHK().getUnitDefinition(), null)));
    // 
    for (int i = 0; i < structureMappings.length; i++) {
        StructureMapping sm = structureMappings[i];
        if (getSimulationContext().getGeometry().getDimension() == 0) {
            StructureMappingParameter sizeParm = sm.getSizeParameter();
            if (sizeParm != null && sizeParm.getExpression() != null) {
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(sizeParm, sm.getGeometryClass()), getIdentifierSubstitutions(sizeParm.getExpression(), sizeParm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
            } else {
                if (sm instanceof MembraneMapping) {
                    MembraneMapping mm = (MembraneMapping) sm;
                    StructureMappingParameter volFrac = mm.getVolumeFractionParameter();
                    if (volFrac != null && volFrac.getExpression() != null) {
                        varHash.addVariable(newFunctionOrConstant(getMathSymbol(volFrac, sm.getGeometryClass()), getIdentifierSubstitutions(volFrac.getExpression(), volFrac.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
                    }
                    StructureMappingParameter surfToVol = mm.getSurfaceToVolumeParameter();
                    if (surfToVol != null && surfToVol.getExpression() != null) {
                        varHash.addVariable(newFunctionOrConstant(getMathSymbol(surfToVol, sm.getGeometryClass()), getIdentifierSubstitutions(surfToVol.getExpression(), surfToVol.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
                    }
                }
            }
        } else {
            Parameter parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitArea);
            if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
            }
            parm = sm.getParameterFromRole(StructureMapping.ROLE_AreaPerUnitVolume);
            if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
            }
            parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitArea);
            if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SurfaceClass) {
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
            }
            parm = sm.getParameterFromRole(StructureMapping.ROLE_VolumePerUnitVolume);
            if (parm != null && parm.getExpression() != null && sm.getGeometryClass() instanceof SubVolume) {
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(parm, sm.getGeometryClass()), getIdentifierSubstitutions(parm.getExpression(), parm.getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass()));
            }
        }
    }
    // 
    // functions
    // 
    enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
        if (scm.getVariable() == null && scm.getDependencyExpression() != null) {
            StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(scm.getSpeciesContext().getStructure());
            Variable dependentVariable = newFunctionOrConstant(getMathSymbol(scm.getSpeciesContext(), sm.getGeometryClass()), getIdentifierSubstitutions(scm.getDependencyExpression(), scm.getSpeciesContext().getUnitDefinition(), sm.getGeometryClass()), sm.getGeometryClass());
            dependentVariable.setDomain(new Domain(sm.getGeometryClass()));
            varHash.addVariable(dependentVariable);
        }
    }
    // 
    for (int i = 0; i < fieldMathMappingParameters.length; i++) {
        if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
            GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass()));
        }
    }
    // 
    // set Variables to MathDescription all at once with the order resolved by "VariableHash"
    // 
    mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
    // 
    if (getSimulationContext().getGeometryContext().getGeometry() != null) {
        try {
            mathDesc.setGeometry(getSimulationContext().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");
    }
    // 
    // create subdomains (volume and surfaces)
    // 
    GeometryClass[] geometryClasses = getSimulationContext().getGeometryContext().getGeometry().getGeometryClasses();
    for (int k = 0; k < geometryClasses.length; k++) {
        if (geometryClasses[k] instanceof SubVolume) {
            SubVolume subVolume = (SubVolume) geometryClasses[k];
            // 
            // get priority of subDomain
            // 
            // now does not have to match spatial feature, *BUT* needs to be unique
            int priority = k;
            // 
            // create subDomain
            // 
            CompartmentSubDomain subDomain = new CompartmentSubDomain(subVolume.getName(), priority);
            mathDesc.addSubDomain(subDomain);
            // 
            // assign boundary condition types
            // 
            StructureMapping[] mappedSMs = getSimulationContext().getGeometryContext().getStructureMappings(subVolume);
            FeatureMapping mappedFM = null;
            for (int i = 0; i < mappedSMs.length; i++) {
                if (mappedSMs[i] instanceof FeatureMapping) {
                    if (mappedFM != null) {
                        lg.warn("WARNING:::: MathMapping.refreshMathDescription() ... assigning boundary condition types not unique");
                    }
                    mappedFM = (FeatureMapping) mappedSMs[i];
                }
            }
            if (mappedFM != null) {
                subDomain.setBoundaryConditionXm(mappedFM.getBoundaryConditionTypeXm());
                subDomain.setBoundaryConditionXp(mappedFM.getBoundaryConditionTypeXp());
                if (getSimulationContext().getGeometry().getDimension() > 1) {
                    subDomain.setBoundaryConditionYm(mappedFM.getBoundaryConditionTypeYm());
                    subDomain.setBoundaryConditionYp(mappedFM.getBoundaryConditionTypeYp());
                }
                if (getSimulationContext().getGeometry().getDimension() > 2) {
                    subDomain.setBoundaryConditionZm(mappedFM.getBoundaryConditionTypeZm());
                    subDomain.setBoundaryConditionZp(mappedFM.getBoundaryConditionTypeZp());
                }
            }
        } else if (geometryClasses[k] instanceof SurfaceClass) {
            SurfaceClass surfaceClass = (SurfaceClass) geometryClasses[k];
            // determine membrane inside and outside subvolume
            // this preserves backward compatibility so that membrane subdomain
            // inside and outside correspond to structure hierarchy when present
            Pair<SubVolume, SubVolume> ret = DiffEquMathMapping.computeBoundaryConditionSource(model, simContext, surfaceClass);
            SubVolume innerSubVolume = ret.one;
            SubVolume outerSubVolume = ret.two;
            // 
            // create subDomain
            // 
            CompartmentSubDomain outerCompartment = mathDesc.getCompartmentSubDomain(outerSubVolume.getName());
            CompartmentSubDomain innerCompartment = mathDesc.getCompartmentSubDomain(innerSubVolume.getName());
            MembraneSubDomain memSubDomain = new MembraneSubDomain(innerCompartment, outerCompartment, surfaceClass.getName());
            mathDesc.addSubDomain(memSubDomain);
        }
    }
    // 
    // create Particle Contexts for all Particle Variables
    // 
    Enumeration<SpeciesContextMapping> enumSCM = getSpeciesContextMappings();
    Expression unitFactor = getUnitFactor(modelUnitSystem.getStochasticSubstanceUnit().divideBy(modelUnitSystem.getVolumeSubstanceUnit()));
    while (enumSCM.hasMoreElements()) {
        SpeciesContextMapping scm = enumSCM.nextElement();
        SpeciesContext sc = scm.getSpeciesContext();
        StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure());
        SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
        if (scm.getVariable() instanceof ParticleVariable && scm.getDependencyExpression() == null) {
            ParticleVariable particleVariable = (ParticleVariable) scm.getVariable();
            // 
            // initial distribution of particles
            // 
            ArrayList<ParticleInitialCondition> particleInitialConditions = new ArrayList<ParticleInitialCondition>();
            ParticleInitialCondition pic = null;
            if (getSimulationContext().isUsingConcentration()) {
                Expression initialDistribution = scs.getInitialConcentrationParameter().getExpression() == null ? null : new Expression(getMathSymbol(scs.getInitialConcentrationParameter(), sm.getGeometryClass()));
                if (particleVariable instanceof VolumeParticleVariable) {
                    initialDistribution = Expression.mult(initialDistribution, unitFactor);
                }
                pic = new ParticleInitialConditionConcentration(initialDistribution);
            } else {
                Expression initialCount = scs.getInitialCountParameter().getExpression() == null ? null : new Expression(getMathSymbol(scs.getInitialCountParameter(), sm.getGeometryClass()));
                if (initialCount == null) {
                    throw new MappingException("initialCount not defined for speciesContext " + scs.getSpeciesContext().getName());
                }
                Expression locationX = new Expression("u");
                Expression locationY = new Expression("u");
                Expression locationZ = new Expression("u");
                pic = new ParticleInitialConditionCount(initialCount, locationX, locationY, locationZ);
            }
            particleInitialConditions.add(pic);
            // 
            // diffusion
            // 
            Expression diffusion = new Expression(getMathSymbol(scs.getDiffusionParameter(), sm.getGeometryClass()));
            Expression driftXExp = null;
            if (scs.getVelocityXParameter().getExpression() != null) {
                driftXExp = new Expression(getMathSymbol(scs.getVelocityXParameter(), sm.getGeometryClass()));
            } else {
                SpatialQuantity[] velX_quantities = scs.getVelocityQuantities(QuantityComponent.X);
                if (velX_quantities.length > 0) {
                    int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(sm.getGeometryClass()).length;
                    if (velX_quantities.length == 1 && numRegions == 1) {
                        driftXExp = new Expression(getMathSymbol(velX_quantities[0], sm.getGeometryClass()));
                    } else {
                        throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
                    }
                }
            }
            Expression driftYExp = null;
            if (scs.getVelocityYParameter().getExpression() != null) {
                driftYExp = new Expression(getMathSymbol(scs.getVelocityYParameter(), sm.getGeometryClass()));
            } else {
                SpatialQuantity[] velY_quantities = scs.getVelocityQuantities(QuantityComponent.Y);
                if (velY_quantities.length > 0) {
                    int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(sm.getGeometryClass()).length;
                    if (velY_quantities.length == 1 && numRegions == 1) {
                        driftYExp = new Expression(getMathSymbol(velY_quantities[0], sm.getGeometryClass()));
                    } else {
                        throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
                    }
                }
            }
            Expression driftZExp = null;
            if (scs.getVelocityZParameter().getExpression() != null) {
                driftZExp = new Expression(getMathSymbol(scs.getVelocityZParameter(), sm.getGeometryClass()));
            } else {
                SpatialQuantity[] velZ_quantities = scs.getVelocityQuantities(QuantityComponent.Z);
                if (velZ_quantities.length > 0) {
                    int numRegions = simContext.getGeometry().getGeometrySurfaceDescription().getGeometricRegions(sm.getGeometryClass()).length;
                    if (velZ_quantities.length == 1 && numRegions == 1) {
                        driftZExp = new Expression(getMathSymbol(velZ_quantities[0], sm.getGeometryClass()));
                    } else {
                        throw new MappingException("multiple advection velocities enabled set for multiple volume domains ");
                    }
                }
            }
            ParticleProperties particleProperties = new ParticleProperties(particleVariable, diffusion, driftXExp, driftYExp, driftZExp, particleInitialConditions);
            GeometryClass myGC = sm.getGeometryClass();
            if (myGC == null) {
                throw new MappingException("Application '" + getSimulationContext().getName() + "'\nGeometry->StructureMapping->(" + sm.getStructure().getTypeName() + ")'" + sm.getStructure().getName() + "' must be mapped to geometry domain.\n(see 'Problems' tab)");
            }
            SubDomain subDomain = mathDesc.getSubDomain(myGC.getName());
            subDomain.addParticleProperties(particleProperties);
        }
    }
    for (ReactionStep reactionStep : reactionSteps) {
        Kinetics kinetics = reactionStep.getKinetics();
        StructureMapping sm = getSimulationContext().getGeometryContext().getStructureMapping(reactionStep.getStructure());
        GeometryClass reactionStepGeometryClass = sm.getGeometryClass();
        SubDomain subdomain = mathDesc.getSubDomain(reactionStepGeometryClass.getName());
        KineticsParameter reactionRateParameter = null;
        if (kinetics instanceof LumpedKinetics) {
            reactionRateParameter = ((LumpedKinetics) kinetics).getLumpedReactionRateParameter();
        } else {
            reactionRateParameter = ((DistributedKinetics) kinetics).getReactionRateParameter();
        }
        // macroscopic_irreversible/Microscopic_irreversible for bimolecular membrane reactions. They will NOT go through MassAction solver.
        if (kinetics.getKineticsDescription().equals(KineticsDescription.Macroscopic_irreversible) || kinetics.getKineticsDescription().equals(KineticsDescription.Microscopic_irreversible)) {
            Expression radiusExp = getIdentifierSubstitutions(reactionStep.getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_Binding_Radius).getExpression(), modelUnitSystem.getBindingRadiusUnit(), reactionStepGeometryClass);
            if (radiusExp != null) {
                Expression expCopy = new Expression(radiusExp);
                try {
                    MassActionSolver.substituteParameters(expCopy, true).evaluateConstant();
                } catch (ExpressionException e) {
                    throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Problem in binding radius of " + reactionStep.getName() + ":  '" + radiusExp.infix() + "', " + e.getMessage()));
                }
            } else {
                throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Binding radius of " + reactionStep.getName() + " is null."));
            }
            List<ParticleVariable> reactantParticles = new ArrayList<ParticleVariable>();
            List<ParticleVariable> productParticles = new ArrayList<ParticleVariable>();
            List<Action> forwardActions = new ArrayList<Action>();
            for (ReactionParticipant rp : reactionStep.getReactionParticipants()) {
                SpeciesContext sc = rp.getSpeciesContext();
                SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
                GeometryClass scGeometryClass = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass();
                String varName = getMathSymbol(sc, scGeometryClass);
                Variable var = mathDesc.getVariable(varName);
                if (var instanceof ParticleVariable) {
                    ParticleVariable particle = (ParticleVariable) var;
                    if (rp instanceof Reactant) {
                        reactantParticles.add(particle);
                        if (!scs.isConstant() && !scs.isForceContinuous()) {
                            for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
                                if (radiusExp != null) {
                                    forwardActions.add(Action.createDestroyAction(particle));
                                }
                            }
                        }
                    } else if (rp instanceof Product) {
                        productParticles.add(particle);
                        if (!scs.isConstant() && !scs.isForceContinuous()) {
                            for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
                                if (radiusExp != null) {
                                    forwardActions.add(Action.createCreateAction(particle));
                                }
                            }
                        }
                    }
                } else {
                    throw new MappingException("particle variable '" + varName + "' not found");
                }
            }
            JumpProcessRateDefinition bindingRadius = new InteractionRadius(radiusExp);
            // get jump process name
            String jpName = TokenMangler.mangleToSName(reactionStep.getName());
            // only for NFSim/Rules for now.
            ProcessSymmetryFactor processSymmetryFactor = null;
            if (forwardActions.size() > 0) {
                ParticleJumpProcess forwardProcess = new ParticleJumpProcess(jpName, reactantParticles, bindingRadius, forwardActions, processSymmetryFactor);
                subdomain.addParticleJumpProcess(forwardProcess);
            }
        } else // other type of reactions
        {
            /* 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;
            // Using the MassActionFunction to write out the math description
            MassActionSolver.MassActionFunction maFunc = null;
            if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction) || kinetics.getKineticsDescription().equals(KineticsDescription.General) || kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
                Expression rateExp = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate).getExpression();
                Parameter forwardRateParameter = null;
                Parameter reverseRateParameter = null;
                if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
                    forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KForward);
                    reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KReverse);
                } else if (kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
                    forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
                    reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
                }
                maFunc = MassActionSolver.solveMassAction(forwardRateParameter, reverseRateParameter, 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();
                    }
                }
            }
            if (maFunc != null) {
                // if the reaction has forward rate (Mass action,HMMs), or don't have either forward or reverse rate (some other rate laws--like general)
                // we process it as forward reaction
                List<ParticleVariable> reactantParticles = new ArrayList<ParticleVariable>();
                List<ParticleVariable> productParticles = new ArrayList<ParticleVariable>();
                List<Action> forwardActions = new ArrayList<Action>();
                List<Action> reverseActions = new ArrayList<Action>();
                List<ReactionParticipant> reactants = maFunc.getReactants();
                List<ReactionParticipant> products = maFunc.getProducts();
                for (ReactionParticipant rp : reactants) {
                    SpeciesContext sc = rp.getSpeciesContext();
                    SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
                    GeometryClass scGeometryClass = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass();
                    String varName = getMathSymbol(sc, scGeometryClass);
                    Variable var = mathDesc.getVariable(varName);
                    if (var instanceof ParticleVariable) {
                        ParticleVariable particle = (ParticleVariable) var;
                        reactantParticles.add(particle);
                        if (!scs.isConstant() && !scs.isForceContinuous()) {
                            for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
                                if (forwardRate != null) {
                                    forwardActions.add(Action.createDestroyAction(particle));
                                }
                                if (reverseRate != null) {
                                    reverseActions.add(Action.createCreateAction(particle));
                                }
                            }
                        }
                    } else {
                        throw new MappingException("particle variable '" + varName + "' not found");
                    }
                }
                for (ReactionParticipant rp : products) {
                    SpeciesContext sc = rp.getSpeciesContext();
                    SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
                    GeometryClass scGeometryClass = getSimulationContext().getGeometryContext().getStructureMapping(sc.getStructure()).getGeometryClass();
                    String varName = getMathSymbol(sc, scGeometryClass);
                    Variable var = mathDesc.getVariable(varName);
                    if (var instanceof ParticleVariable) {
                        ParticleVariable particle = (ParticleVariable) var;
                        productParticles.add(particle);
                        if (!scs.isConstant() && !scs.isForceContinuous()) {
                            for (int i = 0; i < Math.abs(rp.getStoichiometry()); i++) {
                                if (forwardRate != null) {
                                    forwardActions.add(Action.createCreateAction(particle));
                                }
                                if (reverseRate != null) {
                                    reverseActions.add(Action.createDestroyAction(particle));
                                }
                            }
                        }
                    } else {
                        throw new MappingException("particle variable '" + varName + "' not found");
                    }
                }
                // 
                // There are two unit conversions required:
                // 
                // 1) convert entire reaction rate from vcell reaction units to Smoldyn units (molecules/lengthunit^dim/timeunit)
                // (where dim is 2 for membrane reactions and 3 for volume reactions)
                // 
                // for forward rates:
                // 2) convert each reactant from Smoldyn units (molecules/lengthunit^dim) to VCell units
                // (where dim is 2 for membrane reactants and 3 for volume reactants)
                // 
                // or
                // 
                // for reverse rates:
                // 2) convert each product from Smoldyn units (molecules/lengthunit^dim) to VCell units
                // (where dim is 2 for membrane products and 3 for volume products)
                // 
                RationalNumber reactionLocationDim = new RationalNumber(reactionStep.getStructure().getDimension());
                VCUnitDefinition timeUnit = modelUnitSystem.getTimeUnit();
                VCUnitDefinition smoldynReactionSizeUnit = modelUnitSystem.getLengthUnit().raiseTo(reactionLocationDim);
                VCUnitDefinition smoldynSubstanceUnit = modelUnitSystem.getStochasticSubstanceUnit();
                VCUnitDefinition smoldynReactionRateUnit = smoldynSubstanceUnit.divideBy(smoldynReactionSizeUnit).divideBy(timeUnit);
                VCUnitDefinition vcellReactionRateUnit = reactionRateParameter.getUnitDefinition();
                VCUnitDefinition reactionUnitFactor = smoldynReactionRateUnit.divideBy(vcellReactionRateUnit);
                if (forwardRate != null) {
                    VCUnitDefinition smoldynReactantsUnit = modelUnitSystem.getInstance_DIMENSIONLESS();
                    // start with factor to translate entire reaction rate.
                    VCUnitDefinition forwardUnitFactor = reactionUnitFactor;
                    // 
                    for (ReactionParticipant reactant : maFunc.getReactants()) {
                        VCUnitDefinition vcellReactantUnit = reactant.getSpeciesContext().getUnitDefinition();
                        boolean bForceContinuous = simContext.getReactionContext().getSpeciesContextSpec(reactant.getSpeciesContext()).isForceContinuous();
                        VCUnitDefinition smoldynReactantUnit = null;
                        if (bForceContinuous) {
                            // reactant is continuous (vcell units)
                            smoldynReactantUnit = reactant.getSpeciesContext().getUnitDefinition();
                        } else {
                            // reactant is a particle (smoldyn units)
                            RationalNumber reactantLocationDim = new RationalNumber(reactant.getStructure().getDimension());
                            VCUnitDefinition smoldynReactantSize = modelUnitSystem.getLengthUnit().raiseTo(reactantLocationDim);
                            smoldynReactantUnit = smoldynSubstanceUnit.divideBy(smoldynReactantSize);
                        }
                        // keep track of units of all reactants
                        smoldynReactantsUnit = smoldynReactantsUnit.multiplyBy(smoldynReactantUnit);
                        RationalNumber reactantStoichiometry = new RationalNumber(reactant.getStoichiometry());
                        VCUnitDefinition reactantUnitFactor = (vcellReactantUnit.divideBy(smoldynReactantUnit)).raiseTo(reactantStoichiometry);
                        // accumulate unit factors for all reactants
                        forwardUnitFactor = forwardUnitFactor.multiplyBy(reactantUnitFactor);
                    }
                    forwardRate = Expression.mult(forwardRate, getUnitFactor(forwardUnitFactor));
                    VCUnitDefinition smoldynExpectedForwardRateUnit = smoldynReactionRateUnit.divideBy(smoldynReactantsUnit);
                    // get probability
                    Expression exp = getIdentifierSubstitutions(forwardRate, smoldynExpectedForwardRateUnit, reactionStepGeometryClass).flatten();
                    JumpProcessRateDefinition partRateDef = new MacroscopicRateConstant(exp);
                    // create particle jump process
                    String jpName = TokenMangler.mangleToSName(reactionStep.getName());
                    // only for NFSim/Rules for now.
                    ProcessSymmetryFactor processSymmetryFactor = null;
                    if (forwardActions.size() > 0) {
                        ParticleJumpProcess forwardProcess = new ParticleJumpProcess(jpName, reactantParticles, partRateDef, forwardActions, processSymmetryFactor);
                        subdomain.addParticleJumpProcess(forwardProcess);
                    }
                }
                // end of forward rate not null
                if (reverseRate != null) {
                    VCUnitDefinition smoldynProductsUnit = modelUnitSystem.getInstance_DIMENSIONLESS();
                    // start with factor to translate entire reaction rate.
                    VCUnitDefinition reverseUnitFactor = reactionUnitFactor;
                    // 
                    for (ReactionParticipant product : maFunc.getProducts()) {
                        VCUnitDefinition vcellProductUnit = product.getSpeciesContext().getUnitDefinition();
                        boolean bForceContinuous = simContext.getReactionContext().getSpeciesContextSpec(product.getSpeciesContext()).isForceContinuous();
                        VCUnitDefinition smoldynProductUnit = null;
                        if (bForceContinuous) {
                            smoldynProductUnit = product.getSpeciesContext().getUnitDefinition();
                        } else {
                            RationalNumber productLocationDim = new RationalNumber(product.getStructure().getDimension());
                            VCUnitDefinition smoldynProductSize = modelUnitSystem.getLengthUnit().raiseTo(productLocationDim);
                            smoldynProductUnit = smoldynSubstanceUnit.divideBy(smoldynProductSize);
                        }
                        // keep track of units of all products
                        smoldynProductsUnit = smoldynProductsUnit.multiplyBy(smoldynProductUnit);
                        RationalNumber productStoichiometry = new RationalNumber(product.getStoichiometry());
                        VCUnitDefinition productUnitFactor = (vcellProductUnit.divideBy(smoldynProductUnit)).raiseTo(productStoichiometry);
                        // accumulate unit factors for all products
                        reverseUnitFactor = reverseUnitFactor.multiplyBy(productUnitFactor);
                    }
                    reverseRate = Expression.mult(reverseRate, getUnitFactor(reverseUnitFactor));
                    VCUnitDefinition smoldynExpectedReverseRateUnit = smoldynReactionRateUnit.divideBy(smoldynProductsUnit);
                    // get probability
                    Expression exp = getIdentifierSubstitutions(reverseRate, smoldynExpectedReverseRateUnit, reactionStepGeometryClass).flatten();
                    JumpProcessRateDefinition partProbRate = new MacroscopicRateConstant(exp);
                    // get jump process name
                    String jpName = TokenMangler.mangleToSName(reactionStep.getName() + "_reverse");
                    // only for NFSim/Rules for now.
                    ProcessSymmetryFactor processSymmetryFactor = null;
                    if (reverseActions.size() > 0) {
                        ParticleJumpProcess reverseProcess = new ParticleJumpProcess(jpName, productParticles, partProbRate, reverseActions, processSymmetryFactor);
                        subdomain.addParticleJumpProcess(reverseProcess);
                    }
                }
            // end of reverse rate not null
            }
        // end of maFunc not null
        }
    // end of reaction step for loop
    }
    // 
    for (int i = 0; i < fieldMathMappingParameters.length; i++) {
        if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
            GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
            Variable variable = newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass());
            if (mathDesc.getVariable(variable.getName()) == null) {
                mathDesc.addVariable(variable);
            }
        }
    }
    if (!mathDesc.isValid()) {
        lg.warn(mathDesc.getVCML_database());
        throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
    }
    if (lg.isDebugEnabled()) {
        System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string begin ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
        System.out.println(mathDesc.getVCML());
        System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string end ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
    }
}
Also used : MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) LumpedKinetics(cbit.vcell.model.LumpedKinetics) MathDescription(cbit.vcell.math.MathDescription) ArrayList(java.util.ArrayList) Product(cbit.vcell.model.Product) SpeciesContext(cbit.vcell.model.SpeciesContext) StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter) Reactant(cbit.vcell.model.Reactant) ExpressionException(cbit.vcell.parser.ExpressionException) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter) SubVolume(cbit.vcell.geometry.SubVolume) Vector(java.util.Vector) SpatialQuantity(cbit.vcell.mapping.spatial.SpatialObject.SpatialQuantity) SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter) SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter) JumpProcessRateDefinition(cbit.vcell.math.JumpProcessRateDefinition) InteractionRadius(cbit.vcell.math.InteractionRadius) ParticleJumpProcess(cbit.vcell.math.ParticleJumpProcess) ModelParameter(cbit.vcell.model.Model.ModelParameter) VCUnitDefinition(cbit.vcell.units.VCUnitDefinition) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) ParticleInitialCondition(cbit.vcell.math.ParticleProperties.ParticleInitialCondition) ReactionStep(cbit.vcell.model.ReactionStep) ParticleProperties(cbit.vcell.math.ParticleProperties) Kinetics(cbit.vcell.model.Kinetics) DistributedKinetics(cbit.vcell.model.DistributedKinetics) LumpedKinetics(cbit.vcell.model.LumpedKinetics) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) Domain(cbit.vcell.math.Variable.Domain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) ReactionParticipant(cbit.vcell.model.ReactionParticipant) GeometryClass(cbit.vcell.geometry.GeometryClass) Action(cbit.vcell.math.Action) VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable) MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable) ParticleVariable(cbit.vcell.math.ParticleVariable) Variable(cbit.vcell.math.Variable) SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter) SurfaceClass(cbit.vcell.geometry.SurfaceClass) VariableHash(cbit.vcell.math.VariableHash) VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable) MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable) ParticleVariable(cbit.vcell.math.ParticleVariable) MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant) Constant(cbit.vcell.math.Constant) SymbolTableEntry(cbit.vcell.parser.SymbolTableEntry) MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant) RationalNumber(ucar.units_vcell.RationalNumber) ModelUnitSystem(cbit.vcell.model.ModelUnitSystem) Pair(org.vcell.util.Pair) ParticleInitialConditionConcentration(cbit.vcell.math.ParticleProperties.ParticleInitialConditionConcentration) ProcessSymmetryFactor(cbit.vcell.math.ParticleJumpProcess.ProcessSymmetryFactor) Expression(cbit.vcell.parser.Expression) VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable) MathException(cbit.vcell.math.MathException) Model(cbit.vcell.model.Model) StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter) Parameter(cbit.vcell.model.Parameter) KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter) SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter) ModelParameter(cbit.vcell.model.Model.ModelParameter) SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter) MassActionSolver(cbit.vcell.model.MassActionSolver) ParticleInitialConditionCount(cbit.vcell.math.ParticleProperties.ParticleInitialConditionCount)

Example 59 with VCUnitDefinition

use of cbit.vcell.units.VCUnitDefinition in project vcell by virtualcell.

the class ReactionContext method refreshSpeciesContextSpecBoundaryUnits.

/**
 * Insert the method's description here.
 * Creation date: (9/15/2004 4:10:23 PM)
 */
void refreshSpeciesContextSpecBoundaryUnits(StructureMapping[] structureMappings) {
    for (int i = 0; i < fieldSpeciesContextSpecs.length; i++) {
        SpeciesContextSpec scs = fieldSpeciesContextSpecs[i];
        VCUnitDefinition dirichletUnit = scs.getSpeciesContext().getUnitDefinition();
        VCUnitDefinition neumannUnit = scs.computeFluxUnit();
        for (int j = 0; j < structureMappings.length; j++) {
            if (structureMappings[j].getStructure().equals(scs.getSpeciesContext().getStructure())) {
                try {
                    // null;
                    StructureMapping fm = structureMappings[j];
                    // if (structureMappings[j] instanceof FeatureMapping){
                    // fm = (FeatureMapping)structureMappings[j];
                    // }else if (structureMappings[j] instanceof MembraneMapping){
                    // //
                    // // if a membrane, use boundary condition type from the enclosed volume compartment
                    // // (e.g. Plasma membrane uses cytosol boundary condition types).
                    // //
                    // Feature insideFeature = ((MembraneMapping)structureMappings[j]).getMembrane().getInsideFeature();
                    // for (int k = 0; k < structureMappings.length; k++){
                    // if (structureMappings[k].getStructure() == insideFeature){
                    // fm = (FeatureMapping)structureMappings[k];
                    // }
                    // }
                    // }
                    // 
                    // SpeciesContextSpec is a volume species mapped to this resolved Feature
                    // 
                    // 
                    // for the XM direction
                    // 
                    {
                        SpeciesContextSpec.SpeciesContextSpecParameter xmParm = scs.getBoundaryXmParameter();
                        if (fm.getBoundaryConditionTypeXm().isDIRICHLET() && !xmParm.getUnitDefinition().compareEqual(dirichletUnit)) {
                            xmParm.setUnitDefinition(dirichletUnit);
                        } else if (fm.getBoundaryConditionTypeXm().isNEUMANN() && !xmParm.getUnitDefinition().compareEqual(neumannUnit)) {
                            xmParm.setUnitDefinition(neumannUnit);
                        }
                    }
                    // 
                    // for the XP direction
                    // 
                    {
                        SpeciesContextSpec.SpeciesContextSpecParameter xpParm = scs.getBoundaryXpParameter();
                        if (fm.getBoundaryConditionTypeXp().isDIRICHLET() && !xpParm.getUnitDefinition().compareEqual(dirichletUnit)) {
                            xpParm.setUnitDefinition(dirichletUnit);
                        } else if (fm.getBoundaryConditionTypeXp().isNEUMANN() && !xpParm.getUnitDefinition().compareEqual(neumannUnit)) {
                            xpParm.setUnitDefinition(neumannUnit);
                        }
                    }
                    // 
                    // for the YM direction
                    // 
                    {
                        SpeciesContextSpec.SpeciesContextSpecParameter ymParm = scs.getBoundaryYmParameter();
                        if (fm.getBoundaryConditionTypeYm().isDIRICHLET() && !ymParm.getUnitDefinition().compareEqual(dirichletUnit)) {
                            ymParm.setUnitDefinition(dirichletUnit);
                        } else if (fm.getBoundaryConditionTypeYm().isNEUMANN() && !ymParm.getUnitDefinition().compareEqual(neumannUnit)) {
                            ymParm.setUnitDefinition(neumannUnit);
                        }
                    }
                    // 
                    // for the YP direction
                    // 
                    {
                        SpeciesContextSpec.SpeciesContextSpecParameter ypParm = scs.getBoundaryYpParameter();
                        if (fm.getBoundaryConditionTypeYp().isDIRICHLET() && !ypParm.getUnitDefinition().compareEqual(dirichletUnit)) {
                            ypParm.setUnitDefinition(dirichletUnit);
                        } else if (fm.getBoundaryConditionTypeYp().isNEUMANN() && !ypParm.getUnitDefinition().compareEqual(neumannUnit)) {
                            ypParm.setUnitDefinition(neumannUnit);
                        }
                    }
                    // 
                    // for the ZM direction
                    // 
                    {
                        SpeciesContextSpec.SpeciesContextSpecParameter zmParm = scs.getBoundaryZmParameter();
                        if (fm.getBoundaryConditionTypeZm().isDIRICHLET() && !zmParm.getUnitDefinition().compareEqual(dirichletUnit)) {
                            zmParm.setUnitDefinition(dirichletUnit);
                        } else if (fm.getBoundaryConditionTypeZm().isNEUMANN() && !zmParm.getUnitDefinition().compareEqual(neumannUnit)) {
                            zmParm.setUnitDefinition(neumannUnit);
                        }
                    }
                    // 
                    // for the ZP direction
                    // 
                    {
                        SpeciesContextSpec.SpeciesContextSpecParameter zpParm = scs.getBoundaryZpParameter();
                        if (fm.getBoundaryConditionTypeZp().isDIRICHLET() && !zpParm.getUnitDefinition().compareEqual(dirichletUnit)) {
                            zpParm.setUnitDefinition(dirichletUnit);
                        } else if (fm.getBoundaryConditionTypeZp().isNEUMANN() && !zpParm.getUnitDefinition().compareEqual(neumannUnit)) {
                            zpParm.setUnitDefinition(neumannUnit);
                        }
                    }
                } catch (java.beans.PropertyVetoException e) {
                    e.printStackTrace(System.out);
                }
            }
        }
    }
}
Also used : PropertyVetoException(java.beans.PropertyVetoException) VCUnitDefinition(cbit.vcell.units.VCUnitDefinition) SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter)

Example 60 with VCUnitDefinition

use of cbit.vcell.units.VCUnitDefinition in project vcell by virtualcell.

the class ITextWriter method writeMembraneMapping.

protected void writeMembraneMapping(Section simContextSection, SimulationContext simContext) throws DocumentException {
    GeometryContext geoContext = simContext.getGeometryContext();
    if (geoContext == null) {
        return;
    }
    Section memMapSection = null;
    Table memMapTable = null;
    StructureMapping[] structMappings = geoContext.getStructureMappings();
    for (int i = 0; i < structMappings.length; i++) {
        MembraneMapping memMapping = null;
        if (structMappings[i] instanceof FeatureMapping) {
            continue;
        } else {
            memMapping = (MembraneMapping) structMappings[i];
        }
        String structName = memMapping.getStructure().getName();
        String initVoltage = "";
        Expression tempExp = memMapping.getInitialVoltageParameter().getExpression();
        VCUnitDefinition tempUnit = memMapping.getInitialVoltageParameter().getUnitDefinition();
        if (tempExp != null) {
            initVoltage = tempExp.infix();
            if (tempUnit != null) {
                initVoltage += "   " + tempUnit.getSymbolUnicode();
            }
        }
        String spCap = "";
        tempExp = memMapping.getSpecificCapacitanceParameter().getExpression();
        tempUnit = memMapping.getSpecificCapacitanceParameter().getUnitDefinition();
        if (tempExp != null) {
            spCap = tempExp.infix();
            if (tempUnit != null) {
                spCap += "   " + tempUnit.getSymbolUnicode();
            }
        }
        if (memMapTable == null) {
            memMapTable = getTable(4, 100, 1, 3, 3);
            memMapTable.addCell(createCell("Electrical Mapping - Membrane Potential", getBold(DEF_HEADER_FONT_SIZE), 4, 1, Element.ALIGN_CENTER, true));
            memMapTable.addCell(createHeaderCell("Membrane", getBold(), 1));
            memMapTable.addCell(createHeaderCell("Calculate V (T/F)", getBold(), 1));
            memMapTable.addCell(createHeaderCell("V initial", getBold(), 1));
            memMapTable.addCell(createHeaderCell("Specific Capacitance", getBold(), 1));
            memMapTable.endHeaders();
        }
        memMapTable.addCell(createCell(structName, getFont()));
        memMapTable.addCell(createCell((memMapping.getCalculateVoltage() ? " T " : " F "), getFont()));
        memMapTable.addCell(createCell(initVoltage, getFont()));
        memMapTable.addCell(createCell(spCap, getFont()));
    }
    if (memMapTable != null) {
        memMapSection = simContextSection.addSection("Membrane Mapping For " + simContext.getName(), simContextSection.numberDepth() + 1);
        memMapSection.add(memMapTable);
    }
    int[] widths = { 1, 1, 1, 5, 8 };
    Table electTable = null;
    ElectricalStimulus[] electricalStimuli = simContext.getElectricalStimuli();
    for (int j = 0; j < electricalStimuli.length; j++) {
        if (j == 0) {
            electTable = getTable(5, 100, 1, 3, 3);
            electTable.addCell(createCell("Electrical Mapping - Electrical Stimulus", getBold(DEF_HEADER_FONT_SIZE), 5, 1, Element.ALIGN_CENTER, true));
            electTable.addCell(createHeaderCell("Stimulus Name", getBold(), 1));
            electTable.addCell(createHeaderCell("Current Name", getBold(), 1));
            electTable.addCell(createHeaderCell("Clamp Type", getBold(), 1));
            electTable.addCell(createHeaderCell("Voltage/Current Density", getBold(), 1));
            electTable.addCell(createHeaderCell("Clamp Device", getBold(), 1));
            electTable.endHeaders();
        }
        String stimName = electricalStimuli[j].getName();
        String currName = "";
        String clampType = "", expStr = "";
        Expression tempExp = null;
        VCUnitDefinition tempUnit = null;
        if (electricalStimuli[j] instanceof CurrentDensityClampStimulus) {
            CurrentDensityClampStimulus stimulus = (CurrentDensityClampStimulus) electricalStimuli[j];
            LocalParameter currentDensityParameter = stimulus.getCurrentDensityParameter();
            tempExp = currentDensityParameter.getExpression();
            tempUnit = currentDensityParameter.getUnitDefinition();
            clampType = "Current Density (deprecated)";
        } else if (electricalStimuli[j] instanceof TotalCurrentClampStimulus) {
            TotalCurrentClampStimulus stimulus = (TotalCurrentClampStimulus) electricalStimuli[j];
            LocalParameter totalCurrentParameter = stimulus.getCurrentParameter();
            tempExp = totalCurrentParameter.getExpression();
            tempUnit = totalCurrentParameter.getUnitDefinition();
            clampType = "Current";
        } else if (electricalStimuli[j] instanceof VoltageClampStimulus) {
            VoltageClampStimulus stimulus = (VoltageClampStimulus) electricalStimuli[j];
            Parameter voltageParameter = stimulus.getVoltageParameter();
            tempExp = voltageParameter.getExpression();
            tempUnit = voltageParameter.getUnitDefinition();
            clampType = "Voltage";
        }
        if (tempExp != null) {
            expStr = tempExp.infix();
            if (tempUnit != null) {
                expStr += "   " + tempUnit.getSymbolUnicode();
            }
        }
        electTable.addCell(createCell(stimName, getFont()));
        electTable.addCell(createCell(currName, getFont()));
        electTable.addCell(createCell(clampType, getFont()));
        electTable.addCell(createCell(expStr, getFont()));
        // add electrode info
        Electrode electrode = electricalStimuli[j].getElectrode();
        if (electrode == null) {
            electTable.addCell(createCell("N/A", getFont()));
        } else {
            Coordinate c = electrode.getPosition();
            String location = c.getX() + ", " + c.getY() + ", " + c.getZ();
            String featureName = electrode.getFeature().getName();
            electTable.addCell(createCell("(" + location + ") in " + featureName, getFont()));
        }
    }
    if (electTable != null) {
        if (memMapSection == null) {
            memMapSection = simContextSection.addSection("Membrane Mapping For " + simContext.getName(), 1);
        }
        electTable.setWidths(widths);
        memMapSection.add(electTable);
    }
    // add temperature
    Table tempTable = getTable(1, 75, 1, 3, 3);
    tempTable.setAlignment(Table.ALIGN_LEFT);
    tempTable.addCell(createCell("Temperature: " + simContext.getTemperatureKelvin() + " K", getFont()));
    if (memMapSection != null) {
        memMapSection.add(tempTable);
    }
}
Also used : MembraneMapping(cbit.vcell.mapping.MembraneMapping) Table(com.lowagie.text.Table) Electrode(cbit.vcell.mapping.Electrode) CurrentDensityClampStimulus(cbit.vcell.mapping.CurrentDensityClampStimulus) Section(com.lowagie.text.Section) StructureMapping(cbit.vcell.mapping.StructureMapping) TotalCurrentClampStimulus(cbit.vcell.mapping.TotalCurrentClampStimulus) LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter) ElectricalStimulus(cbit.vcell.mapping.ElectricalStimulus) VCUnitDefinition(cbit.vcell.units.VCUnitDefinition) FeatureMapping(cbit.vcell.mapping.FeatureMapping) Expression(cbit.vcell.parser.Expression) Coordinate(org.vcell.util.Coordinate) VoltageClampStimulus(cbit.vcell.mapping.VoltageClampStimulus) Parameter(cbit.vcell.model.Parameter) SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter) LocalParameter(cbit.vcell.mapping.ParameterContext.LocalParameter) GeometryContext(cbit.vcell.mapping.GeometryContext)

Aggregations

VCUnitDefinition (cbit.vcell.units.VCUnitDefinition)113 Expression (cbit.vcell.parser.Expression)73 ModelUnitSystem (cbit.vcell.model.ModelUnitSystem)36 ExpressionException (cbit.vcell.parser.ExpressionException)26 PropertyVetoException (java.beans.PropertyVetoException)24 LocalParameter (cbit.vcell.mapping.ParameterContext.LocalParameter)21 Element (org.jdom.Element)20 KineticsParameter (cbit.vcell.model.Kinetics.KineticsParameter)19 ArrayList (java.util.ArrayList)17 ModelParameter (cbit.vcell.model.Model.ModelParameter)16 Model (cbit.vcell.model.Model)14 Parameter (cbit.vcell.model.Parameter)14 SpeciesContext (cbit.vcell.model.SpeciesContext)13 SymbolTableEntry (cbit.vcell.parser.SymbolTableEntry)13 Membrane (cbit.vcell.model.Membrane)12 Structure (cbit.vcell.model.Structure)12 StructureMapping (cbit.vcell.mapping.StructureMapping)11 MathException (cbit.vcell.math.MathException)10 ReactionStep (cbit.vcell.model.ReactionStep)10 Feature (cbit.vcell.model.Feature)9