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

Example 6 with ReactionParticipant

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

the class BioModelEditorReactionTableModel method setValueAt.

public void setValueAt(Object value, int row, int column) {
    if (getModel() == null || value == null) {
        return;
    }
    try {
        ModelProcess modelProcess = getValueAt(row);
        if (modelProcess != null) {
            switch(column) {
                case COLUMN_NAME:
                    {
                        String inputValue = ((String) value);
                        inputValue = inputValue.trim();
                        modelProcess.setName(inputValue);
                        break;
                    }
                case COLUMN_EQUATION:
                    {
                        String inputValue = (String) value;
                        inputValue = inputValue.trim();
                        if (modelProcess instanceof ReactionStep) {
                            ReactionStep reactionStep = (ReactionStep) modelProcess;
                            ReactionParticipant[] rpArray = ModelProcessEquation.parseReaction(reactionStep, getModel(), inputValue);
                            for (ReactionParticipant rp : rpArray) {
                                SpeciesContext speciesContext = rp.getSpeciesContext();
                                if (bioModel.getModel().getSpeciesContext(speciesContext.getName()) == null) {
                                    bioModel.getModel().addSpecies(speciesContext.getSpecies());
                                    bioModel.getModel().addSpeciesContext(speciesContext);
                                }
                            }
                            reactionStep.setReactionParticipants(rpArray);
                        } else if (modelProcess instanceof ReactionRule) {
                            ReactionRule oldReactionRule = (ReactionRule) modelProcess;
                            // when editing an existing reaction rule
                            ReactionRule newReactionRule = (ReactionRule) RbmUtils.parseReactionRule(inputValue, oldReactionRule.getStructure(), bioModel);
                            if (newReactionRule != null) {
                                oldReactionRule.setProductPatterns(newReactionRule.getProductPatterns(), false, false);
                                oldReactionRule.setReactantPatterns(newReactionRule.getReactantPatterns(), false, false);
                            // String name = oldReactionRule.getName();
                            // RbmKineticLaw kl = oldReactionRule.getKineticLaw();
                            // Structure st = oldReactionRule.getStructure();
                            // getModel().getRbmModelContainer().removeReactionRule(oldReactionRule);
                            // newReactionRule.setName(name);
                            // newReactionRule.setKineticLaw(kl);
                            // newReactionRule.setStructure(st);
                            // getModel().getRbmModelContainer().addReactionRule(newReactionRule);
                            }
                        }
                        break;
                    }
                case COLUMN_STRUCTURE:
                    {
                        Structure s = (Structure) value;
                        modelProcess.setStructure(s);
                        break;
                    }
            }
        } else {
            switch(column) {
                case COLUMN_EQUATION:
                    {
                        if (getModel().getNumStructures() == 1) {
                            String inputValue = ((String) value);
                            inputValue = inputValue.trim();
                            if (inputValue.contains("(") && inputValue.contains(")")) {
                                ReactionRule reactionRule = (ReactionRule) RbmUtils.parseReactionRule(inputValue, getModel().getStructure(0), bioModel);
                                getModel().getRbmModelContainer().addReactionRule(reactionRule);
                            } else {
                                if (BioModelEditorRightSideTableModel.ADD_NEW_HERE_REACTION_TEXT.equals(inputValue)) {
                                    return;
                                }
                                ReactionStep reactionStep = getModel().createSimpleReaction(getModel().getStructure(0));
                                ReactionParticipant[] rpArray = ModelProcessEquation.parseReaction(reactionStep, getModel(), inputValue);
                                for (ReactionParticipant rp : rpArray) {
                                    SpeciesContext speciesContext = rp.getSpeciesContext();
                                    if (bioModel.getModel().getSpeciesContext(speciesContext.getName()) == null) {
                                        bioModel.getModel().addSpecies(speciesContext.getSpecies());
                                        bioModel.getModel().addSpeciesContext(speciesContext);
                                    }
                                }
                                reactionStep.setReactionParticipants(rpArray);
                            }
                        }
                        break;
                    }
            }
        }
    } catch (Exception e) {
        e.printStackTrace(System.out);
        DialogUtils.showErrorDialog(ownerTable, e.getMessage(), e);
    }
}

Example 7 with ReactionParticipant

use of cbit.vcell.model.ReactionParticipant 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 8 with ReactionParticipant

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

the class StructureAnalyzer method refreshTotalMatrices.

/**
 * This method was created in VisualAge.
 */
private void refreshTotalMatrices() throws Exception {
    // System.out.println("StructureAnalyzer.refreshTotalMatrices()");
    // 
    // 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();
// }
}

Example 9 with ReactionParticipant

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

the class StructureAnalyzer method refreshFastMatrices.

/**
 * This method was created in VisualAge.
 */
private void refreshFastMatrices() throws Exception {
    // System.out.println("StructureAnalyzer.refreshFastMatrices()");
    // 
    // update scheme matrix for fast system
    // 
    fastSchemeMatrix = new RationalNumberMatrix(fastSpeciesContextMappings.length, fastReactionSteps.length);
    for (int i = 0; i < fastSpeciesContextMappings.length; i++) {
        for (int j = 0; j < fastReactionSteps.length; j++) {
            fastSchemeMatrix.set_elem(i, j, fastReactionSteps[j].getStoichiometry(fastSpeciesContextMappings[i].getSpeciesContext()));
        }
    }
    // 
    for (int i = 0; i < fastSpeciesContextMappings.length; i++) {
        SpeciesContextMapping scm = fastSpeciesContextMappings[i];
        SpeciesContext sc = scm.getSpeciesContext();
        // 
        // collect fast rate expression
        // 
        Expression exp = new Expression(0.0);
        for (int j = 0; j < fastReactionSteps.length; j++) {
            int stoichiometry = fastReactionSteps[j].getStoichiometry(sc);
            ReactionSpec reactionSpec = mathMapping.getSimulationContext().getReactionContext().getReactionSpec(fastReactionSteps[j]);
            if (stoichiometry != 0) {
                if (!reactionSpec.isFast()) {
                    throw new Exception("expected only fast rates");
                }
                if (reactionSpec.isExcluded()) {
                    throw new Exception("expected only included rates");
                }
                ReactionParticipant[] rps = fastReactionSteps[j].getReactionParticipants();
                ReactionParticipant rp0 = null;
                for (ReactionParticipant rp : rps) {
                    if (rp.getSpeciesContext() == sc) {
                        rp0 = rp;
                        break;
                    }
                }
                // 
                if (rp0 != null) {
                    Expression fastRateExpression = getCorrectedRateExpression(fastReactionSteps[j], rp0, RateType.ConcentrationRate).renameBoundSymbols(mathMapping.getNameScope());
                    exp = Expression.add(exp, fastRateExpression);
                }
            }
        }
        // exp.bindExpression(mathMapping);
        scm.setFastRate(exp.flatten());
    }
    // System.out.println("StructureAnalyzer.refreshFastMatrices(), scheme matrix:");
    // fastSchemeMatrix.show();
    // 
    // update null space matrix
    // 
    fastNullSpaceMatrix = fastSchemeMatrix.findNullSpace();
// if (fastNullSpaceMatrix==null){
// System.out.println("fast system has full rank");
// }else{
// System.out.println("StructureAnalyzer.refreshFastMatrices(), nullSpace matrix:");
// fastNullSpaceMatrix.show();
// }
}

Example 10 with ReactionParticipant

use of cbit.vcell.model.ReactionParticipant 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;
}