Search in sources :

Example 21 with FluxReaction

use of cbit.vcell.model.FluxReaction 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 = 1 / 602.0;
    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;
    }
}
Also used : LumpedKinetics(cbit.vcell.model.LumpedKinetics) Product(cbit.vcell.model.Product) FluxReaction(cbit.vcell.model.FluxReaction) SpeciesContext(cbit.vcell.model.SpeciesContext) SpeciesContextSpec(cbit.vcell.mapping.SpeciesContextSpec) Reactant(cbit.vcell.model.Reactant) Feature(cbit.vcell.model.Feature) FeatureMapping(cbit.vcell.mapping.FeatureMapping) Membrane(cbit.vcell.model.Membrane) GeometryContext(cbit.vcell.mapping.GeometryContext) Vector(java.util.Vector) Species(cbit.vcell.model.Species) DistributedKinetics(cbit.vcell.model.DistributedKinetics) SimpleReaction(cbit.vcell.model.SimpleReaction) StructureTopology(cbit.vcell.model.Model.StructureTopology) ReactionSpec(cbit.vcell.mapping.ReactionSpec) Expression(cbit.vcell.parser.Expression) ReactionStep(cbit.vcell.model.ReactionStep) ReactionParticipant(cbit.vcell.model.ReactionParticipant)

Example 22 with FluxReaction

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

the class StochMathMapping_4_8 method refreshMathDescription.

/**
 * set up a math description based on current simulationContext.
 */
private void refreshMathDescription() throws MappingException, MatrixException, MathException, ExpressionException, ModelException {
    // use local variable instead of using getter all the time.
    SimulationContext simContext = getSimulationContext();
    // local structure mapping list
    StructureMapping[] structureMappings = simContext.getGeometryContext().getStructureMappings();
    // We have to check if all the reactions are able to tranform to stochastic jump processes before generating the math.
    String stochChkMsg = simContext.getModel().isValidForStochApp();
    if (!(stochChkMsg.equals(""))) {
        throw new ModelException("Problem updating math description: " + simContext.getName() + "\n" + stochChkMsg);
    }
    // All sizes must be set for new ODE models and ratios must be set for old ones.
    simContext.checkValidity();
    // 
    // verify that all structures are mapped to subvolumes and all subvolumes are mapped to a structure
    // 
    Structure[] structures = simContext.getGeometryContext().getModel().getStructures();
    for (int i = 0; i < structures.length; i++) {
        StructureMapping sm = simContext.getGeometryContext().getStructureMapping(structures[i]);
        if (sm == null || (sm instanceof FeatureMapping && getSubVolume(((FeatureMapping) sm)) == null)) {
            throw new MappingException("model structure '" + structures[i].getName() + "' not mapped to a geometry subVolume");
        }
        if (sm != null && (sm instanceof MembraneMapping) && ((MembraneMapping) sm).getVolumeFractionParameter() != null) {
            Expression volFractExp = ((MembraneMapping) sm).getVolumeFractionParameter().getExpression();
            try {
                if (volFractExp != null) {
                    double volFract = volFractExp.evaluateConstant();
                    if (volFract >= 1.0) {
                        throw new MappingException("model structure '" + (getSimulationContext().getModel().getStructureTopology().getInsideFeature(((MembraneMapping) sm).getMembrane()).getName() + "' has volume fraction >= 1.0"));
                    }
                }
            } catch (ExpressionException e) {
                e.printStackTrace(System.out);
            }
        }
    }
    SubVolume[] subVolumes = simContext.getGeometryContext().getGeometry().getGeometrySpec().getSubVolumes();
    for (int i = 0; i < subVolumes.length; i++) {
        if (getStructures(subVolumes[i]) == null || getStructures(subVolumes[i]).length == 0) {
            throw new MappingException("geometry subVolume '" + subVolumes[i].getName() + "' not mapped from a model structure");
        }
    }
    // 
    // gather only those reactionSteps that are not "excluded"
    // 
    ReactionSpec[] reactionSpecs = simContext.getReactionContext().getReactionSpecs();
    Vector<ReactionStep> rsList = new Vector<ReactionStep>();
    for (int i = 0; i < reactionSpecs.length; i++) {
        if (reactionSpecs[i].isExcluded() == false) {
            rsList.add(reactionSpecs[i].getReactionStep());
        }
    }
    ReactionStep[] reactionSteps = new ReactionStep[rsList.size()];
    rsList.copyInto(reactionSteps);
    // 
    for (int i = 0; i < reactionSteps.length; i++) {
        Kinetics.UnresolvedParameter[] unresolvedParameters = reactionSteps[i].getKinetics().getUnresolvedParameters();
        if (unresolvedParameters != null && unresolvedParameters.length > 0) {
            StringBuffer buffer = new StringBuffer();
            for (int j = 0; j < unresolvedParameters.length; j++) {
                if (j > 0) {
                    buffer.append(", ");
                }
                buffer.append(unresolvedParameters[j].getName());
            }
            throw new MappingException(reactionSteps[i].getDisplayType() + " '" + reactionSteps[i].getName() + "' contains unresolved identifier(s): " + buffer);
        }
    }
    // 
    // create new MathDescription (based on simContext's previous MathDescription if possible)
    // 
    MathDescription oldMathDesc = simContext.getMathDescription();
    mathDesc = null;
    if (oldMathDesc != null) {
        if (oldMathDesc.getVersion() != null) {
            mathDesc = new MathDescription(oldMathDesc.getVersion());
        } else {
            mathDesc = new MathDescription(oldMathDesc.getName());
        }
    } else {
        mathDesc = new MathDescription(simContext.getName() + "_generated");
    }
    // 
    // temporarily place all variables in a hashtable (before binding) and discarding duplicates
    // 
    VariableHash varHash = new VariableHash();
    // 
    // conversion factors
    // 
    Model model = simContext.getModel();
    ModelUnitSystem modelUnitSystem = model.getUnitSystem();
    varHash.addVariable(new Constant(getMathSymbol(model.getKMOLE(), null), getIdentifierSubstitutions(model.getKMOLE().getExpression(), model.getKMOLE().getUnitDefinition(), null)));
    varHash.addVariable(new Constant(getMathSymbol(model.getN_PMOLE(), null), getIdentifierSubstitutions(model.getN_PMOLE().getExpression(), model.getN_PMOLE().getUnitDefinition(), null)));
    varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT().getExpression(), model.getFARADAY_CONSTANT().getUnitDefinition(), null)));
    varHash.addVariable(new Constant(getMathSymbol(model.getFARADAY_CONSTANT_NMOLE(), null), getIdentifierSubstitutions(model.getFARADAY_CONSTANT_NMOLE().getExpression(), model.getFARADAY_CONSTANT_NMOLE().getUnitDefinition(), null)));
    varHash.addVariable(new Constant(getMathSymbol(model.getGAS_CONSTANT(), null), getIdentifierSubstitutions(model.getGAS_CONSTANT().getExpression(), model.getGAS_CONSTANT().getUnitDefinition(), null)));
    varHash.addVariable(new Constant(getMathSymbol(model.getTEMPERATURE(), null), getIdentifierSubstitutions(new Expression(simContext.getTemperatureKelvin()), model.getTEMPERATURE().getUnitDefinition(), null)));
    Enumeration<SpeciesContextMapping> enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = enum1.nextElement();
        if (scm.getVariable() instanceof StochVolVariable) {
            varHash.addVariable(scm.getVariable());
        }
    }
    // 
    // add rate term for all reactions
    // add current source terms for each reaction step in a membrane
    // 
    /*for (int i = 0; i < reactionSteps.length; i++){
			boolean bAllReactionParticipantsFixed = true;
			ReactionParticipant rp_Array[] = reactionSteps[i].getReactionParticipants();
			for (int j = 0; j < rp_Array.length; j++) {
				SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(rp_Array[j].getSpeciesContext());
				if (!(rp_Array[j] instanceof Catalyst) && !scs.isConstant()){
					bAllReactionParticipantsFixed = false;  // found at least one reactionParticipant that is not fixed and needs this rate
				}
			}
			StructureMapping sm = simContext.getGeometryContext().getStructureMapping(reactionSteps[i].getStructure());
		}---don't think it's useful, isn't it?*/
    // deals with model parameters
    ModelParameter[] modelParameters = simContext.getModel().getModelParameters();
    for (int j = 0; j < modelParameters.length; j++) {
        Expression expr = getSubstitutedExpr(modelParameters[j].getExpression(), true, false);
        expr = getIdentifierSubstitutions(expr, modelParameters[j].getUnitDefinition(), null);
        varHash.addVariable(newFunctionOrConstant(getMathSymbol(modelParameters[j], null), expr));
    }
    // added July 2009, ElectricalStimulusParameter electric mapping tab
    ElectricalStimulus[] elecStimulus = simContext.getElectricalStimuli();
    if (elecStimulus.length > 0) {
        throw new MappingException("Modles with electrophysiology are not supported for stochastic applications.");
    }
    for (int j = 0; j < structureMappings.length; j++) {
        if (structureMappings[j] instanceof MembraneMapping) {
            MembraneMapping memMapping = (MembraneMapping) structureMappings[j];
            Parameter initialVoltageParm = memMapping.getInitialVoltageParameter();
            try {
                Expression exp = initialVoltageParm.getExpression();
                exp.evaluateConstant();
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(memMapping.getMembrane().getMembraneVoltage(), memMapping), getIdentifierSubstitutions(memMapping.getInitialVoltageParameter().getExpression(), memMapping.getInitialVoltageParameter().getUnitDefinition(), memMapping)));
            } catch (ExpressionException e) {
                e.printStackTrace(System.out);
                throw new MappingException("Membrane initial voltage: " + initialVoltageParm.getName() + " cannot be evaluated as constant.");
            }
        }
    }
    // 
    for (int j = 0; j < reactionSteps.length; j++) {
        ReactionStep rs = reactionSteps[j];
        if (simContext.getReactionContext().getReactionSpec(rs).isExcluded()) {
            continue;
        }
        if (rs.getKinetics() instanceof LumpedKinetics) {
            throw new RuntimeException("Lumped Kinetics not yet supported for Stochastic Math Generation");
        }
        Kinetics.KineticsParameter[] parameters = rs.getKinetics().getKineticsParameters();
        StructureMapping sm = simContext.getGeometryContext().getStructureMapping(rs.getStructure());
        if (parameters != null) {
            for (int i = 0; i < parameters.length; i++) {
                if ((parameters[i].getRole() == Kinetics.ROLE_CurrentDensity) && (parameters[i].getExpression() == null || parameters[i].getExpression().isZero())) {
                    continue;
                }
                // don't add rate, we'll do it later when creating the jump processes
                if (parameters[i].getRole() != Kinetics.ROLE_ReactionRate) {
                    Expression expr = getSubstitutedExpr(parameters[i].getExpression(), true, false);
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(parameters[i], sm), getIdentifierSubstitutions(expr, parameters[i].getUnitDefinition(), sm)));
                }
            }
        }
    }
    // the parameter "Size" is already put into mathsymbolmapping in refreshSpeciesContextMapping()
    for (int i = 0; i < structureMappings.length; i++) {
        StructureMapping sm = structureMappings[i];
        StructureMapping.StructureMappingParameter parm = sm.getParameterFromRole(StructureMapping.ROLE_Size);
        if (parm.getExpression() != null) {
            try {
                double value = parm.getExpression().evaluateConstant();
                varHash.addVariable(new Constant(getMathSymbol(parm, sm), new Expression(value)));
            } catch (ExpressionException e) {
                // varHash.addVariable(new Function(getMathSymbol0(parm,sm),getIdentifierSubstitutions(parm.getExpression(),parm.getUnitDefinition(),sm)));
                e.printStackTrace(System.out);
                throw new MappingException("Size of structure:" + sm.getNameScope().getName() + " cannot be evaluated as constant.");
            }
        }
    }
    // 
    // species initial values (either function or constant)
    // 
    SpeciesContextSpec[] speciesContextSpecs = simContext.getReactionContext().getSpeciesContextSpecs();
    for (int i = 0; i < speciesContextSpecs.length; i++) {
        // can be concentration or amount
        SpeciesContextSpec.SpeciesContextSpecParameter initParam = null;
        Expression iniExp = null;
        StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
        if (speciesContextSpecs[i].getInitialConcentrationParameter() != null && speciesContextSpecs[i].getInitialConcentrationParameter().getExpression() != null) {
            // use concentration, need to set up amount functions
            initParam = speciesContextSpecs[i].getInitialConcentrationParameter();
            iniExp = initParam.getExpression();
            iniExp = getSubstitutedExpr(iniExp, true, !speciesContextSpecs[i].isConstant());
            // now create the appropriate function or Constant for the speciesContextSpec.
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(initParam, sm), getIdentifierSubstitutions(iniExp, initParam.getUnitDefinition(), sm)));
            // add function for initial amount
            SpeciesContextSpec.SpeciesContextSpecParameter initAmountParam = speciesContextSpecs[i].getInitialCountParameter();
            Expression iniAmountExp = getExpressionConcToAmt(new Expression(initParam, getNameScope()), speciesContextSpecs[i].getSpeciesContext());
            // iniAmountExp.bindExpression(this);
            varHash.addVariable(new Function(getMathSymbol(initAmountParam, sm), getIdentifierSubstitutions(iniAmountExp, initAmountParam.getUnitDefinition(), sm), nullDomain));
        } else if (speciesContextSpecs[i].getInitialCountParameter() != null && speciesContextSpecs[i].getInitialCountParameter().getExpression() != null) {
            // use amount
            initParam = speciesContextSpecs[i].getInitialCountParameter();
            iniExp = initParam.getExpression();
            iniExp = getSubstitutedExpr(iniExp, false, !speciesContextSpecs[i].isConstant());
            // now create the appropriate function or Constant for the speciesContextSpec.
            varHash.addVariable(newFunctionOrConstant(getMathSymbol(initParam, sm), getIdentifierSubstitutions(iniExp, initParam.getUnitDefinition(), sm)));
        }
        // add spConcentration (concentration of species) to varHash as function or constant
        SpeciesConcentrationParameter spConcParam = getSpeciesConcentrationParameter(speciesContextSpecs[i].getSpeciesContext());
        varHash.addVariable(newFunctionOrConstant(getMathSymbol(spConcParam, sm), getIdentifierSubstitutions(spConcParam.getExpression(), spConcParam.getUnitDefinition(), sm)));
    }
    // 
    // constant species (either function or constant)
    // 
    enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
        if (scm.getVariable() instanceof Constant) {
            varHash.addVariable(scm.getVariable());
        }
    }
    // 
    if (simContext.getGeometryContext().getGeometry() != null) {
        try {
            mathDesc.setGeometry(simContext.getGeometryContext().getGeometry());
        } catch (java.beans.PropertyVetoException e) {
            e.printStackTrace(System.out);
            throw new MappingException("failure setting geometry " + e.getMessage());
        }
    } else {
        throw new MappingException("geometry must be defined");
    }
    // 
    // functions: species which is not a variable, but has dependency expression
    // 
    enum1 = getSpeciesContextMappings();
    while (enum1.hasMoreElements()) {
        SpeciesContextMapping scm = (SpeciesContextMapping) enum1.nextElement();
        if (scm.getVariable() == null && scm.getDependencyExpression() != null) {
            StructureMapping sm = simContext.getGeometryContext().getStructureMapping(scm.getSpeciesContext().getStructure());
            Expression exp = scm.getDependencyExpression();
            exp.bindExpression(this);
            SpeciesCountParameter spCountParam = getSpeciesCountParameter(scm.getSpeciesContext());
            varHash.addVariable(new Function(getMathSymbol(spCountParam, sm), getIdentifierSubstitutions(exp, spCountParam.getUnitDefinition(), sm), nullDomain));
        }
    }
    // 
    // create subDomains
    // 
    SubDomain subDomain = null;
    subVolumes = simContext.getGeometryContext().getGeometry().getGeometrySpec().getSubVolumes();
    for (int j = 0; j < subVolumes.length; j++) {
        SubVolume subVolume = (SubVolume) subVolumes[j];
        // 
        // get priority of subDomain
        // 
        int priority;
        Feature spatialFeature = getResolvedFeature(subVolume);
        if (spatialFeature == null) {
            if (simContext.getGeometryContext().getGeometry().getDimension() > 0) {
                throw new MappingException("no compartment (in Physiology) is mapped to subdomain '" + subVolume.getName() + "' (in Geometry)");
            } else {
                priority = CompartmentSubDomain.NON_SPATIAL_PRIORITY;
            }
        } else {
            // now does not have to match spatial feature, *BUT* needs to be unique
            priority = j;
        }
        subDomain = new CompartmentSubDomain(subVolume.getName(), priority);
        mathDesc.addSubDomain(subDomain);
    }
    // ReactionSpec[] reactionSpecs = simContext.getReactionContext().getReactionSpecs();---need to take a look here!
    for (int i = 0; i < reactionSpecs.length; i++) {
        if (reactionSpecs[i].isExcluded()) {
            continue;
        }
        // get the reaction
        ReactionStep reactionStep = reactionSpecs[i].getReactionStep();
        Kinetics kinetics = reactionStep.getKinetics();
        // the structure where reaction happens
        StructureMapping sm = simContext.getGeometryContext().getStructureMapping(reactionStep.getStructure());
        // create symbol table for jump process based on reactionStep and structure mapping
        // final ReactionStep finalRS = reactionStep;
        // final StructureMapping finalSM = sm;
        // SymbolTable symTable = new SymbolTable(){
        // public SymbolTableEntry getEntry(String identifierString) throws ExpressionBindingException {
        // SymbolTableEntry ste = finalRS.getEntry(identifierString);
        // if(ste == null)
        // {
        // ste = finalSM.getEntry(identifierString);
        // }
        // return ste;
        // }
        // };
        // Different ways to deal with simple reactions and flux reactions
        // probability parameter from modelUnitSystem
        VCUnitDefinition probabilityParamUnit = modelUnitSystem.getStochasticSubstanceUnit().divideBy(modelUnitSystem.getTimeUnit());
        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)) {
                forwardRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KForward).getExpression();
                reverseRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KReverse).getExpression();
            } else if (kinetics.getKineticsDescription().equals(KineticsDescription.General)) {
                Expression rateExp = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate).getExpression();
                MassActionSolver.MassActionFunction maFunc = MassActionSolver.solveMassAction(null, null, 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 (kinetics.getKineticsDescription().getName().compareTo(KineticsDescription.HMM_irreversible.getName())==0)
			    {
				    forwardRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_Km).getExpression();
				}
			    else if (kinetics.getKineticsDescription().getName().compareTo(KineticsDescription.HMM_reversible.getName())==0)
			    {
					forwardRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KmFwd).getExpression();
					reverseRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_KmRev).getExpression();
				}*/
            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 mass actions
                exp = getProbabilityRate(reactionStep, true);
                // bind symbol table before substitute identifiers in the reaction step
                exp.bindExpression(this);
                MathMapping_4_8.ProbabilityParameter probParm = null;
                try {
                    probParm = addProbabilityParameter("P_" + jpName, exp, MathMapping_4_8.PARAMETER_ROLE_P, probabilityParamUnit, reactionSpecs[i]);
                } catch (PropertyVetoException pve) {
                    pve.printStackTrace();
                    throw new MappingException(pve.getMessage());
                }
                // add probability to function or constant
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(probParm, sm), getIdentifierSubstitutions(exp, probabilityParamUnit, sm)));
                JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probParm, sm)));
                // 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 = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression("-" + String.valueOf(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 = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", 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()) + "_reverse";
                Expression exp = null;
                // reactions are mass actions
                exp = getProbabilityRate(reactionStep, false);
                // bind symbol table before substitute identifiers in the reaction step
                exp.bindExpression(this);
                MathMapping_4_8.ProbabilityParameter probRevParm = null;
                try {
                    probRevParm = addProbabilityParameter("P_" + jpName, exp, MathMapping_4_8.PARAMETER_ROLE_P_reverse, probabilityParamUnit, reactionSpecs[i]);
                } catch (PropertyVetoException pve) {
                    pve.printStackTrace();
                    throw new MappingException(pve.getMessage());
                }
                // add probability to function or constant
                varHash.addVariable(newFunctionOrConstant(getMathSymbol(probRevParm, sm), getIdentifierSubstitutions(exp, probabilityParamUnit, sm)));
                JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probRevParm, sm)));
                // 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 = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", 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 = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression("-" + String.valueOf(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)) {
                Expression fluxRate = kinetics.getKineticsParameterFromRole(Kinetics.ROLE_ReactionRate).getExpression();
                // we have to pass the math description para to flux solver, coz somehow math description in simulation context is not updated.
                MassActionSolver.MassActionFunction fluxFunc = MassActionSolver.solveMassAction(null, null, fluxRate, (FluxReaction) reactionStep);
                // create jump process for forward flux if it exists.
                if (fluxFunc.getForwardRate() != null && !fluxFunc.getForwardRate().isZero()) {
                    // jump process name
                    // +"_reverse";
                    String jpName = TokenMangler.mangleToSName(reactionStep.getName());
                    // we do it here instead of fluxsolver, coz we need to use getMathSymbol0(), structuremapping...etc.
                    Expression rate = fluxFunc.getForwardRate();
                    // get species expression (depend on structure, if mem: Species/mem_Size, if vol: species*KMOLE/vol_size)
                    SpeciesContext scOut = fluxFunc.getReactants().get(0).getSpeciesContext();
                    Expression speciesFactor = null;
                    if (scOut.getStructure() instanceof Feature) {
                        Expression exp1 = new Expression(1.0 / 602.0);
                        Expression exp2 = new Expression(scOut.getStructure().getStructureSize(), getNameScope());
                        speciesFactor = Expression.div(Expression.invert(exp1), exp2);
                    } else {
                        throw new MappingException("Species involved in a flux have to be volume species.");
                    }
                    Expression speciesExp = Expression.mult(speciesFactor, new Expression(scOut, getNameScope()));
                    // get probability expression by adding factor to rate (rate: rate*size_mem/KMOLE)
                    Expression expr1 = Expression.mult(rate, speciesExp);
                    Expression numeratorExpr = Expression.mult(expr1, new Expression(sm.getStructure().getStructureSize(), getNameScope()));
                    Expression exp = new Expression(1.0 / 602.0);
                    Expression probExp = Expression.mult(numeratorExpr, exp);
                    // bind symbol table before substitute identifiers in the reaction step
                    probExp.bindExpression(reactionStep);
                    MathMapping_4_8.ProbabilityParameter probParm = null;
                    try {
                        probParm = addProbabilityParameter("P_" + jpName, probExp, MathMapping_4_8.PARAMETER_ROLE_P, probabilityParamUnit, reactionSpecs[i]);
                    } catch (PropertyVetoException pve) {
                        pve.printStackTrace();
                        throw new MappingException(pve.getMessage());
                    }
                    // add probability to function or constant
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(probParm, sm), getIdentifierSubstitutions(probExp, probabilityParamUnit, sm)));
                    JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probParm, sm)));
                    // actions
                    Action action = null;
                    SpeciesContext sc = fluxFunc.getReactants().get(0).getSpeciesContext();
                    if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
                        action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression(-1));
                        jp.addAction(action);
                    }
                    sc = fluxFunc.getProducts().get(0).getSpeciesContext();
                    if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
                        action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression(1));
                        jp.addAction(action);
                    }
                    subDomain.addJumpProcess(jp);
                }
                if (fluxFunc.getReverseRate() != null && !fluxFunc.getReverseRate().isZero()) {
                    // jump process name
                    String jpName = TokenMangler.mangleToSName(reactionStep.getName()) + "_reverse";
                    Expression rate = fluxFunc.getReverseRate();
                    // get species expression (depend on structure, if mem: Species/mem_Size, if vol: species*KMOLE/vol_size)
                    SpeciesContext scIn = fluxFunc.getProducts().get(0).getSpeciesContext();
                    Expression speciesFactor = null;
                    if (scIn.getStructure() instanceof Feature) {
                        Expression exp1 = new Expression(1.0 / 602.0);
                        Expression exp2 = new Expression(scIn.getStructure().getStructureSize(), getNameScope());
                        speciesFactor = Expression.div(Expression.invert(exp1), exp2);
                    } else {
                        throw new MappingException("Species involved in a flux have to be volume species.");
                    }
                    Expression speciesExp = Expression.mult(speciesFactor, new Expression(scIn, getNameScope()));
                    // get probability expression by adding factor to rate (rate: rate*size_mem/KMOLE)
                    Expression expr1 = Expression.mult(rate, speciesExp);
                    Expression numeratorExpr = Expression.mult(expr1, new Expression(sm.getStructure().getStructureSize(), getNameScope()));
                    Expression exp = new Expression(1.0 / 602.0);
                    Expression probRevExp = Expression.mult(numeratorExpr, exp);
                    // bind symbol table before substitute identifiers in the reaction step
                    probRevExp.bindExpression(reactionStep);
                    MathMapping_4_8.ProbabilityParameter probRevParm = null;
                    try {
                        probRevParm = addProbabilityParameter("P_" + jpName, probRevExp, MathMapping_4_8.PARAMETER_ROLE_P_reverse, probabilityParamUnit, reactionSpecs[i]);
                    } catch (PropertyVetoException pve) {
                        pve.printStackTrace();
                        throw new MappingException(pve.getMessage());
                    }
                    // add probability to function or constant
                    varHash.addVariable(newFunctionOrConstant(getMathSymbol(probRevParm, sm), getIdentifierSubstitutions(probRevExp, probabilityParamUnit, sm)));
                    JumpProcess jp = new JumpProcess(jpName, new Expression(getMathSymbol(probRevParm, sm)));
                    // actions
                    Action action = null;
                    SpeciesContext sc = fluxFunc.getReactants().get(0).getSpeciesContext();
                    if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
                        action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression(1));
                        jp.addAction(action);
                    }
                    sc = fluxFunc.getProducts().get(0).getSpeciesContext();
                    if (!simContext.getReactionContext().getSpeciesContextSpec(sc).isConstant()) {
                        SpeciesCountParameter spCountParam = getSpeciesCountParameter(sc);
                        action = new Action(varHash.getVariable(getMathSymbol(spCountParam, sm)), "inc", new Expression(-1));
                        jp.addAction(action);
                    }
                    subDomain.addJumpProcess(jp);
                }
            }
        }
    // end of if (simplereaction)...else if(fluxreaction)
    }
    // end of reaction step loop
    // 
    // set Variables to MathDescription all at once with the order resolved by "VariableHash"
    // 
    mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
    // set up variable initial conditions in subDomain
    SpeciesContextSpec[] scSpecs = simContext.getReactionContext().getSpeciesContextSpecs();
    for (int i = 0; i < speciesContextSpecs.length; i++) {
        // get stochastic variable by name
        SpeciesCountParameter spCountParam = getSpeciesCountParameter(speciesContextSpecs[i].getSpeciesContext());
        StructureMapping sm = simContext.getGeometryContext().getStructureMapping(speciesContextSpecs[i].getSpeciesContext().getStructure());
        String varName = getMathSymbol(spCountParam, sm);
        if (scSpecs[i].isConstant()) {
            continue;
        }
        StochVolVariable var = (StochVolVariable) mathDesc.getVariable(varName);
        // stochastic use initial number of particles
        SpeciesContextSpec.SpeciesContextSpecParameter initParm = scSpecs[i].getInitialCountParameter();
        // stochastic variables initial expression.
        if (initParm != null) {
            VarIniCondition varIni = new VarIniCount(var, new Expression(getMathSymbol(initParm, sm)));
            subDomain.addVarIniCondition(varIni);
        }
    }
    if (!mathDesc.isValid()) {
        throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
    }
}
Also used : VarIniCondition(cbit.vcell.math.VarIniCondition) MembraneMapping(cbit.vcell.mapping.MembraneMapping) LumpedKinetics(cbit.vcell.model.LumpedKinetics) MathDescription(cbit.vcell.math.MathDescription) SpeciesContextMapping(cbit.vcell.mapping.SpeciesContextMapping) Product(cbit.vcell.model.Product) FluxReaction(cbit.vcell.model.FluxReaction) SpeciesContext(cbit.vcell.model.SpeciesContext) SpeciesContextSpec(cbit.vcell.mapping.SpeciesContextSpec) Feature(cbit.vcell.model.Feature) Reactant(cbit.vcell.model.Reactant) ExpressionException(cbit.vcell.parser.ExpressionException) MappingException(cbit.vcell.mapping.MappingException) PropertyVetoException(java.beans.PropertyVetoException) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter) SubVolume(cbit.vcell.geometry.SubVolume) Vector(java.util.Vector) ModelException(cbit.vcell.model.ModelException) ReactionSpec(cbit.vcell.mapping.ReactionSpec) PropertyVetoException(java.beans.PropertyVetoException) ModelParameter(cbit.vcell.model.Model.ModelParameter) VCUnitDefinition(cbit.vcell.units.VCUnitDefinition) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) ReactionStep(cbit.vcell.model.ReactionStep) Kinetics(cbit.vcell.model.Kinetics) LumpedKinetics(cbit.vcell.model.LumpedKinetics) ReactionParticipant(cbit.vcell.model.ReactionParticipant) Action(cbit.vcell.math.Action) VariableHash(cbit.vcell.math.VariableHash) Constant(cbit.vcell.math.Constant) StructureMapping(cbit.vcell.mapping.StructureMapping) Function(cbit.vcell.math.Function) FeatureMapping(cbit.vcell.mapping.FeatureMapping) JumpProcess(cbit.vcell.math.JumpProcess) Structure(cbit.vcell.model.Structure) StochVolVariable(cbit.vcell.math.StochVolVariable) ModelUnitSystem(cbit.vcell.model.ModelUnitSystem) SimpleReaction(cbit.vcell.model.SimpleReaction) VarIniCount(cbit.vcell.math.VarIniCount) SimulationContext(cbit.vcell.mapping.SimulationContext) ElectricalStimulus(cbit.vcell.mapping.ElectricalStimulus) Expression(cbit.vcell.parser.Expression) Model(cbit.vcell.model.Model) ProxyParameter(cbit.vcell.model.ProxyParameter) Parameter(cbit.vcell.model.Parameter) KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter) ModelParameter(cbit.vcell.model.Model.ModelParameter)

Example 23 with FluxReaction

use of cbit.vcell.model.FluxReaction 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_4_8.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++) {
            int stoichiometry = reactionSteps[j].getStoichiometry(sc);
            totalSchemeMatrix.set_elem(i, j, stoichiometry);
            if (stoichiometry != 0) {
                if (!(reactionSteps[j] instanceof DummyReactionStep) && !reactionSpecs[j].isFast() && !reactionSpecs[j].isExcluded()) {
                    ReactionParticipant[] rps1 = reactionSteps[j].getReactionParticipants();
                    ReactionParticipant rp0 = null;
                    for (ReactionParticipant rp : rps1) {
                        if (rp.getSpeciesContext() == sc) {
                            rp0 = rp;
                            break;
                        }
                    }
                    Structure structure = reactionSteps[j].getStructure();
                    // 
                    if (rp0 != null) {
                        Expression reactRateExp = getReactionRateExpression(reactionSteps[j], rp0).renameBoundSymbols(mathMapping_4_8.getNameScope());
                        if ((structure instanceof Membrane) && (sc.getStructure() != structure)) {
                            Membrane membrane = (Membrane) structure;
                            MembraneMapping membraneMapping = (MembraneMapping) mathMapping_4_8.getSimulationContext().getGeometryContext().getStructureMapping(membrane);
                            Parameter fluxCorrectionParameter = mathMapping_4_8.getFluxCorrectionParameter(membraneMapping, (Feature) sc.getStructure());
                            Expression fluxCorrection = new Expression(fluxCorrectionParameter, mathMapping_4_8.getNameScope());
                            if (reactionSteps[j] instanceof FluxReaction) {
                                exp = Expression.add(exp, Expression.mult(fluxCorrection, reactRateExp));
                            // Expression.add(exp,new Expression(fluxCorrectionParameterSymbolName+"*"+expInfix));
                            } else if (reactionSteps[j] instanceof SimpleReaction) {
                                ModelUnitSystem unitSystem = mathMapping_4_8.getSimulationContext().getModel().getUnitSystem();
                                Expression unitFactor = mathMapping_4_8.getUnitFactor(unitSystem.getVolumeSubstanceUnit().divideBy(unitSystem.getMembraneSubstanceUnit()));
                                exp = Expression.add(exp, Expression.mult(fluxCorrection, unitFactor, reactRateExp));
                            // exp = Expression.add(exp,new Expression(fluxCorrectionParameterSymbolName+"*"+ReservedSymbol.KMOLE.getName()+"*"+expInfix));
                            } else {
                                throw new RuntimeException("Internal Error: expected ReactionStep " + reactionSteps[j] + " to be of type SimpleReaction or FluxReaction");
                            }
                        } else {
                            exp = Expression.add(exp, reactRateExp);
                        }
                    }
                }
            }
        }
        // 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();
// }
}
Also used : MembraneMapping(cbit.vcell.mapping.MembraneMapping) SimpleReaction(cbit.vcell.model.SimpleReaction) SpeciesContextMapping(cbit.vcell.mapping.SpeciesContextMapping) ReactionSpec(cbit.vcell.mapping.ReactionSpec) RationalNumberMatrix(cbit.vcell.matrix.RationalNumberMatrix) DummyReactionStep(cbit.vcell.mapping.DummyReactionStep) EventDummyReactionStep(cbit.vcell.mapping.EventDummyReactionStep) ParticleDummyReactionStep(cbit.vcell.mapping.ParticleDummyReactionStep) DiffusionDummyReactionStep(cbit.vcell.mapping.DiffusionDummyReactionStep) HybridDummyReactionStep(cbit.vcell.mapping.HybridDummyReactionStep) FluxReaction(cbit.vcell.model.FluxReaction) SpeciesContext(cbit.vcell.model.SpeciesContext) Expression(cbit.vcell.parser.Expression) Membrane(cbit.vcell.model.Membrane) Parameter(cbit.vcell.model.Parameter) Structure(cbit.vcell.model.Structure) ReactionParticipant(cbit.vcell.model.ReactionParticipant) ModelUnitSystem(cbit.vcell.model.ModelUnitSystem)

Example 24 with FluxReaction

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

the class MembraneStructureAnalyzer method refreshResolvedFluxes.

/**
 * This method was created in VisualAge.
 */
private void refreshResolvedFluxes() throws Exception {
    // System.out.println("MembraneStructureAnalyzer.refreshResolvedFluxes()");
    ModelUnitSystem unitSystem = mathMapping.getSimulationContext().getModel().getUnitSystem();
    GeometryContext geoContext = mathMapping.getSimulationContext().getGeometryContext();
    Vector<ResolvedFlux> resolvedFluxList = new Vector<ResolvedFlux>();
    // 
    // for each reaction, get all fluxReactions associated with this membrane
    // 
    Vector<FluxReaction> fluxList = new Vector<FluxReaction>();
    ReactionSpec[] reactionSpecs = mathMapping.getSimulationContext().getReactionContext().getReactionSpecs();
    for (int j = 0; j < reactionSpecs.length; j++) {
        if (reactionSpecs[j].isExcluded()) {
            continue;
        }
        ReactionStep rs = reactionSpecs[j].getReactionStep();
        if (rs.getStructure() != null && geoContext.getStructureMapping(rs.getStructure()).getGeometryClass() == surfaceClass) {
            if (rs instanceof FluxReaction) {
                fluxList.addElement((FluxReaction) rs);
            }
        }
    }
    // 
    for (int i = 0; i < fluxList.size(); i++) {
        FluxReaction fr = fluxList.elementAt(i);
        ReactionParticipant[] reactionParticipants = fr.getReactionParticipants();
        for (int j = 0; j < reactionParticipants.length; j++) {
            if (!(reactionParticipants[j] instanceof Reactant) && !(reactionParticipants[j] instanceof Product)) {
                continue;
            }
            ResolvedFlux rf = null;
            SpeciesContext speciesContext = reactionParticipants[j].getSpeciesContext();
            for (int k = 0; k < resolvedFluxList.size(); k++) {
                ResolvedFlux rf_tmp = resolvedFluxList.elementAt(k);
                if (rf_tmp.getSpeciesContext() == reactionParticipants[j].getSpeciesContext()) {
                    rf = rf_tmp;
                }
            }
            // 
            // if speciesContext is not "fixed" and is mapped to a volume, add flux to ResolvedFlux
            // 
            StructureMapping structureMapping = mathMapping.getSimulationContext().getGeometryContext().getStructureMapping(reactionParticipants[j].getStructure());
            if (structureMapping.getGeometryClass() == surfaceClass) {
                // flux within surface
                continue;
            }
            if (structureMapping.getGeometryClass() instanceof SubVolume && surfaceClass.isAdjacentTo((SubVolume) structureMapping.getGeometryClass())) {
                SpeciesContextSpec speciesContextSpec = mathMapping.getSimulationContext().getReactionContext().getSpeciesContextSpec(speciesContext);
                if (!speciesContextSpec.isConstant()) {
                    if (rf == null) {
                        VCUnitDefinition speciesFluxUnit = speciesContext.getUnitDefinition().multiplyBy(unitSystem.getLengthUnit()).divideBy(unitSystem.getTimeUnit());
                        rf = new ResolvedFlux(speciesContext, speciesFluxUnit);
                        resolvedFluxList.addElement(rf);
                    }
                    FeatureMapping featureMapping = (FeatureMapping) structureMapping;
                    Expression insideFluxCorrection = Expression.invert(new Expression(featureMapping.getVolumePerUnitVolumeParameter(), mathMapping.getNameScope()));
                    // 
                    if (fr.getKinetics() instanceof DistributedKinetics) {
                        KineticsParameter reactionRateParameter = ((DistributedKinetics) fr.getKinetics()).getReactionRateParameter();
                        Expression correctedReactionRate = Expression.mult(new Expression(reactionRateParameter, mathMapping.getNameScope()), insideFluxCorrection);
                        if (reactionParticipants[j] instanceof Product) {
                            if (rf.getFluxExpression().isZero()) {
                                rf.setFluxExpression(correctedReactionRate.flatten());
                            } else {
                                rf.setFluxExpression(Expression.add(rf.getFluxExpression(), correctedReactionRate.flatten()));
                            }
                        } else if (reactionParticipants[j] instanceof Reactant) {
                            if (rf.getFluxExpression().isZero()) {
                                rf.setFluxExpression(Expression.negate(correctedReactionRate).flatten());
                            } else {
                                rf.setFluxExpression(Expression.add(rf.getFluxExpression(), Expression.negate(correctedReactionRate).flatten()));
                            }
                        } else {
                            throw new RuntimeException("expected either FluxReactant or FluxProduct");
                        }
                    } 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.getFluxExpression().bindExpression(mathMapping);
                }
            }
        }
    }
    // 
    for (int i = 0; i < reactionSpecs.length; i++) {
        if (reactionSpecs[i].isExcluded()) {
            continue;
        }
        ReactionStep rs = reactionSpecs[i].getReactionStep();
        if (rs.getStructure() != null && geoContext.getStructureMapping(rs.getStructure()).getGeometryClass() == surfaceClass) {
            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 rpSCS = mathMapping.getSimulationContext().getReactionContext().getSpeciesContextSpec(rp_Array[k].getSpeciesContext());
                        StructureMapping rpSM = mathMapping.getSimulationContext().getGeometryContext().getStructureMapping(rp_Array[k].getStructure());
                        if (rpSM.getGeometryClass() instanceof SubVolume && !rpSCS.isConstant()) {
                            // 
                            // 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.getSpeciesContext() == rp_Array[k].getSpeciesContext()) {
                                    rf = rf_tmp;
                                }
                            }
                            if (rf == null) {
                                VCUnitDefinition speciesFluxUnit = rp_Array[k].getSpeciesContext().getUnitDefinition().multiplyBy(unitSystem.getLengthUnit()).divideBy(unitSystem.getTimeUnit());
                                rf = new ResolvedFlux(rp_Array[k].getSpeciesContext(), speciesFluxUnit);
                                resolvedFluxList.addElement(rf);
                            }
                            if (rpSM.getGeometryClass() instanceof SubVolume && surfaceClass.isAdjacentTo((SubVolume) rpSM.getGeometryClass())) {
                                // 
                                // for binding on inside or outside, add to ResolvedFlux.flux
                                // 
                                Expression fluxRateExpression = getCorrectedRateExpression(sr, rp_Array[k], RateType.ResolvedFluxRate).renameBoundSymbols(mathMapping.getNameScope());
                                if (rf.getFluxExpression().isZero()) {
                                    rf.setFluxExpression(fluxRateExpression);
                                } else {
                                    rf.setFluxExpression(Expression.add(rf.getFluxExpression(), fluxRateExpression));
                                }
                                rf.getFluxExpression().bindExpression(mathMapping);
                            } else {
                                String structureName = ((rs.getStructure() != null) ? (rs.getStructure().getName()) : ("<null>"));
                                throw new Exception("In Application '" + mathMapping.getSimulationContext().getName() + "', SpeciesContext '" + rp_Array[k].getSpeciesContext().getName() + "' is not mapped adjacent to structure '" + structureName + "' but reacts there");
                            }
                        }
                    }
                }
            }
        }
    }
    // 
    if (resolvedFluxList.size() > 0) {
        resolvedFluxes = new ResolvedFlux[resolvedFluxList.size()];
        resolvedFluxList.copyInto(resolvedFluxes);
    } else {
        resolvedFluxes = null;
    }
}
Also used : LumpedKinetics(cbit.vcell.model.LumpedKinetics) Product(cbit.vcell.model.Product) FluxReaction(cbit.vcell.model.FluxReaction) SpeciesContext(cbit.vcell.model.SpeciesContext) Reactant(cbit.vcell.model.Reactant) KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter) SubVolume(cbit.vcell.geometry.SubVolume) Vector(java.util.Vector) ModelUnitSystem(cbit.vcell.model.ModelUnitSystem) DistributedKinetics(cbit.vcell.model.DistributedKinetics) SimpleReaction(cbit.vcell.model.SimpleReaction) VCUnitDefinition(cbit.vcell.units.VCUnitDefinition) Expression(cbit.vcell.parser.Expression) ReactionStep(cbit.vcell.model.ReactionStep) ReactionParticipant(cbit.vcell.model.ReactionParticipant)

Example 25 with FluxReaction

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

the class ReactionCartoonFull method refreshAll.

@Override
protected void refreshAll(boolean reallocateShapes) {
    try {
        if (getModel() == null || getStructureSuite() == null) {
            return;
        }
        System.out.println("ReactionCartoonFull, RefreshAll()");
        for (Structure structure : structureSuite.getStructures()) {
            Diagram diagram = getModel().getDiagram(structure);
            if (diagram != null) {
                // Maintain consistency between rule participant nodes, signatures and
                // species pattern when a molecule is being modified.
                rebindAll(diagram);
            }
        }
        // calculate species context weight (number of reactions for which it's a participant)
        Map<SpeciesContext, Integer> scWeightMap = new HashMap<>();
        // all the species contexts that are catalysts
        Set<SpeciesContext> scCatalystSet = new HashSet<>();
        // calculate species context length (number of species patterns it contains, 1 if has no species patterns)
        for (ReactionStep rs : getModel().getReactionSteps()) {
            ReactionParticipant[] rpList = rs.getReactionParticipants();
            for (int i = 0; i < rpList.length; i++) {
                ReactionParticipant rp = rpList[i];
                SpeciesContext sc = rp.getSpeciesContext();
                // int increment = rp.getStoichiometry();
                int increment = 1;
                if (rp instanceof Catalyst) {
                    scCatalystSet.add(sc);
                }
                if (scWeightMap.containsKey(sc)) {
                    int weight = scWeightMap.get(sc);
                    weight += increment;
                    scWeightMap.put(sc, weight);
                } else {
                    scWeightMap.put(sc, increment);
                }
            }
        }
        Set<Shape> unwantedShapes = new HashSet<Shape>();
        Set<RuleParticipantSignature> unwantedSignatures = new HashSet<RuleParticipantSignature>();
        unwantedShapes.addAll(getShapes());
        unwantedSignatures.addAll(ruleParticipantSignatures);
        ContainerContainerShape containerShape = (ContainerContainerShape) getShapeFromModelObject(getModel());
        List<ReactionContainerShape> reactionContainerShapeList = new ArrayList<ReactionContainerShape>();
        List<Structure> structureList = new ArrayList<Structure>(getStructureSuite().getStructures());
        // create all ReactionContainerShapes (one for each Structure)
        for (Structure structure : structureList) {
            if (structure instanceof Membrane) {
                Membrane membrane = (Membrane) structure;
                ReactionContainerShape membraneShape = (ReactionContainerShape) getShapeFromModelObject(membrane);
                if (membraneShape == null) {
                    membraneShape = new ReactionContainerShape(membrane, structureSuite, this);
                    addShape(membraneShape);
                    membrane.getMembraneVoltage().removePropertyChangeListener(this);
                    membrane.getMembraneVoltage().addPropertyChangeListener(this);
                } else {
                    membraneShape.setStructureSuite(structureSuite);
                }
                membrane.removePropertyChangeListener(this);
                membrane.addPropertyChangeListener(this);
                membraneShape.refreshLabel();
                unwantedShapes.remove(membraneShape);
                reactionContainerShapeList.add(membraneShape);
            } else if (structure instanceof Feature) {
                Feature feature = (Feature) structure;
                ReactionContainerShape featureShape = (ReactionContainerShape) getShapeFromModelObject(feature);
                if (featureShape == null) {
                    featureShape = new ReactionContainerShape(feature, structureSuite, this);
                    addShape(featureShape);
                } else {
                    featureShape.setStructureSuite(structureSuite);
                }
                feature.removePropertyChangeListener(this);
                feature.addPropertyChangeListener(this);
                featureShape.refreshLabel();
                unwantedShapes.remove(featureShape);
                reactionContainerShapeList.add(featureShape);
            }
        }
        if (containerShape == null) {
            containerShape = new ContainerContainerShape(this, getModel(), reactionContainerShapeList);
            addShape(containerShape);
        } else {
            containerShape.setReactionContainerShapeList(reactionContainerShapeList);
        }
        containerShape.refreshLabel();
        unwantedShapes.remove(containerShape);
        // add all species context shapes within the structures
        for (Structure structure : getStructureSuite().getStructures()) {
            ReactionContainerShape reactionContainerShape = (ReactionContainerShape) getShapeFromModelObject(structure);
            structure.removePropertyChangeListener(this);
            structure.addPropertyChangeListener(this);
            for (SpeciesContext structSpeciesContext : getModel().getSpeciesContexts(structure)) {
                SpeciesContextShape ss = (SpeciesContextShape) getShapeFromModelObject(structSpeciesContext);
                if (ss == null) {
                    ss = new SpeciesContextShape(structSpeciesContext, this);
                    ss.truncateLabelName(false);
                    structSpeciesContext.getSpecies().removePropertyChangeListener(this);
                    structSpeciesContext.getSpecies().addPropertyChangeListener(this);
                    reactionContainerShape.addChildShape(ss);
                    addShape(ss);
                    ss.getSpaceManager().setRelPos(reactionContainerShape.getRandomPosition());
                }
                if (speciesSizeOption == SpeciesSizeOptions.weight) {
                    // this number sets the diameter of the shape
                    Integer weight = scWeightMap.get(structSpeciesContext);
                    if (weight != null) {
                        // we cap the diameter of the shape to something reasonable
                        weight = Math.min(weight, 16);
                    }
                    ss.setFilters(highlightCatalystOption ? scCatalystSet.contains(structSpeciesContext) : false, weight);
                } else if (speciesSizeOption == SpeciesSizeOptions.length) {
                    Integer length = null;
                    if (structSpeciesContext.getSpeciesPattern() != null && !structSpeciesContext.getSpeciesPattern().getMolecularTypePatterns().isEmpty()) {
                        length = structSpeciesContext.getSpeciesPattern().getMolecularTypePatterns().size() * 2;
                        length = Math.min(length, 16);
                    }
                    ss.setFilters(highlightCatalystOption ? scCatalystSet.contains(structSpeciesContext) : false, length);
                } else {
                    ss.setFilters(highlightCatalystOption ? scCatalystSet.contains(structSpeciesContext) : false, null);
                }
                structSpeciesContext.removePropertyChangeListener(this);
                structSpeciesContext.addPropertyChangeListener(this);
                ss.refreshLabel();
                unwantedShapes.remove(ss);
            }
        }
        // add all reactionSteps that are in this structure (ReactionContainerShape), and draw the lines
        getModel().removePropertyChangeListener(this);
        getModel().addPropertyChangeListener(this);
        // 
        for (ReactionRule rr : getModel().getRbmModelContainer().getReactionRuleList()) {
            rr.removePropertyChangeListener(this);
            rr.addPropertyChangeListener(this);
            Structure structure = rr.getStructure();
            if (getStructureSuite().areReactionsShownFor(structure)) {
                ReactionContainerShape reactionContainerShape = (ReactionContainerShape) getShapeFromModelObject(structure);
                ReactionRuleFullDiagramShape rrShape = (ReactionRuleFullDiagramShape) getShapeFromModelObject(rr);
                if (rrShape == null) {
                    rrShape = new ReactionRuleFullDiagramShape(rr, this);
                    addShape(rrShape);
                    rrShape.getSpaceManager().setRelPos(reactionContainerShape.getRandomPosition());
                    reactionContainerShape.addChildShape(rrShape);
                    rrShape.getSpaceManager().setRelPos(reactionContainerShape.getRandomPosition());
                }
                rrShape.refreshLabel();
                unwantedShapes.remove(rrShape);
                // 
                // add reaction participants as edges and SignatureShapes as needed
                // 
                List<ReactionRuleParticipant> participants = rr.getReactionRuleParticipants();
                List<RuleParticipantEdgeDiagramShape> ruleEdges = new ArrayList<>();
                for (ReactionRuleParticipant participant : participants) {
                    participant.getSpeciesPattern().removePropertyChangeListener(this);
                    participant.getSpeciesPattern().addPropertyChangeListener(this);
                    Structure speciesStructure = participant.getStructure();
                    Structure reactionStructure = rr.getStructure();
                    if (getStructureSuite().getStructures().contains(speciesStructure) && getStructureSuite().areReactionsShownFor(reactionStructure)) {
                        // 
                        // find existing RuleParticipantSignatureShape in cartoon
                        // 
                        RuleParticipantLongSignature ruleParticipantLongSignature = null;
                        for (RuleParticipantSignature signature : ruleParticipantSignatures) {
                            if (signature instanceof RuleParticipantLongSignature && signature.getStructure() == participant.getStructure() && signature.compareByCriteria(participant.getSpeciesPattern(), GroupingCriteria.full)) {
                                ruleParticipantLongSignature = (RuleParticipantLongSignature) signature;
                                break;
                            }
                        }
                        for (RuleParticipantSignature signature : ruleParticipantSignatures) {
                            if (signature instanceof RuleParticipantShortSignature && signature.getStructure() == participant.getStructure()) {
                                System.out.println("ReactionCartoonFull, refreshAll(), RuleParticipantShortSignature");
                            }
                        }
                        // 
                        // if didn't find signature in cartoons list of signatures, then create one (and create a shape for it).
                        // 
                        RuleParticipantSignatureFullDiagramShape signatureShape = null;
                        if (ruleParticipantLongSignature == null) {
                            ruleParticipantLongSignature = RuleParticipantLongSignature.fromReactionRuleParticipant(participant, this);
                            ruleParticipantSignatures.add(ruleParticipantLongSignature);
                            signatureShape = new RuleParticipantSignatureFullDiagramShape(ruleParticipantLongSignature, this);
                            addShape(signatureShape);
                            ReactionContainerShape participantContainerShape = (ReactionContainerShape) getShapeFromModelObject(participant.getStructure());
                            signatureShape.getSpaceManager().setRelPos(participantContainerShape.getRandomPosition());
                            participantContainerShape.addChildShape(signatureShape);
                            signatureShape.getSpaceManager().setRelPos(participantContainerShape.getRandomPosition());
                        } else {
                            signatureShape = (RuleParticipantSignatureFullDiagramShape) getShapeFromModelObject(ruleParticipantLongSignature);
                        }
                        unwantedShapes.remove(signatureShape);
                        unwantedSignatures.remove(ruleParticipantLongSignature);
                        signatureShape.refreshLabel();
                        signatureShape.setVisible(true);
                        // 
                        // add edge for ReactionRuleParticipant if not already present.
                        // 
                        RuleParticipantEdgeDiagramShape ruleParticipantShape = (RuleParticipantEdgeDiagramShape) getShapeFromModelObject(participant);
                        if (ruleParticipantShape == null || ruleParticipantShape.getRuleParticipantSignatureShape() != signatureShape) {
                            if (participant instanceof ReactantPattern && signatureShape.isVisible()) {
                                ruleParticipantShape = new ReactantPatternEdgeDiagramShape((ReactantPattern) participant, rrShape, signatureShape, this);
                            } else if (participant instanceof ProductPattern && signatureShape.isVisible()) {
                                ruleParticipantShape = new ProductPatternEdgeDiagramShape((ProductPattern) participant, rrShape, signatureShape, this);
                            } else {
                                throw new RuntimeException("unsupported ReactionRuleParticipant " + participant.getClass());
                            }
                            addShape(ruleParticipantShape);
                        }
                        if (!containerShape.getChildren().contains(ruleParticipantShape)) {
                            containerShape.addChildShape(ruleParticipantShape);
                        }
                        unwantedShapes.remove(ruleParticipantShape);
                        ruleParticipantShape.refreshLabel();
                        // all the edges for this rule
                        ruleEdges.add(ruleParticipantShape);
                    }
                }
                // a product edge (a closed loop) between the rule diagram shape and the signature diagram shape
                for (RuleParticipantEdgeDiagramShape ours : ruleEdges) {
                    // reset them all
                    ours.setSibling(false);
                }
                for (RuleParticipantEdgeDiagramShape ours : ruleEdges) {
                    for (RuleParticipantEdgeDiagramShape theirs : ruleEdges) {
                        if (ours == theirs) {
                            // don't compare with self
                            continue;
                        }
                        if (ours.getRuleParticipantSignatureShape() == theirs.getRuleParticipantSignatureShape()) {
                            ours.setSibling(true);
                            theirs.setSibling(true);
                        }
                    }
                }
            }
        }
        ruleParticipantSignatures.removeAll(unwantedSignatures);
        for (ReactionStep reactionStep : getModel().getReactionSteps()) {
            reactionStep.removePropertyChangeListener(this);
            reactionStep.addPropertyChangeListener(this);
            Structure structure = reactionStep.getStructure();
            if (getStructureSuite().areReactionsShownFor(structure)) {
                ReactionContainerShape reactionContainerShape = (ReactionContainerShape) getShapeFromModelObject(structure);
                if (reactionContainerShape == null) {
                    System.out.println("Reaction container shape is null for structure " + structure + " for reaction step " + reactionStep);
                }
                ReactionStepShape reactionStepShape = (ReactionStepShape) getShapeFromModelObject(reactionStep);
                if (reactionStepShape == null) {
                    if (reactionStep instanceof SimpleReaction) {
                        reactionStepShape = new SimpleReactionShape((SimpleReaction) reactionStep, this);
                    } else if (reactionStep instanceof FluxReaction) {
                        reactionStepShape = new FluxReactionShape((FluxReaction) reactionStep, this);
                    } else {
                        throw new RuntimeException("unknown type of ReactionStep '" + reactionStep.getClass().toString());
                    }
                    addShape(reactionStepShape);
                    reactionStepShape.getSpaceManager().setRelPos(reactionContainerShape.getRandomPosition());
                    reactionContainerShape.addChildShape(reactionStepShape);
                    reactionStepShape.getSpaceManager().setRelPos(reactionContainerShape.getRandomPosition());
                }
                reactionStepShape.refreshLabel();
                unwantedShapes.remove(reactionStepShape);
                // add reaction participants as edges
                for (ReactionParticipant participant : reactionStep.getReactionParticipants()) {
                    participant.removePropertyChangeListener(this);
                    participant.addPropertyChangeListener(this);
                    Structure speciesStructure = participant.getStructure();
                    Structure reactionStructure = reactionStep.getStructure();
                    if (getStructureSuite().getStructures().contains(speciesStructure) && getStructureSuite().areReactionsShownFor(reactionStructure)) {
                        SpeciesContext speciesContext = getModel().getSpeciesContext(participant.getSpecies(), speciesStructure);
                        // add speciesContextShapes that are not in this structure, but are referenced from the reactionParticipants
                        // these are only when reactionParticipants are from features that are outside of the membrane being displayed
                        SpeciesContextShape speciesContextShape = (SpeciesContextShape) getShapeFromModelObject(speciesContext);
                        if (speciesContextShape == null) {
                            speciesContextShape = new SpeciesContextShape(speciesContext, this);
                            speciesContextShape.truncateLabelName(false);
                            reactionContainerShape.addChildShape(speciesContextShape);
                            addShape(speciesContextShape);
                            speciesContextShape.getSpaceManager().setRelPos(reactionContainerShape.getRandomPosition());
                        }
                        speciesContextShape.refreshLabel();
                        unwantedShapes.remove(speciesContextShape);
                        ReactionParticipantShape reactionParticipantShape = (ReactionParticipantShape) getShapeFromModelObject(participant);
                        if (reactionParticipantShape == null) {
                            if (participant instanceof Reactant) {
                                reactionParticipantShape = new ReactantShape((Reactant) participant, reactionStepShape, speciesContextShape, this);
                            } else if (participant instanceof Product) {
                                reactionParticipantShape = new ProductShape((Product) participant, reactionStepShape, speciesContextShape, this);
                            } else if (participant instanceof Catalyst) {
                                reactionParticipantShape = new CatalystShape((Catalyst) participant, reactionStepShape, speciesContextShape, this);
                            } else {
                                throw new RuntimeException("unsupported ReactionParticipant " + participant.getClass());
                            }
                            addShape(reactionParticipantShape);
                        }
                        if (!containerShape.getChildren().contains(reactionParticipantShape)) {
                            containerShape.addChildShape(reactionParticipantShape);
                        }
                        unwantedShapes.remove(reactionParticipantShape);
                        reactionParticipantShape.refreshLabel();
                    }
                }
            }
        }
        for (Shape unwantedShape : unwantedShapes) {
            removeShape(unwantedShape);
        }
        // update diagrams
        for (Structure structure : structureSuite.getStructures()) {
            Diagram diagram = getModel().getDiagram(structure);
            if (diagram != null) {
                applyDefaults(diagram);
            }
        }
        fireGraphChanged(new GraphEvent(this));
    } catch (Throwable e) {
        handleException(e);
    }
}
Also used : HashMap(java.util.HashMap) ArrayList(java.util.ArrayList) Product(cbit.vcell.model.Product) FluxReaction(cbit.vcell.model.FluxReaction) SpeciesContext(cbit.vcell.model.SpeciesContext) Feature(cbit.vcell.model.Feature) Reactant(cbit.vcell.model.Reactant) HashSet(java.util.HashSet) ReactionStep(cbit.vcell.model.ReactionStep) ReactionRuleParticipant(cbit.vcell.model.ReactionRuleParticipant) ReactionParticipant(cbit.vcell.model.ReactionParticipant) Catalyst(cbit.vcell.model.Catalyst) RuleParticipantSignature(cbit.vcell.model.RuleParticipantSignature) Shape(cbit.gui.graph.Shape) GraphEvent(cbit.gui.graph.GraphEvent) RuleParticipantShortSignature(cbit.vcell.model.RuleParticipantShortSignature) Membrane(cbit.vcell.model.Membrane) Structure(cbit.vcell.model.Structure) ReactantPattern(cbit.vcell.model.ReactantPattern) RuleParticipantLongSignature(cbit.vcell.model.RuleParticipantLongSignature) SimpleReaction(cbit.vcell.model.SimpleReaction) ReactionRule(cbit.vcell.model.ReactionRule) ProductPattern(cbit.vcell.model.ProductPattern) Point(java.awt.Point) Diagram(cbit.vcell.model.Diagram)

Aggregations

FluxReaction (cbit.vcell.model.FluxReaction)35 SimpleReaction (cbit.vcell.model.SimpleReaction)26 SpeciesContext (cbit.vcell.model.SpeciesContext)17 Structure (cbit.vcell.model.Structure)17 Membrane (cbit.vcell.model.Membrane)16 ReactionParticipant (cbit.vcell.model.ReactionParticipant)16 ArrayList (java.util.ArrayList)14 ReactionStep (cbit.vcell.model.ReactionStep)13 Shape (cbit.gui.graph.Shape)12 Feature (cbit.vcell.model.Feature)12 Reactant (cbit.vcell.model.Reactant)12 Expression (cbit.vcell.parser.Expression)11 Product (cbit.vcell.model.Product)10 Point (java.awt.Point)10 PropertyVetoException (java.beans.PropertyVetoException)9 KineticsParameter (cbit.vcell.model.Kinetics.KineticsParameter)8 RuleParticipantSignature (cbit.vcell.model.RuleParticipantSignature)8 Catalyst (cbit.vcell.model.Catalyst)7 ModelUnitSystem (cbit.vcell.model.ModelUnitSystem)7 LumpedKinetics (cbit.vcell.model.LumpedKinetics)6