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

Example 11 with Reactant

use of cbit.vcell.model.Reactant 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 12 with Reactant

use of cbit.vcell.model.Reactant 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 13 with Reactant

use of cbit.vcell.model.Reactant 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 14 with Reactant

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

the class RbmNetworkGenerator method generateModel.

public static void generateModel(BioModel bioModel, String netfile) throws Exception {
    Model model = bioModel.getModel();
    Map<String, SpeciesContext> speciesMap = new HashMap<String, SpeciesContext>();
    Map<String, ReactionStep> reactionMap = new HashMap<String, ReactionStep>();
    List<ReactionLine> reactionLineList = new ArrayList<ReactionLine>();
    BufferedReader br = new BufferedReader(new StringReader(netfile));
    int reversibleCount = 0;
    int reactionCount = 0;
    while (true) {
        String line = br.readLine();
        if (line == null) {
            break;
        }
        line = line.trim();
        if (line.equals(BEGIN_PARAMETERS)) {
            while (true) {
                String line2 = br.readLine();
                line2 = line2.trim();
                if (line2.length() == 0) {
                    continue;
                }
                if (line2.equals(END_PARAMETERS)) {
                    break;
                }
                StringTokenizer st = new StringTokenizer(line2);
                String token1 = st.nextToken();
                String token2 = st.nextToken();
                String token3 = st.nextToken();
                ModelParameter mp = model.new ModelParameter(token2, new Expression(token3), Model.ROLE_UserDefined, bioModel.getModel().getUnitSystem().getInstance_TBD());
                model.addModelParameter(mp);
            }
        } else if (line.equals(BEGIN_SPECIES)) {
            while (true) {
                String line2 = br.readLine();
                line2 = line2.trim();
                if (line2.length() == 0) {
                    continue;
                }
                if (line2.equals(END_SPECIES)) {
                    break;
                }
                StringTokenizer st = new StringTokenizer(line2);
                // no
                String token1 = st.nextToken();
                // pattern
                String token2 = st.nextToken();
                // initial condition
                String token3 = st.nextToken();
                String newname = token2.replaceAll("\\.", "_");
                newname = newname.replaceAll("[\\(,][a-zA-Z]\\w*", "");
                newname = newname.replaceAll("~|!\\d*", "");
                newname = newname.replaceAll("\\(\\)", "");
                newname = newname.replaceAll("\\)", "");
                SpeciesContext sc = model.createSpeciesContext(model.getStructure(0));
                sc.setName(newname);
                bioModel.getVCMetaData().setFreeTextAnnotation(sc, token2);
                bioModel.getVCMetaData().setFreeTextAnnotation(sc.getSpecies(), token2);
                speciesMap.put(token1, sc);
            }
        } else if (line.equals(BEGIN_REACTIONS)) {
            while (true) {
                String line2 = br.readLine();
                line2 = line2.trim();
                if (line2.length() == 0) {
                    continue;
                }
                if (line2.equals(END_REACTIONS)) {
                    break;
                }
                ++reactionCount;
                StringTokenizer st = new StringTokenizer(line2);
                String token1 = st.nextToken();
                // reactants
                String token2 = st.nextToken();
                // products
                String token3 = st.nextToken();
                // rate
                String token4 = st.nextToken();
                String token5 = st.nextToken();
                boolean bFoundReversible = false;
                Expression rate = new Expression(token4);
                for (ReactionLine rl : reactionLineList) {
                    if (token2.equals(rl.products) && token3.equals(rl.reactants) && token5.equals(rl.ruleLabel + "r")) {
                        ReactionStep rs = reactionMap.get(rl.no);
                        ((MassActionKinetics) rs.getKinetics()).getReverseRateParameter().setExpression(rate);
                        reactionLineList.remove(rl);
                        bFoundReversible = true;
                        break;
                    }
                }
                if (bFoundReversible) {
                    ++reversibleCount;
                    continue;
                }
                ReactionLine rl = new ReactionLine(token1, token2, token3, token5);
                reactionLineList.add(rl);
                SimpleReaction reaction = model.createSimpleReaction(model.getStructure(0));
                reactionMap.put(token1, reaction);
                reaction.setModel(model);
                bioModel.getVCMetaData().setFreeTextAnnotation(reaction, line2);
                MassActionKinetics kinetics = new MassActionKinetics(reaction);
                reaction.setKinetics(kinetics);
                st = new StringTokenizer(token2, ",");
                while (st.hasMoreTokens()) {
                    String t = st.nextToken();
                    SpeciesContext sc = speciesMap.get(t);
                    if (sc != null) {
                        boolean bExists = false;
                        for (ReactionParticipant rp : reaction.getReactionParticipants()) {
                            if (rp instanceof Reactant && rp.getSpeciesContext() == sc) {
                                rp.setStoichiometry(rp.getStoichiometry() + 1);
                                bExists = true;
                                break;
                            }
                        }
                        if (!bExists) {
                            reaction.addReactant(sc, 1);
                        }
                    }
                }
                st = new StringTokenizer(token3, ",");
                while (st.hasMoreTokens()) {
                    String t = st.nextToken();
                    SpeciesContext sc = speciesMap.get(t);
                    if (sc != null) {
                        boolean bExists = false;
                        for (ReactionParticipant rp : reaction.getReactionParticipants()) {
                            if (rp instanceof Product && rp.getSpeciesContext() == sc) {
                                rp.setStoichiometry(rp.getStoichiometry() + 1);
                                bExists = true;
                                break;
                            }
                        }
                        if (!bExists) {
                            reaction.addProduct(sc, 1);
                        }
                    }
                }
                kinetics.getForwardRateParameter().setExpression(rate);
            }
        }
    }
    System.out.println(model.getNumSpecies() + " species added");
    System.out.println(model.getNumReactions() + " reactions added");
    System.out.println(reversibleCount + " reversible reactions found");
    if (reactionCount != model.getNumReactions() + reversibleCount) {
        throw new RuntimeException("Reactions are not imported correctly!");
    }
}

Example 15 with Reactant

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

the class ReactionParticipantShape method getCurve.

@Override
protected final CubicCurve2D.Double getCurve() {
    // TODO is this the best place for layout?
    refreshLayoutSelf();
    // default behavior of control points is for direction at ends to follow secant between end-points.
    if (lastCurve_Start == null || !lastCurve_Start.equals(start) || lastCurve_End == null || !lastCurve_End.equals(end)) {
        lastp1ctrl = new Point2D.Double((1.0 - FRACT_WEIGHT) * start.getX() + FRACT_WEIGHT * end.getX(), (1.0 - FRACT_WEIGHT) * start.getY() + FRACT_WEIGHT * end.getY());
    }
    Point2D.Double p2ctrl = new Point2D.Double(FRACT_WEIGHT * start.getX() + (1.0 - FRACT_WEIGHT) * end.getX(), FRACT_WEIGHT * start.getY() + (1.0 - FRACT_WEIGHT) * end.getY());
    // check for siblings in the reaction (like a+a-> or a->a) and calculate a correction factor
    // so that the edges won't overlap
    correctionFactor = 0;
    if (endShape instanceof ReactionStepShape) {
        // index in the array of siblings
        int myPosition = 0;
        // including myself
        int numSiblingsReactant = 0;
        int numSiblingsProduct = 0;
        // in total
        int numSiblings = 0;
        int numOthers = 0;
        ReactionStepShape reactionStepShape = (ReactionStepShape) endShape;
        ReactionStep rs = reactionStepShape.getReactionStep();
        for (ReactionParticipant rp : rs.getReactionParticipants()) {
            if (rp == reactionParticipant) {
                // myself
                if (rp instanceof Reactant) {
                    myPosition = numSiblings;
                    numSiblingsReactant++;
                    numSiblings++;
                } else {
                    myPosition = numSiblings;
                    numSiblingsProduct++;
                    numSiblings++;
                }
            } else if (rp.getSpeciesContext().getName().equals(reactionParticipant.getSpeciesContext().getName())) {
                if (rp instanceof Reactant) {
                    numSiblingsReactant++;
                    numSiblings++;
                } else {
                    numSiblingsProduct++;
                    numSiblings++;
                }
            } else {
                numOthers++;
            }
        }
        if (numSiblings > 1) {
            if (!(numSiblingsReactant == 1 && numSiblingsProduct == 1 && numOthers >= 1)) {
                double offset = numSiblings / 2 - 0.5;
                correctionFactor = (myPosition - offset) * 0.08;
            } else {
                // TODO: comment out the next 2 lines to have the "old style" behavior for a+b->a
                double offset = numSiblings / 2 - 0.5;
                correctionFactor = (myPosition - offset) * 0.08;
            }
        }
    }
    // calculate tangent direction at "reactionStep"
    double tangentX = 0.0;
    double tangentY = 0.0;
    if (endShape instanceof ReactionStepShape) {
        ReactionStepShape reactionStepShape = (ReactionStepShape) endShape;
        for (Shape shape : graphModel.getShapes()) {
            if (shape instanceof ReactionParticipantShape && ((ReactionParticipantShape) shape).endShape == reactionStepShape) {
                ReactionParticipantShape rpShape = (ReactionParticipantShape) shape;
                double dx = rpShape.start.getX() - rpShape.end.getX();
                double dy = rpShape.start.getY() - rpShape.end.getY();
                double len = dx * dx + dy * dy;
                if (shape instanceof ProductShape) {
                    ProductShape ps = (ProductShape) shape;
                    tangentX += (ps.start.getX() - ps.end.getX()) / len;
                    tangentY += (ps.start.getY() - ps.end.getY()) / len;
                } else if (shape instanceof ReactantShape) {
                    ReactantShape rs = (ReactantShape) shape;
                    tangentX -= (rs.start.getX() - rs.end.getX()) / len;
                    tangentY -= (rs.start.getY() - rs.end.getY()) / len;
                }
            }
        }
    }
    double tangentLength = Math.sqrt(tangentX * tangentX + tangentY * tangentY);
    if (tangentLength != 0) {
        tangentX = tangentX * CONTROL_WEIGHT / tangentLength;
        tangentY = tangentY * CONTROL_WEIGHT / tangentLength;
    }
    // tangentY = 0.0;
    if (this instanceof CatalystShape) {
        // choose side based on inner product with displacement vector between catalyst and reactionStep
        if (((start.getX() - end.getX()) * tangentY - (start.getY() - end.getY()) * tangentX) > 0) {
            p2ctrl.setLocation(end.getX() + tangentY, end.getY() - tangentX);
        } else {
            p2ctrl.setLocation(end.getX() - tangentY, end.getY() + tangentX);
        }
    } else if (this instanceof ProductShape) {
        p2ctrl.setLocation(end.getX() + tangentX, end.getY() + tangentY);
    } else if (this instanceof ReactantShape) {
        p2ctrl.setLocation(end.getX() - tangentX, end.getY() - tangentY);
    }
    if (lastCurve != null && lastCurve_Start != null && lastCurve_Start.equals(start) && lastCurve_End != null && lastCurve_End.equals(end) && lastp2ctrl != null && lastp2ctrl.equals(p2ctrl)) {
    // Do Nothing
    } else {
        lastCurve = new CubicCurve2D.Double(start.getX(), start.getY(), lastp1ctrl.getX() * (1 + correctionFactor), lastp1ctrl.getY() * (1 + correctionFactor), p2ctrl.getX(), p2ctrl.getY(), end.getX(), end.getY());
        lastCurve_Start = new Point(start);
        lastCurve_End = new Point(end);
        lastp2ctrl = p2ctrl;
    }
    return lastCurve;
}