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

Example 76 with ReactionStep

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

the class StochMathMapping_4_8 method getProbabilityRate.

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

Example 77 with ReactionStep

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

the class MathMapping_4_8 method refreshSpeciesContextMappings.

/**
 * This method was created in VisualAge.
 */
private void refreshSpeciesContextMappings() throws ExpressionException, MappingException, MathException {
    // 
    // create a SpeciesContextMapping for each speciesContextSpec.
    // 
    // set initialExpression from SpeciesContextSpec.
    // set diffusing
    // set variable (only if "Constant" or "Function", else leave it as null)
    // 
    speciesContextMappingList.removeAllElements();
    SpeciesContextSpec[] speciesContextSpecs = simContext.getReactionContext().getSpeciesContextSpecs();
    for (int i = 0; i < speciesContextSpecs.length; i++) {
        SpeciesContextSpec scs = speciesContextSpecs[i];
        SpeciesContextMapping scm = new SpeciesContextMapping(scs.getSpeciesContext());
        scm.setPDERequired(simContext.isPDERequired(scs.getSpeciesContext()));
        scm.setHasEventAssignment(simContext.hasEventAssignment(scs.getSpeciesContext()));
        scm.setHasHybridReaction(false);
        for (ReactionSpec reactionSpec : getSimulationContext().getReactionContext().getReactionSpecs()) {
            if (!reactionSpec.isExcluded() && reactionSpec.hasHybrid(getSimulationContext(), scs.getSpeciesContext())) {
                scm.setHasHybridReaction(true);
            }
        }
        // scm.setAdvecting(isAdvectionRequired(scs.getSpeciesContext()));
        if (scs.isConstant()) {
            Expression initCond = scs.getInitialConditionParameter() == null ? null : scs.getInitialConditionParameter().getExpression();
            scm.setDependencyExpression(initCond);
        // //
        // // determine if a Function is necessary
        // //
        // boolean bNeedFunction = false;
        // if (initCond.getSymbols()!=null){
        // bNeedFunction = true;
        // }
        // if (bNeedFunction){
        // scm.setVariable(new Function(scm.getSpeciesContext().getName(),initCond));
        // }else{
        // scm.setVariable(new Constant(scm.getSpeciesContext().getName(),initCond));
        // }
        }
        // 
        // test if participant in fast reaction step, request elimination if possible
        // 
        scm.setFastParticipant(false);
        ReactionSpec[] reactionSpecs = simContext.getReactionContext().getReactionSpecs();
        for (int j = 0; j < reactionSpecs.length; j++) {
            ReactionSpec reactionSpec = reactionSpecs[j];
            if (reactionSpec.isExcluded()) {
                continue;
            }
            ReactionStep rs = reactionSpec.getReactionStep();
            if (rs instanceof SimpleReaction && rs.countNumReactionParticipants(scs.getSpeciesContext()) > 0) {
                if (reactionSpec.isFast()) {
                    scm.setFastParticipant(true);
                }
            }
        }
        speciesContextMappingList.addElement(scm);
    }
}

Example 78 with ReactionStep

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

the class MembraneStructureAnalyzer method refreshResolvedFluxes.

/**
 * This method was created in VisualAge.
 */
void refreshResolvedFluxes() throws Exception {
    // System.out.println("MembraneStructureAnalyzer.refreshResolvedFluxes()");
    GeometryContext geoContext = mathMapping_4_8.getSimulationContext().getGeometryContext();
    StructureTopology structTopology = mathMapping_4_8.getSimulationContext().getModel().getStructureTopology();
    Vector<ResolvedFlux> resolvedFluxList = new Vector<ResolvedFlux>();
    // 
    // for each reaction, get all fluxReactions associated with this membrane
    // 
    Vector<ReactionStep> fluxList = new Vector<ReactionStep>();
    ReactionSpec[] reactionSpecs = mathMapping_4_8.getSimulationContext().getReactionContext().getReactionSpecs();
    for (int j = 0; j < reactionSpecs.length; j++) {
        if (reactionSpecs[j].isExcluded()) {
            continue;
        }
        ReactionStep rs = reactionSpecs[j].getReactionStep();
        if (rs.getStructure() == getMembrane()) {
            if (rs instanceof FluxReaction) {
                fluxList.addElement(rs);
            }
        }
    }
    // 
    for (int i = 0; i < fluxList.size(); i++) {
        FluxReaction fr = (FluxReaction) fluxList.elementAt(i);
        Species fluxCarrier = null;
        for (ReactionParticipant rp : fr.getReactionParticipants()) {
            if (rp instanceof Reactant || rp instanceof Product) {
                if (fluxCarrier == null) {
                    fluxCarrier = rp.getSpecies();
                } else {
                    if (fluxCarrier != rp.getSpecies()) {
                        throw new Exception("Flux reaction '" + fr.getName() + "' with multiple species not allowed in VCell 4.8.");
                    }
                }
            }
        }
        if (fluxCarrier == null) {
            continue;
        }
        ResolvedFlux rf = null;
        for (int j = 0; j < resolvedFluxList.size(); j++) {
            ResolvedFlux rf_tmp = (ResolvedFlux) resolvedFluxList.elementAt(j);
            if (rf_tmp.getSpecies() == fluxCarrier) {
                rf = rf_tmp;
            }
        }
        // 
        // if "inside" speciesContext is not "fixed", add flux to ResolvedFlux
        // 
        SpeciesContext insideSpeciesContext = mathMapping_4_8.getSimulationContext().getModel().getSpeciesContext(fluxCarrier, structTopology.getInsideFeature(getMembrane()));
        SpeciesContextSpec insideSpeciesContextSpec = mathMapping_4_8.getSimulationContext().getReactionContext().getSpeciesContextSpec(insideSpeciesContext);
        // if (!insideSpeciesContextSpec.isConstant()){
        if (bNoFluxIfFixed || !insideSpeciesContextSpec.isConstant()) {
            if (bNoFluxIfFixed && insideSpeciesContextSpec.isConstant()) {
                bNoFluxIfFixedExercised = true;
            }
            if (rf == null) {
                rf = new ResolvedFlux(fluxCarrier, fr.getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate).getUnitDefinition());
                resolvedFluxList.addElement(rf);
            }
            FeatureMapping insideFeatureMapping = (FeatureMapping) geoContext.getStructureMapping((structTopology.getInsideFeature((Membrane) fr.getStructure())));
            Expression residualVolumeFraction = mathMapping_4_8.getResidualVolumeFraction(insideFeatureMapping).renameBoundSymbols(mathMapping_4_8.getNameScope());
            Expression insideFluxCorrection = Expression.invert(residualVolumeFraction);
            // 
            if (bResolvedFluxCorrectionBug && !residualVolumeFraction.compareEqual(new Expression(1.0))) {
                bResolvedFluxCorrectionBugExercised = true;
                System.out.println("MembraneStructureAnalyzer.refreshResolvedFluxes() ... 'ResolvedFluxCorrection' bug compatability mode");
                insideFluxCorrection = new Expression(1.0);
            }
            // 
            if (fr.getKinetics() instanceof DistributedKinetics) {
                Expression reactionRateParameter = new Expression(((DistributedKinetics) fr.getKinetics()).getReactionRateParameter(), mathMapping_4_8.getNameScope());
                if (rf.inFluxExpression.isZero()) {
                    rf.inFluxExpression = Expression.mult(reactionRateParameter, insideFluxCorrection).flatten();
                } else {
                    rf.inFluxExpression = Expression.add(rf.inFluxExpression, Expression.mult(reactionRateParameter, insideFluxCorrection).flatten());
                }
            } else if (fr.getKinetics() instanceof LumpedKinetics) {
                throw new RuntimeException("Lumped Kinetics for fluxes not yet supported");
            } else {
                throw new RuntimeException("unexpected Kinetic type in MembraneStructureAnalyzer.refreshResolvedFluxes()");
            }
        // rf.inFlux.bindExpression(mathMapping);
        }
        SpeciesContext outsideSpeciesContext = mathMapping_4_8.getSimulationContext().getModel().getSpeciesContext(fluxCarrier, structTopology.getOutsideFeature(getMembrane()));
        SpeciesContextSpec outsideSpeciesContextSpec = mathMapping_4_8.getSimulationContext().getReactionContext().getSpeciesContextSpec(outsideSpeciesContext);
        // if (!outsideSpeciesContextSpec.isConstant()){
        if (bNoFluxIfFixed || !outsideSpeciesContextSpec.isConstant()) {
            if (bNoFluxIfFixed && outsideSpeciesContextSpec.isConstant()) {
                bNoFluxIfFixedExercised = true;
            }
            if (rf == null) {
                rf = new ResolvedFlux(fluxCarrier, fr.getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate).getUnitDefinition());
                resolvedFluxList.addElement(rf);
            }
            FeatureMapping outsideFeatureMapping = (FeatureMapping) geoContext.getStructureMapping(structTopology.getOutsideFeature((Membrane) fr.getStructure()));
            Expression residualVolumeFraction = mathMapping_4_8.getResidualVolumeFraction(outsideFeatureMapping).renameBoundSymbols(mathMapping_4_8.getNameScope());
            Expression outsideFluxCorrection = Expression.invert(residualVolumeFraction);
            // 
            if (bResolvedFluxCorrectionBug && !residualVolumeFraction.compareEqual(new Expression(1.0))) {
                bResolvedFluxCorrectionBugExercised = true;
                System.out.println("MembraneStructureAnalyzer.refreshResolvedFluxes() ... 'ResolvedFluxCorrection' bug compatability mode");
                outsideFluxCorrection = new Expression(1.0);
            }
            // 
            if (fr.getKinetics() instanceof DistributedKinetics) {
                Expression reactionRateParameter = new Expression(((DistributedKinetics) fr.getKinetics()).getReactionRateParameter(), mathMapping_4_8.getNameScope());
                if (rf.outFluxExpression.isZero()) {
                    rf.outFluxExpression = Expression.mult(Expression.negate(reactionRateParameter), outsideFluxCorrection).flatten();
                } else {
                    rf.outFluxExpression = Expression.add(rf.outFluxExpression, Expression.mult(Expression.negate(reactionRateParameter), outsideFluxCorrection).flatten());
                }
            } else if (fr.getKinetics() instanceof LumpedKinetics) {
                throw new RuntimeException("Lumped Kinetics not yet supported for Flux Reaction: " + fr.getName());
            } else {
                throw new RuntimeException("unexpected Kinetics type for Flux Reaction " + fr.getName());
            }
        // rf.outFlux.bindExpression(mathMapping);
        }
    }
    // 
    // for each reaction, incorporate all reactionSteps involving binding with volumetric species
    // 
    double kMoleValue = Model.reservedConstantsMap.get(ReservedSymbolRole.KMOLE);
    for (int i = 0; i < reactionSpecs.length; i++) {
        if (reactionSpecs[i].isExcluded()) {
            continue;
        }
        ReactionStep rs = reactionSpecs[i].getReactionStep();
        if (rs.getStructure() == getMembrane()) {
            if (rs instanceof SimpleReaction) {
                SimpleReaction sr = (SimpleReaction) rs;
                ReactionParticipant[] rp_Array = sr.getReactionParticipants();
                for (int k = 0; k < rp_Array.length; k++) {
                    if (rp_Array[k] instanceof Reactant || rp_Array[k] instanceof Product) {
                        SpeciesContextSpec scs = mathMapping_4_8.getSimulationContext().getReactionContext().getSpeciesContextSpec(rp_Array[k].getSpeciesContext());
                        // if (rp_Array[k].getStructure() instanceof Feature && !scs.isConstant()){
                        if (rp_Array[k].getStructure() instanceof Feature && (bNoFluxIfFixed || !scs.isConstant())) {
                            if (bNoFluxIfFixed && scs.isConstant()) {
                                bNoFluxIfFixedExercised = true;
                            }
                            // 
                            // for each Reactant or Product binding to this membrane...
                            // 
                            // 
                            // get ResolvedFlux for this species
                            // 
                            ResolvedFlux rf = null;
                            for (int j = 0; j < resolvedFluxList.size(); j++) {
                                ResolvedFlux rf_tmp = (ResolvedFlux) resolvedFluxList.elementAt(j);
                                if (rf_tmp.getSpecies() == rp_Array[k].getSpecies()) {
                                    rf = rf_tmp;
                                }
                            }
                            if (rf == null) {
                                rf = new ResolvedFlux(rp_Array[k].getSpecies(), sr.getKinetics().getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate).getUnitDefinition());
                                resolvedFluxList.addElement(rf);
                            }
                            Expression reactionRateExpression = getReactionRateExpression(sr, rp_Array[k]).renameBoundSymbols(mathMapping_4_8.getNameScope());
                            if (rp_Array[k].getStructure() == structTopology.getInsideFeature(getMembrane())) {
                                // 
                                // for binding on inside, add to ResolvedFlux.inFlux
                                // 
                                FeatureMapping insideFeatureMapping = (FeatureMapping) geoContext.getStructureMapping(structTopology.getInsideFeature(getMembrane()));
                                Expression residualVolumeFraction = mathMapping_4_8.getResidualVolumeFraction(insideFeatureMapping).renameBoundSymbols(mathMapping_4_8.getNameScope());
                                Expression insideFluxCorrection = Expression.div(new Expression(kMoleValue), residualVolumeFraction).flatten();
                                // 
                                if (bResolvedFluxCorrectionBug && !residualVolumeFraction.compareEqual(new Expression(1.0))) {
                                    bResolvedFluxCorrectionBugExercised = true;
                                    System.out.println("MembraneStructureAnalyzer.refreshResolvedFluxes() ... 'ResolvedFluxCorrection' bug compatability mode");
                                    insideFluxCorrection = new Expression(kMoleValue);
                                }
                                if (rf.inFluxExpression.isZero()) {
                                    rf.inFluxExpression = Expression.mult(insideFluxCorrection, reactionRateExpression);
                                } else {
                                    rf.inFluxExpression = Expression.add(rf.inFluxExpression, Expression.mult(insideFluxCorrection, reactionRateExpression));
                                }
                            // rf.inFlux.bindExpression(mathMapping);
                            } else if (rp_Array[k].getStructure() == structTopology.getOutsideFeature(getMembrane())) {
                                // 
                                // for binding on outside, add to ResolvedFlux.outFlux
                                // 
                                FeatureMapping outsideFeatureMapping = (FeatureMapping) geoContext.getStructureMapping(structTopology.getOutsideFeature(getMembrane()));
                                Expression residualVolumeFraction = mathMapping_4_8.getResidualVolumeFraction(outsideFeatureMapping).renameBoundSymbols(mathMapping_4_8.getNameScope());
                                Expression outsideFluxCorrection = Expression.div(new Expression(kMoleValue), residualVolumeFraction).flatten();
                                // 
                                if (bResolvedFluxCorrectionBug && !residualVolumeFraction.compareEqual(new Expression(1.0))) {
                                    bResolvedFluxCorrectionBugExercised = true;
                                    System.out.println("MembraneStructureAnalyzer.refreshResolvedFluxes() ... 'ResolvedFluxCorrection' bug compatability mode");
                                    outsideFluxCorrection = new Expression(kMoleValue);
                                }
                                if (rf.outFluxExpression.isZero()) {
                                    rf.outFluxExpression = Expression.mult(outsideFluxCorrection, reactionRateExpression);
                                } else {
                                    rf.outFluxExpression = Expression.add(rf.outFluxExpression, Expression.mult(outsideFluxCorrection, reactionRateExpression));
                                }
                            // rf.outFlux.bindExpression(mathMapping);
                            } else {
                                throw new Exception("SpeciesContext " + rp_Array[k].getSpeciesContext().getName() + " doesn't border membrane " + getMembrane().getName() + " but reacts there");
                            }
                        }
                    }
                }
            }
        }
    }
    // 
    if (resolvedFluxList.size() > 0) {
        resolvedFluxes = new ResolvedFlux[resolvedFluxList.size()];
        resolvedFluxList.copyInto(resolvedFluxes);
    } else {
        resolvedFluxes = null;
    }
}

Example 79 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();
    Model model = mathMapping_4_8.getSimulationContext().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_4_8.getSpeciesContextMapping(sc);
            SpeciesContextSpec scs = mathMapping_4_8.getSimulationContext().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_4_8.getSpeciesContextMapping(sc);
            SpeciesContextSpec scs = mathMapping_4_8.getSimulationContext().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 = mathMapping_4_8.getSimulationContext().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 (scm.isPDERequired()) {
            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()));
        }
        SimulationContext simContext = mathMapping_4_8.getSimulationContext();
        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 80 with ReactionStep

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

the class StructureAnalyzer method refreshFastSpeciesContextMappings.

/**
 * This method was created in VisualAge.
 */
private void refreshFastSpeciesContextMappings() {
    // System.out.println("StructureAnalyzer.refreshFastSpeciesContextMappings()");
    // GeometryContext geoContext = mathMapping.getSimulationContext().getGeometryContext();
    Vector<SpeciesContextMapping> scFastList = new Vector<SpeciesContextMapping>();
    // 
    if (structures == null) {
        return;
    }
    for (int i = 0; i < structures.length; i++) {
        SpeciesContext[] speciesContexts = mathMapping_4_8.getSimulationContext().getReactionContext().getModel().getSpeciesContexts(structures[i]);
        for (int j = 0; j < speciesContexts.length; j++) {
            SpeciesContext sc = speciesContexts[j];
            SpeciesContextMapping scm = mathMapping_4_8.getSpeciesContextMapping(sc);
            SpeciesContextSpec scs = mathMapping_4_8.getSimulationContext().getReactionContext().getSpeciesContextSpec(sc);
            if (scm.getDependencyExpression() == null && scs.isConstant() == false) {
                // right now all independent vars
                scFastList.addElement(scm);
            }
        }
    }
    if (scFastList.size() > 0) {
        fastSpeciesContextMappings = new SpeciesContextMapping[scFastList.size()];
        scFastList.copyInto(fastSpeciesContextMappings);
    // for (int i=0;i<fastSpeciesContextMappings.length;i++){
    // System.out.println("fastSpeciesContextMappings["+i+"] = "+fastSpeciesContextMappings[i].getSpeciesContext().getName());
    // }
    } else {
        fastSpeciesContextMappings = null;
    }
    // System.out.println("StructureAnalyzer.refreshFastSpeciesContextMapping(), fastSpeciesContextMappings.length = "+scFastList.size());
    // 
    // for each reaction, get all reactionSteps associated with these structures
    // 
    Vector<ReactionStep> rsFastList = new Vector<ReactionStep>();
    ReactionSpec[] reactionSpecs = mathMapping_4_8.getSimulationContext().getReactionContext().getReactionSpecs();
    for (int i = 0; i < reactionSpecs.length; i++) {
        ReactionStep rs = reactionSpecs[i].getReactionStep();
        if (reactionSpecs[i].isExcluded()) {
            continue;
        }
        for (int j = 0; j < structures.length; j++) {
            if (rs.getStructure() == structures[j]) {
                if (reactionSpecs[i].isFast()) {
                    rsFastList.addElement(rs);
                }
            }
        }
    }
    if (rsFastList.size() > 0) {
        fastReactionSteps = new ReactionStep[rsFastList.size()];
        rsFastList.copyInto(fastReactionSteps);
    } else {
        fastReactionSteps = null;
        fastSpeciesContextMappings = null;
    }
// System.out.println("StructureAnalyzer.refreshFastSpeciesContextMapping(), reactionSteps.length = "+scFastList.size());
}