Examples with ReactionStep cbit.vcell.model.ReactionStep used on opensource projects

Example 36 with ReactionStep

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

the class StochMathMapping method addJumpProcesses.

private void addJumpProcesses(VariableHash varHash, GeometryClass geometryClass, SubDomain subDomain) throws ExpressionException, ModelException, MappingException, MathException {
    // set up jump processes
    // get all the reactions from simulation context
    // ReactionSpec[] reactionSpecs = simContext.getReactionContext().getReactionSpecs();---need to take a look here!
    ModelUnitSystem modelUnitSystem = getSimulationContext().getModel().getUnitSystem();
    ReactionSpec[] reactionSpecs = getSimulationContext().getReactionContext().getReactionSpecs();
    for (ReactionSpec reactionSpec : reactionSpecs) {
        if (reactionSpec.isExcluded()) {
            continue;
        }
        // get the reaction
        ReactionStep reactionStep = reactionSpec.getReactionStep();
        Kinetics kinetics = reactionStep.getKinetics();
        // probability parameter from modelUnitSystem
        VCUnitDefinition probabilityParamUnit = modelUnitSystem.getStochasticSubstanceUnit().divideBy(modelUnitSystem.getTimeUnit());
        // Different ways to deal with simple reactions and flux reactions
        if (// simple reactions
        reactionStep instanceof SimpleReaction) {
            // check the reaction rate law to see if we need to decompose a reaction(reversible) into two jump processes.
            // rate constants are important in calculating the probability rate.
            // for Mass Action, we use KForward and KReverse,
            // for General Kinetics we parse reaction rate J to see if it is in Mass Action form.
            Expression forwardRate = null;
            Expression reverseRate = null;
            if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction) || kinetics.getKineticsDescription().equals(KineticsDescription.General)) {
                Expression rateExp = new Expression(kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate), reactionStep.getNameScope());
                Parameter forwardRateParameter = null;
                Parameter reverseRateParameter = null;
                if (kinetics.getKineticsDescription().equals(KineticsDescription.MassAction)) {
                    forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KForward);
                    reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KReverse);
                }
                MassActionSolver.MassActionFunction 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();
                    }
                }
            } else // if it's macro/microscopic kinetics, we'll have them set up as reactions with only forward rate.
            if (kinetics.getKineticsDescription().equals(KineticsDescription.Macroscopic_irreversible) || kinetics.getKineticsDescription().equals(KineticsDescription.Microscopic_irreversible)) {
                Expression Kon = getIdentifierSubstitutions(new Expression(reactionStep.getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_KOn), getNameScope()), reactionStep.getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_Binding_Radius).getUnitDefinition(), geometryClass);
                if (Kon != null) {
                    Expression KonCopy = new Expression(Kon);
                    try {
                        MassActionSolver.substituteParameters(KonCopy, true).evaluateConstant();
                        forwardRate = new Expression(Kon);
                    } catch (ExpressionException e) {
                        throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Problem with Kon parameter in " + reactionStep.getName() + ":  '" + KonCopy.infix() + "', " + e.getMessage()));
                    }
                } else {
                    throw new MathException(VCellErrorMessages.getMassActionSolverMessage(reactionStep.getName(), "Kon parameter of " + reactionStep.getName() + " is null."));
                }
            }
            boolean isForwardRatePresent = false;
            boolean isReverseRatePresent = false;
            if (forwardRate != null) {
                isForwardRatePresent = true;
            }
            if (reverseRate != null) {
                isReverseRatePresent = true;
            }
            // we process it as forward reaction
            if ((isForwardRatePresent)) /*|| ((forwardRate == null) && (reverseRate == null))*/
            {
                // get jump process name
                String jpName = TokenMangler.mangleToSName(reactionStep.getName());
                // get probability
                Expression exp = null;
                // reactions are of mass action form
                exp = getProbabilityRate(reactionStep, forwardRate, true);
                ProbabilityParameter probParm = null;
                try {
                    probParm = addProbabilityParameter(PARAMETER_PROBABILITYRATE_PREFIX + jpName, exp, PARAMETER_ROLE_P, probabilityParamUnit, reactionStep);
                } catch (PropertyVetoException pve) {
                    pve.printStackTrace();
                    throw new MappingException(pve.getMessage());
                }
                // add probability to function or constant
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(probParm, geometryClass), getIdentifierSubstitutions(exp, probabilityParamUnit, geometryClass), geometryClass));
                JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probParm, geometryClass)));
                // actions
                ReactionParticipant[] reacPart = reactionStep.getReactionParticipants();
                for (int j = 0; j < reacPart.length; j++) {
                    Action action = null;
                    SpeciesCountParameter spCountParam = getSpeciesCountParameter(reacPart[j].getSpeciesContext());
                    if (reacPart[j] instanceof Reactant) {
                        // check if the reactant is a constant. If the species is a constant, there will be no action taken on this species
                        if (// not a constant
                        !simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
                            int stoi = ((Reactant) reacPart[j]).getStoichiometry();
                            action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(-stoi));
                            jp.addAction(action);
                        }
                    } else if (reacPart[j] instanceof Product) {
                        // check if the product is a constant. If the product is a constant, there will be no action taken on this species
                        if (// not a constant
                        !simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
                            int stoi = ((Product) reacPart[j]).getStoichiometry();
                            action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(stoi));
                            jp.addAction(action);
                        }
                    }
                }
                // add jump process to compartment subDomain
                subDomain.addJumpProcess(jp);
            }
            if (// one more jump process for a reversible reaction
            isReverseRatePresent) {
                // get jump process name
                String jpName = TokenMangler.mangleToSName(reactionStep.getName()) + PARAMETER_PROBABILITY_RATE_REVERSE_SUFFIX;
                Expression exp = null;
                // reactions are mass actions
                exp = getProbabilityRate(reactionStep, reverseRate, false);
                ProbabilityParameter probRevParm = null;
                try {
                    probRevParm = addProbabilityParameter(PARAMETER_PROBABILITYRATE_PREFIX + jpName, exp, PARAMETER_ROLE_P_reverse, probabilityParamUnit, reactionStep);
                } catch (PropertyVetoException pve) {
                    pve.printStackTrace();
                    throw new MappingException(pve.getMessage());
                }
                // add probability to function or constant
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(probRevParm, geometryClass), getIdentifierSubstitutions(exp, probabilityParamUnit, geometryClass), geometryClass));
                JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probRevParm, geometryClass)));
                // actions
                ReactionParticipant[] reacPart = reactionStep.getReactionParticipants();
                for (int j = 0; j < reacPart.length; j++) {
                    Action action = null;
                    SpeciesCountParameter spCountParam = getSpeciesCountParameter(reacPart[j].getSpeciesContext());
                    if (reacPart[j] instanceof Reactant) {
                        // check if the reactant is a constant. If the species is a constant, there will be no action taken on this species
                        if (// not a constant
                        !simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
                            int stoi = ((Reactant) reacPart[j]).getStoichiometry();
                            action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(stoi));
                            jp.addAction(action);
                        }
                    } else if (reacPart[j] instanceof Product) {
                        // check if the product is a constant. If the product is a constant, there will be no action taken on this species
                        if (// not a constant
                        !simContext.getReactionContext().getSpeciesContextSpec(reacPart[j].getSpeciesContext()).isConstant()) {
                            int stoi = ((Product) reacPart[j]).getStoichiometry();
                            action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(-stoi));
                            jp.addAction(action);
                        }
                    }
                }
                // add jump process to compartment subDomain
                subDomain.addJumpProcess(jp);
            }
        // end of if(isForwardRateNonZero), if(isReverseRateNonRate)
        } else if (// flux reactions
        reactionStep instanceof FluxReaction) {
            // we could set jump processes for general flux rate in forms of p1*Sout + p2*Sin
            if (kinetics.getKineticsDescription().equals(KineticsDescription.General) || kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
                Expression fluxRate = new Expression(kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate), reactionStep.getNameScope());
                // we have to pass the math description para to flux solver, coz somehow math description in simulation context is not updated.
                // forward and reverse rate parameters may be null
                Parameter forwardRateParameter = null;
                Parameter reverseRateParameter = null;
                if (kinetics.getKineticsDescription().equals(KineticsDescription.GeneralPermeability)) {
                    forwardRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
                    reverseRateParameter = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Permeability);
                }
                MassActionSolver.MassActionFunction fluxFunc = MassActionSolver.solveMassAction(forwardRateParameter, reverseRateParameter, fluxRate, (FluxReaction) reactionStep);
                // create jump process for forward flux if it exists.
                Expression rsStructureSize = new Expression(reactionStep.getStructure().getStructureSize(), getNameScope());
                VCUnitDefinition probRateUnit = modelUnitSystem.getStochasticSubstanceUnit().divideBy(modelUnitSystem.getAreaUnit()).divideBy(modelUnitSystem.getTimeUnit());
                Expression rsRateUnitFactor = getUnitFactor(probRateUnit.divideBy(modelUnitSystem.getFluxReactionUnit()));
                if (fluxFunc.getForwardRate() != null && !fluxFunc.getForwardRate().isZero()) {
                    Expression rate = fluxFunc.getForwardRate();
                    // get species expression (depend on structure, if mem: Species/mem_Size, if vol: species*KMOLE/vol_size)
                    if (fluxFunc.getReactants().size() != 1) {
                        throw new MappingException("Flux " + reactionStep.getName() + " should have only one reactant.");
                    }
                    SpeciesContext scReactant = fluxFunc.getReactants().get(0).getSpeciesContext();
                    Expression scConcExpr = new Expression(getSpeciesConcentrationParameter(scReactant), getNameScope());
                    Expression probExp = Expression.mult(rate, rsRateUnitFactor, rsStructureSize, scConcExpr);
                    // jump process name
                    // +"_reverse";
                    String jpName = TokenMangler.mangleToSName(reactionStep.getName());
                    ProbabilityParameter probParm = null;
                    try {
                        probParm = addProbabilityParameter(PARAMETER_PROBABILITYRATE_PREFIX + jpName, probExp, PARAMETER_ROLE_P, probabilityParamUnit, reactionStep);
                    } catch (PropertyVetoException pve) {
                        pve.printStackTrace();
                        throw new MappingException(pve.getMessage());
                    }
                    // add probability to function or constant
                    String ms = getMathSymbol(probParm, geometryClass);
                    Expression is = getIdentifierSubstitutions(probExp, probabilityParamUnit, geometryClass);
                    Variable nfoc = newFunctionOrConstant(ms, is, geometryClass);
                    varHash.addVariable(nfoc);
                    JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probParm, geometryClass)));
                    // actions
                    Action action = null;
                    SpeciesContext sc = fluxFunc.getReactants().get(0).getSpeciesContext();
                    if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
                        action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(-1));
                        jp.addAction(action);
                    }
                    sc = fluxFunc.getProducts().get(0).getSpeciesContext();
                    if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
                        action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(1));
                        jp.addAction(action);
                    }
                    subDomain.addJumpProcess(jp);
                }
                // create jump process for reverse flux if it exists.
                if (fluxFunc.getReverseRate() != null && !fluxFunc.getReverseRate().isZero()) {
                    // jump process name
                    String jpName = TokenMangler.mangleToSName(reactionStep.getName()) + PARAMETER_PROBABILITY_RATE_REVERSE_SUFFIX;
                    Expression rate = fluxFunc.getReverseRate();
                    // get species expression (depend on structure, if mem: Species/mem_Size, if vol: species*KMOLE/vol_size)
                    if (fluxFunc.getProducts().size() != 1) {
                        throw new MappingException("Flux " + reactionStep.getName() + " should have only one product.");
                    }
                    SpeciesContext scProduct = fluxFunc.getProducts().get(0).getSpeciesContext();
                    Expression scConcExpr = new Expression(getSpeciesConcentrationParameter(scProduct), getNameScope());
                    Expression probExp = Expression.mult(rate, rsRateUnitFactor, rsStructureSize, scConcExpr);
                    ProbabilityParameter probRevParm = null;
                    try {
                        probRevParm = addProbabilityParameter(PARAMETER_PROBABILITYRATE_PREFIX + jpName, probExp, PARAMETER_ROLE_P_reverse, probabilityParamUnit, reactionStep);
                    } catch (PropertyVetoException pve) {
                        pve.printStackTrace();
                        throw new MappingException(pve.getMessage());
                    }
                    // add probability to function or constant
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(probRevParm, geometryClass), getIdentifierSubstitutions(probExp, probabilityParamUnit, geometryClass), geometryClass));
                    JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probRevParm, geometryClass)));
                    // actions
                    Action action = null;
                    SpeciesContext sc = fluxFunc.getReactants().get(0).getSpeciesContext();
                    if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
                        action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(1));
                        jp.addAction(action);
                    }
                    sc = fluxFunc.getProducts().get(0).getSpeciesContext();
                    if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
                        action = Action.createIncrementAction(varHash.getVariable(getMathSymbol(spCountParam, geometryClass)), new Expression(-1));
                        jp.addAction(action);
                    }
                    subDomain.addJumpProcess(jp);
                }
            }
        }
    // end of if (simplereaction)...else if(fluxreaction)
    }
// end of reaction step loop
}

Example 37 with ReactionStep

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

the class RateRule method gatherIssues.

public void gatherIssues(IssueContext issueContext, List<Issue> issueList, Set<String> alreadyIssue) {
    // if we find an issue with the current rule we post it and stop looking for more
    if (fieldName == null || fieldName.isEmpty()) {
        String msg = typeName + " Name is missing";
        issueList.add(new Issue(this, issueContext, IssueCategory.Identifiers, msg, Issue.Severity.ERROR));
        return;
    }
    if (rateRuleVar == null || rateRuleVar.getName() == null || rateRuleVar.getName().isEmpty()) {
        String msg = typeName + " Variable is missing";
        issueList.add(new Issue(this, issueContext, IssueCategory.Identifiers, msg, Issue.Severity.ERROR));
        return;
    }
    if (rateRuleVar instanceof Structure.StructureSize) {
        String msg = Structure.StructureSize.typeName + " Variable is not supported at this time";
        issueList.add(new Issue(this, issueContext, IssueCategory.Identifiers, msg, Issue.Severity.ERROR));
        return;
    }
    if (!(rateRuleVar instanceof SymbolTableEntry)) {
        String msg = typeName + " Variable is not a SymbolTableEntry";
        issueList.add(new Issue(this, issueContext, IssueCategory.Identifiers, msg, Issue.Severity.ERROR));
        return;
    }
    if (rateRuleExpression == null) {
        String msg = typeName + " Expression is missing";
        issueList.add(new Issue(this, issueContext, IssueCategory.Identifiers, msg, Issue.Severity.ERROR));
        return;
    }
    if (rateRuleExpression != null) {
        String[] symbols = rateRuleExpression.getSymbols();
        if (symbols != null && symbols.length > 0) {
            for (String symbol : symbols) {
                SymbolTableEntry ste = simulationContext.getEntry(symbol);
                if (ste == null) {
                    String msg = "Missing Symbol '" + symbol + "' in Expression for " + typeName + " '" + fieldName + "'.";
                    issueList.add(new Issue(this, issueContext, IssueCategory.Identifiers, msg, Issue.Severity.ERROR));
                    return;
                }
            }
        }
    }
    // TODO: should also look into interaction with events!
    if (simulationContext.getRateRules() != null) {
        for (RateRule rr : simulationContext.getRateRules()) {
            if (rr == this) {
                continue;
            }
            if (rr.getRateRuleVar() == null) {
                // another rate rule variable may be still not initialized
                continue;
            }
            String ruleVariableName = rateRuleVar.getName();
            if (!alreadyIssue.contains(ruleVariableName) && ruleVariableName.equals(rr.getRateRuleVar().getName())) {
                String msg = typeName + " Variable '" + rateRuleVar.getName() + " is duplicated in " + rr.getDisplayName();
                issueList.add(new Issue(this, issueContext, IssueCategory.Identifiers, msg, Issue.Severity.ERROR));
                alreadyIssue.add(ruleVariableName);
                return;
            }
        }
    }
    // we do the following check for the assignemnt rule too, so there's no point in complaining here too about the same problem
    // TODO: uncomment to raise an issue here as well
    // if(simulationContext.getAssignmentRules() != null) {
    // for(AssignmentRule ar : simulationContext.getAssignmentRules()) {
    // if(ar.getAssignmentRuleVar() == null) {
    // continue;		// we may be in the middle of creating the assignment rule and the variable is still missing
    // }
    // if(rateRuleVar.getName().equals(ar.getAssignmentRuleVar().getName())) {
    // String msg = typeName + " Variable '" + rateRuleVar.getName() + "' is duplicated as an " + AssignmentRule.typeName + " Variable";
    // issueList.add(new Issue(this, issueContext, IssueCategory.Identifiers, msg, Issue.Severity.ERROR));
    // return;
    // }
    // }
    // }
    // TODO: the following code belongs to SpeciesContextSpec rather than here (so that the Issue navigates to the Specifications / Species panel
    // rate rule variable must not be also a reactant or a product in a reaction
    ReactionContext rc = getSimulationContext().getReactionContext();
    ReactionSpec[] rsArray = rc.getReactionSpecs();
    if (rateRuleVar instanceof SpeciesContext) {
        for (ReactionSpec rs : rsArray) {
            if (rs.isExcluded()) {
                // we don't care about reactions which are excluded
                continue;
            }
            ReactionStep step = rs.getReactionStep();
            for (ReactionParticipant p : step.getReactionParticipants()) {
                if (p instanceof Product || p instanceof Reactant) {
                    SpeciesContext sc = p.getSpeciesContext();
                    SpeciesContextSpec scs = rc.getSpeciesContextSpec(sc);
                    if (!scs.isConstant() && sc.getName().equals(rateRuleVar.getName())) {
                        String msg = "'" + rateRuleVar.getName() + "' must be clamped to be both a rate rule variable and a participant in reaction '" + step.getDisplayName() + "'.";
                        issueList.add(new Issue(this, issueContext, IssueCategory.Identifiers, msg, msg, Issue.Severity.ERROR));
                        return;
                    }
                }
            }
        }
    }
    if (sbmlName != null && sbmlName.isEmpty()) {
        String message = "SbmlName cannot be an empty string.";
        issueList.add(new Issue(this, issueContext, IssueCategory.Identifiers, message, Issue.Severity.ERROR));
    }
    if (sbmlId != null && sbmlId.isEmpty()) {
        String message = "SbmlId cannot be an empty string.";
        issueList.add(new Issue(this, issueContext, IssueCategory.Identifiers, message, Issue.Severity.ERROR));
    }
}

Example 38 with ReactionStep

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

the class SimulationContext method gatherIssues.

public void gatherIssues(IssueContext issueContext, List<Issue> issueVector, boolean bIgnoreMathDescription) {
    // issueContext = issueContext.newChildContext(ContextType.SimContext, this);
    if (applicationType.equals(Application.RULE_BASED_STOCHASTIC)) {
        for (ReactionRuleSpec rrs : getReactionContext().getReactionRuleSpecs()) {
            if (rrs.isExcluded()) {
                continue;
            }
            ReactionRule rr = rrs.getReactionRule();
            if (rr.getReactantPatterns().size() > 2) {
                String message = "NFSim doesn't support more than 2 reactants within a reaction rule.";
                issueVector.add(new Issue(rr, issueContext, IssueCategory.Identifiers, message, Issue.Severity.WARNING));
            }
            if (rr.isReversible() && rr.getProductPatterns().size() > 2) {
                String message = "NFSim doesn't support more than 2 products within a reversible reaction rule.";
                issueVector.add(new Issue(rr, issueContext, IssueCategory.Identifiers, message, Issue.Severity.WARNING));
            }
        }
        for (ReactionSpec rrs : getReactionContext().getReactionSpecs()) {
            if (rrs.isExcluded()) {
                continue;
            }
            ReactionStep rs = rrs.getReactionStep();
            if (rs.getNumReactants() > 2) {
                String message = "NFSim doesn't support more than 2 reactants within a reaction step.";
                issueVector.add(new Issue(rs, issueContext, IssueCategory.Identifiers, message, Issue.Severity.WARNING));
            }
            if (rs.isReversible() && rs.getNumProducts() > 2) {
                String message = "NFSim doesn't support more than 2 products within a reversible reaction step.";
                issueVector.add(new Issue(rs, issueContext, IssueCategory.Identifiers, message, Issue.Severity.WARNING));
            }
        }
        // we give warning when we have plain reactions with participants with patterns;
        // making rules from these may result in inconsistent interpretation for the constant rates
        boolean isParticipantWithPattern = false;
        for (ReactionSpec rrs : getReactionContext().getReactionSpecs()) {
            if (rrs.isExcluded()) {
                continue;
            }
            ReactionStep rs = rrs.getReactionStep();
            for (Reactant r : rs.getReactants()) {
                if (r.getSpeciesContext().hasSpeciesPattern()) {
                    isParticipantWithPattern = true;
                    break;
                }
            }
            if (isParticipantWithPattern) {
                break;
            }
            for (Product p : rs.getProducts()) {
                if (p.getSpeciesContext().hasSpeciesPattern()) {
                    isParticipantWithPattern = true;
                    break;
                }
            }
            if (isParticipantWithPattern) {
                break;
            }
        }
        if (isParticipantWithPattern) {
            String message = SimulationContext.rateWarning2;
            String tooltip = SimulationContext.rateWarning;
            issueVector.add(new Issue(this, issueContext, IssueCategory.Identifiers, message, tooltip, Issue.Severity.WARNING));
        }
        for (Structure struct : getModel().getStructures()) {
            String name = struct.getName();
            if (!name.equals(TokenMangler.fixTokenStrict(name))) {
                String msg = "'" + name + "' not legal identifier for rule-based stochastic applications, try '" + TokenMangler.fixTokenStrict(name) + "'.";
                issueVector.add(new Issue(struct, issueContext, IssueCategory.Identifiers, msg, Issue.Severity.ERROR));
            }
        }
    }
    if (fieldBioEvents != null) {
        for (BioEvent bioEvent : fieldBioEvents) {
            bioEvent.gatherIssues(issueContext, issueVector);
        }
    }
    if (spatialObjects != null) {
        for (SpatialObject spatialObject : spatialObjects) {
            spatialObject.gatherIssues(issueContext, issueVector);
        }
    }
    if (spatialProcesses != null) {
        for (SpatialProcess spatialProcess : spatialProcesses) {
            spatialProcess.gatherIssues(issueContext, issueVector);
        }
    }
    if (fieldRateRules != null) {
        // to avoid duplicated Issues for the same problem
        Set<String> alreadyIssue = new HashSet<>();
        for (RateRule rr : fieldRateRules) {
            rr.gatherIssues(issueContext, issueVector, alreadyIssue);
        }
    }
    if (fieldAssignmentRules != null) {
        Set<String> alreadyIssue = new HashSet<>();
        for (AssignmentRule ar : fieldAssignmentRules) {
            ar.gatherIssues(issueContext, issueVector, alreadyIssue);
        }
    }
    if (applicationType.equals(Application.NETWORK_DETERMINISTIC) && getModel().getRbmModelContainer().getMolecularTypeList().size() > 0) {
        // we're going to use network transformer to flatten (or we already did)
        if (isInsufficientIterations()) {
            issueVector.add(new Issue(this, issueContext, IssueCategory.RbmNetworkConstraintsBad, IssueInsufficientIterations, Issue.Severity.WARNING));
        }
        if (isInsufficientMaxMolecules()) {
            issueVector.add(new Issue(this, issueContext, IssueCategory.RbmNetworkConstraintsBad, IssueInsufficientMolecules, Issue.Severity.WARNING));
        }
    }
    Geometry geo = getGeometryContext().getGeometry();
    int dimension = geo.getDimension();
    if (dimension != 0) {
        for (ReactionSpec rrs : getReactionContext().getReactionSpecs()) {
            if (rrs.isExcluded()) {
                continue;
            }
            ReactionStep rs = rrs.getReactionStep();
            if (rs.getStructure() instanceof Membrane) {
                continue;
            }
            // for spatial applications
            // we look for reactions where reactants and products are in more than one compartments
            // if such reactions are not on a membrane, we issue a warning
            Set<Structure> features = new HashSet<>();
            for (Reactant r : rs.getReactants()) {
                Structure struct = r.getStructure();
                if (struct instanceof Feature) {
                    features.add(struct);
                }
            }
            for (Product p : rs.getProducts()) {
                Structure struct = p.getStructure();
                if (struct instanceof Feature) {
                    features.add(struct);
                }
            }
            if (features.size() > 1) {
                String message = "Reaction must be situated on a membrane (spatial application present)";
                String tooltip = "Spatial application '" + getName() + "' requires that reactions between compartments must be situated on a membrane.";
                issueVector.add(new Issue(rs, issueContext, IssueCategory.Identifiers, message, tooltip, Issue.Severity.WARNING));
            }
        }
    }
    getReactionContext().gatherIssues(issueContext, issueVector);
    getGeometryContext().gatherIssues(issueContext, issueVector);
    if (fieldAnalysisTasks != null) {
        for (AnalysisTask analysisTask : fieldAnalysisTasks) {
            analysisTask.gatherIssues(issueContext, issueVector);
        }
    }
    getOutputFunctionContext().gatherIssues(issueContext, issueVector);
    getMicroscopeMeasurement().gatherIssues(issueContext, issueVector);
    if (getMathDescription() != null && !bIgnoreMathDescription) {
        getMathDescription().gatherIssues(issueContext, issueVector);
    }
    if (networkConstraints == null) {
    // issueVector.add(new Issue(this, issueContext, IssueCategory.RbmNetworkConstraintsBad, "Network Constraints is null", Issue.Severity.ERROR));
    } else {
        networkConstraints.gatherIssues(issueContext, issueVector);
    }
}

Example 39 with ReactionStep

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

the class StructureAnalyzer method refreshTotalSpeciesContextMappings.

/**
 * This method was created in VisualAge.
 */
private void refreshTotalSpeciesContextMappings() throws java.beans.PropertyVetoException {
    if (structures == null) {
        return;
    }
    // System.out.println("StructureAnalyzer.refreshSpeciesContextMappings()");
    // GeometryContext geoContext = mathMapping.getSimulationContext().getGeometryContext();
    SimulationContext simContext = mathMapping.getSimulationContext();
    Model model = simContext.getReactionContext().getModel();
    // 
    // note, the order of species is specified such that the first species have priority
    // when the null space is solved for dependent variables.  So the order of priority
    // for elimination are as follows:
    // 
    // 1) Species involved with fast reactions.
    // 2) Species not involved with fast reactions.
    // 
    Vector<SpeciesContextMapping> scmList = new Vector<SpeciesContextMapping>();
    // 
    for (int i = 0; i < structures.length; i++) {
        SpeciesContext[] speciesContexts = model.getSpeciesContexts(structures[i]);
        for (int j = 0; j < speciesContexts.length; j++) {
            SpeciesContext sc = speciesContexts[j];
            SpeciesContextMapping scm = mathMapping.getSpeciesContextMapping(sc);
            SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(sc);
            if (scm.isFastParticipant() && !scs.isConstant()) {
                scmList.addElement(scm);
            }
        }
    }
    // 
    for (int i = 0; i < structures.length; i++) {
        SpeciesContext[] speciesContexts = model.getSpeciesContexts(structures[i]);
        for (int j = 0; j < speciesContexts.length; j++) {
            SpeciesContext sc = speciesContexts[j];
            SpeciesContextMapping scm = mathMapping.getSpeciesContextMapping(sc);
            SpeciesContextSpec scs = simContext.getReactionContext().getSpeciesContextSpec(sc);
            if (!scm.isFastParticipant() && !scs.isConstant()) {
                scmList.addElement(scm);
            }
        }
    }
    if (scmList.size() > 0) {
        speciesContextMappings = new SpeciesContextMapping[scmList.size()];
        scmList.copyInto(speciesContextMappings);
        for (int i = 0; i < speciesContextMappings.length; i++) {
            speciesContextMappings[i].setRate(new Expression(0.0));
        // System.out.println("speciesContextMappings["+i+"] = "+speciesContextMappings[i].getSpeciesContext().getName());
        }
    } else {
        speciesContextMappings = null;
    }
    // System.out.println("StructureAnalyzer.refreshTotalSpeciesContextMapping(), speciesContextMappings.length = "+scmList.size());
    // 
    // get all reactionSteps associated with these structures
    // 
    Vector<ReactionStep> rsList = new Vector<ReactionStep>();
    ReactionSpec[] allReactionSpecs = simContext.getReactionContext().getReactionSpecs();
    for (int i = 0; i < allReactionSpecs.length; i++) {
        if (allReactionSpecs[i].isExcluded()) {
            continue;
        }
        ReactionStep rs = allReactionSpecs[i].getReactionStep();
        for (int j = 0; j < structures.length; j++) {
            if (rs.getStructure() == structures[j]) {
                rsList.addElement(rs);
            }
        }
    }
    // 
    for (int i = 0; i < scmList.size(); i++) {
        SpeciesContextMapping scm = (SpeciesContextMapping) scmList.elementAt(i);
        if (!simContext.isUsingMassConservationModelReduction()) {
            // 
            // break mass conservation on all species (disables model reduction).
            // 
            rsList.addElement(new DisableModelReductionReactionStep("DisableModelReductionReactionStep" + i, model, scm.getSpeciesContext().getStructure(), scm.getSpeciesContext()));
        } else {
            // 
            if (scm.isPDERequired() || simContext.getReactionContext().getSpeciesContextSpec(scm.getSpeciesContext()).isWellMixed()) {
                rsList.addElement(new DiffusionDummyReactionStep("DiffusionDummyReactionStep" + i, model, scm.getSpeciesContext().getStructure(), scm.getSpeciesContext()));
            }
            if (scm.hasEventAssignment()) {
                rsList.addElement(new EventDummyReactionStep("EventDummyReactionStep" + i, model, scm.getSpeciesContext().getStructure(), scm.getSpeciesContext()));
            }
            if (scm.hasHybridReaction()) {
                rsList.addElement(new HybridDummyReactionStep("HybridDummyReactionStep" + i, model, scm.getSpeciesContext().getStructure(), scm.getSpeciesContext()));
            }
            if (simContext.isStoch() && simContext.getGeometry().getDimension() > 0 && !simContext.getReactionContext().getSpeciesContextSpec(scm.getSpeciesContext()).isForceContinuous()) {
                rsList.addElement(new ParticleDummyReactionStep("ParticleDummyReactionStep" + i, model, scm.getSpeciesContext().getStructure(), scm.getSpeciesContext()));
            }
        }
    }
    if (rsList.size() > 0) {
        reactionSteps = new ReactionStep[rsList.size()];
        rsList.copyInto(reactionSteps);
    } else {
        reactionSteps = null;
    }
// System.out.println("StructureAnalyzer.refreshTotalSpeciesContextMapping(), reactionSteps.length = "+scmList.size());
}

Example 40 with ReactionStep

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

the class StructureAnalyzer method refreshTotalMatrices.

/**
 * This method was created in VisualAge.
 */
private void refreshTotalMatrices() throws Exception {
    // 
    if (reactionSteps == null) {
        totalNullSpaceMatrix = new RationalNumberMatrix(speciesContextMappings.length, speciesContextMappings.length);
        totalNullSpaceMatrix.identity();
        return;
    }
    // 
    // update scheme matrix for full system (slow and fast)
    // 
    ReactionSpec[] reactionSpecs = new ReactionSpec[reactionSteps.length];
    for (int j = 0; j < reactionSteps.length; j++) {
        reactionSpecs[j] = mathMapping.getSimulationContext().getReactionContext().getReactionSpec(reactionSteps[j]);
    }
    // 
    // initialize rate expressions for speciesContext's due to scheme matrix
    // 
    totalSchemeMatrix = new RationalNumberMatrix(speciesContextMappings.length, reactionSteps.length);
    for (int i = 0; i < speciesContextMappings.length; i++) {
        SpeciesContextMapping scm = speciesContextMappings[i];
        SpeciesContext sc = scm.getSpeciesContext();
        // 
        // collect slow rate expression (fast handled by FastSystem)
        // 
        Expression exp = new Expression(0.0);
        for (int j = 0; j < reactionSteps.length; j++) {
            ReactionStep reactionStep = reactionSteps[j];
            int stoichiometry = reactionStep.getStoichiometry(sc);
            totalSchemeMatrix.set_elem(i, j, stoichiometry);
            if (stoichiometry != 0) {
                if (!(reactionStep instanceof DummyReactionStep) && !reactionSpecs[j].isFast() && !reactionSpecs[j].isExcluded()) {
                    ReactionParticipant[] rps1 = reactionStep.getReactionParticipants();
                    ReactionParticipant rp0 = null;
                    for (ReactionParticipant rp : rps1) {
                        if (rp.getSpeciesContext() == sc) {
                            rp0 = rp;
                            break;
                        }
                    }
                    if (rp0 != null) {
                        Expression distributedReactRateExp = getCorrectedRateExpression(reactionStep, rp0, RateType.ConcentrationRate);
                        exp = Expression.add(exp, distributedReactRateExp);
                    }
                }
            }
        }
        // exp.bindExpression(mathMapping);
        scm.setRate(exp.flatten());
    }
    // 
    if (totalSchemeMatrix.getNumRows() > 1) {
        totalNullSpaceMatrix = (RationalMatrix) totalSchemeMatrix.findNullSpace();
    } else {
        totalNullSpaceMatrix = null;
    }
// if (totalNullSpaceMatrix==null){
// System.out.println("total system has full rank");
// }else{
// System.out.println("StructureAnalyzer.refreshTotalMatrices(), nullSpace matrix:");
// totalNullSpaceMatrix.show();
// }
}