Examples with SubDomain cbit.vcell.math.SubDomain used on opensource projects

Example 36 with SubDomain

use of cbit.vcell.math.SubDomain in project vcell by virtualcell.

the class MathDescriptionTreeModel method createBaseTree.

/**
 * Insert the method's description here.
 * Creation date: (11/28/00 1:06:51 PM)
 * @return cbit.vcell.desktop.BioModelNode
 * @param docManager cbit.vcell.clientdb.DocumentManager
 */
private BioModelNode createBaseTree() {
    if (getMathDescription() == null) {
        return new BioModelNode(" ", false);
    }
    // 
    // make root node
    // 
    BioModelNode rootNode = new BioModelNode("math description", true);
    // 
    // create subTree of Constants
    // 
    BioModelNode constantRootNode = new BioModelNode("constants", true);
    Enumeration<Constant> enum1 = getMathDescription().getConstants();
    while (enum1.hasMoreElements()) {
        Constant constant = enum1.nextElement();
        BioModelNode constantNode = new BioModelNode(constant, false);
        constantRootNode.add(constantNode);
    }
    if (constantRootNode.getChildCount() > 0) {
        rootNode.add(constantRootNode);
    }
    // 
    // create subTree of Functions
    // 
    BioModelNode functionRootNode = new BioModelNode("functions", true);
    Enumeration<Function> enum2 = getMathDescription().getFunctions();
    while (enum2.hasMoreElements()) {
        Function function = enum2.nextElement();
        BioModelNode functionNode = new BioModelNode(function, false);
        functionRootNode.add(functionNode);
    }
    if (functionRootNode.getChildCount() > 0) {
        rootNode.add(functionRootNode);
    }
    // 
    // create subTree of VolumeSubDomains and MembraneSubDomains
    // 
    BioModelNode volumeRootNode = new BioModelNode("volume domains", true);
    BioModelNode membraneRootNode = new BioModelNode("membrane domains", true);
    Enumeration<SubDomain> enum3 = getMathDescription().getSubDomains();
    while (enum3.hasMoreElements()) {
        SubDomain subDomain = enum3.nextElement();
        if (subDomain instanceof cbit.vcell.math.CompartmentSubDomain) {
            CompartmentSubDomain volumeSubDomain = (CompartmentSubDomain) subDomain;
            BioModelNode volumeSubDomainNode = new BioModelNode(volumeSubDomain, true);
            if (// stochastic subtree
            getMathDescription().isNonSpatialStoch()) {
                // add stoch variable initial conditions
                BioModelNode varIniConditionNode = new BioModelNode("variable_initial_conditions", true);
                for (VarIniCondition varIni : volumeSubDomain.getVarIniConditions()) {
                    BioModelNode varIniNode = new BioModelNode(varIni, false);
                    varIniConditionNode.add(varIniNode);
                }
                volumeSubDomainNode.add(varIniConditionNode);
                // add jump processes
                for (JumpProcess jp : volumeSubDomain.getJumpProcesses()) {
                    BioModelNode jpNode = new BioModelNode(jp, true);
                    // add probability rate.
                    String probRate = "P_" + jp.getName();
                    BioModelNode prNode = new BioModelNode("probability_rate = " + probRate, false);
                    jpNode.add(prNode);
                    // add Actions
                    Enumeration<Action> actions = Collections.enumeration(jp.getActions());
                    while (actions.hasMoreElements()) {
                        Action action = actions.nextElement();
                        BioModelNode actionNode = new BioModelNode(action, false);
                        jpNode.add(actionNode);
                    }
                    volumeSubDomainNode.add(jpNode);
                }
            } else // non-stochastic subtree
            {
                // 
                // add equation children
                // 
                Enumeration<Equation> eqnEnum = volumeSubDomain.getEquations();
                while (eqnEnum.hasMoreElements()) {
                    Equation equation = eqnEnum.nextElement();
                    BioModelNode equationNode = new BioModelNode(equation, false);
                    volumeSubDomainNode.add(equationNode);
                }
                // 
                // add fast system
                // 
                FastSystem fastSystem = volumeSubDomain.getFastSystem();
                if (fastSystem != null) {
                    BioModelNode fsNode = new BioModelNode(fastSystem, true);
                    Enumeration<FastInvariant> enumFI = fastSystem.getFastInvariants();
                    while (enumFI.hasMoreElements()) {
                        FastInvariant fi = enumFI.nextElement();
                        fsNode.add(new BioModelNode(fi, false));
                    }
                    Enumeration<FastRate> enumFR = fastSystem.getFastRates();
                    while (enumFR.hasMoreElements()) {
                        FastRate fr = enumFR.nextElement();
                        fsNode.add(new BioModelNode(fr, false));
                    }
                    volumeSubDomainNode.add(fsNode);
                }
            }
            volumeRootNode.add(volumeSubDomainNode);
        } else if (subDomain instanceof MembraneSubDomain) {
            MembraneSubDomain membraneSubDomain = (MembraneSubDomain) subDomain;
            BioModelNode membraneSubDomainNode = new BioModelNode(membraneSubDomain, true);
            // 
            // add equation children
            // 
            Enumeration<Equation> eqnEnum = membraneSubDomain.getEquations();
            while (eqnEnum.hasMoreElements()) {
                Equation equation = eqnEnum.nextElement();
                BioModelNode equationNode = new BioModelNode(equation, false);
                membraneSubDomainNode.add(equationNode);
            }
            // 
            // add jump condition children
            // 
            Enumeration<JumpCondition> jcEnum = membraneSubDomain.getJumpConditions();
            while (jcEnum.hasMoreElements()) {
                JumpCondition jumpCondition = jcEnum.nextElement();
                BioModelNode jcNode = new BioModelNode(jumpCondition, false);
                membraneSubDomainNode.add(jcNode);
            }
            // 
            // add fast system
            // 
            FastSystem fastSystem = membraneSubDomain.getFastSystem();
            if (fastSystem != null) {
                BioModelNode fsNode = new BioModelNode(fastSystem, true);
                Enumeration<FastInvariant> enumFI = fastSystem.getFastInvariants();
                while (enumFI.hasMoreElements()) {
                    FastInvariant fi = enumFI.nextElement();
                    fsNode.add(new BioModelNode(fi, false));
                }
                Enumeration<FastRate> enumFR = fastSystem.getFastRates();
                while (enumFR.hasMoreElements()) {
                    FastRate fr = enumFR.nextElement();
                    fsNode.add(new BioModelNode(fr, false));
                }
                membraneSubDomainNode.add(fsNode);
            }
            membraneRootNode.add(membraneSubDomainNode);
        }
    }
    if (volumeRootNode.getChildCount() > 0) {
        rootNode.add(volumeRootNode);
    }
    if (membraneRootNode.getChildCount() > 0) {
        rootNode.add(membraneRootNode);
    }
    return rootNode;
}

Example 37 with SubDomain

use of cbit.vcell.math.SubDomain 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());
    }
}

Example 38 with SubDomain

use of cbit.vcell.math.SubDomain in project vcell by virtualcell.

the class ParticleMathMapping method combineHybrid.

private void combineHybrid() throws MappingException, ExpressionException, MatrixException, MathException, ModelException {
    ArrayList<SpeciesContext> continuousSpecies = new ArrayList<SpeciesContext>();
    ArrayList<ParticleVariable> continuousSpeciesParticleVars = new ArrayList<ParticleVariable>();
    ArrayList<SpeciesContext> stochSpecies = new ArrayList<SpeciesContext>();
    // 
    // categorize speciesContexts as continuous and stochastic
    // 
    SpeciesContextSpec[] scsArray = getSimulationContext().getReactionContext().getSpeciesContextSpecs();
    continuousSpecies = new ArrayList<SpeciesContext>();
    stochSpecies = new ArrayList<SpeciesContext>();
    for (SpeciesContextSpec speciesContextSpec : scsArray) {
        if (!getSimulationContext().isStoch() || speciesContextSpec.isForceContinuous()) {
            continuousSpecies.add(speciesContextSpec.getSpeciesContext());
            Variable variable = getMathSymbolMapping().getVariable(speciesContextSpec.getSpeciesContext());
            if (variable instanceof ParticleVariable) {
                continuousSpeciesParticleVars.add((ParticleVariable) variable);
            }
        } else {
            stochSpecies.add(speciesContextSpec.getSpeciesContext());
        }
    }
    if (continuousSpecies.isEmpty()) {
        return;
    }
    // 
    // create continuous mathDescription ... add stochastic variables and processes to the continuous Math and use this.
    // 
    DiffEquMathMapping mathMapping = new DiffEquMathMapping(getSimulationContext(), callback, networkGenerationRequirements);
    mathMapping.refresh(null);
    MathDescription contMathDesc = mathMapping.getMathDescription();
    // 
    // get list of all continuous variables
    // 
    HashMap<String, Variable> allContinuousVars = new HashMap<String, Variable>();
    Enumeration<Variable> enumVar = contMathDesc.getVariables();
    while (enumVar.hasMoreElements()) {
        Variable var = enumVar.nextElement();
        allContinuousVars.put(var.getName(), var);
    }
    // 
    // replace those continuous variables and equations for stochastic speciesContexts
    // with the particleVariables and particleProperties
    // (ParticleJumpProcesses removed later)
    // 
    ModelUnitSystem unitSystem = getSimulationContext().getModel().getUnitSystem();
    for (SpeciesContext stochSpeciesContext : stochSpecies) {
        Variable contVar = mathMapping.getMathSymbolMapping().getVariable(stochSpeciesContext);
        Variable stochVar = getMathSymbolMapping().getVariable(stochSpeciesContext);
        allContinuousVars.put(stochVar.getName(), stochVar);
        // 
        // replace continuous "concentration" VolVariable/MemVariable for this particle with a Function for concentration
        // 
        allContinuousVars.remove(contVar);
        VCUnitDefinition sizeUnit = unitSystem.getLengthUnit().raiseTo(new RationalNumber(stochSpeciesContext.getStructure().getDimension()));
        VCUnitDefinition stochasticDensityUnit = unitSystem.getStochasticSubstanceUnit().divideBy(sizeUnit);
        VCUnitDefinition continuousDensityUnit = unitSystem.getConcentrationUnit(stochSpeciesContext.getStructure());
        if (stochasticDensityUnit.isEquivalent(continuousDensityUnit)) {
            allContinuousVars.put(contVar.getName(), new Function(contVar.getName(), new Expression(stochVar, getNameScope()), contVar.getDomain()));
        } else {
            Expression conversionFactorExp = getUnitFactor(continuousDensityUnit.divideBy(stochasticDensityUnit));
            allContinuousVars.put(contVar.getName(), new Function(contVar.getName(), Expression.mult(new Expression(stochVar, getNameScope()), conversionFactorExp), contVar.getDomain()));
        }
        // 
        // remove continuous equation
        // 
        Enumeration<SubDomain> contSubDomains = contMathDesc.getSubDomains();
        while (contSubDomains.hasMoreElements()) {
            SubDomain contSubDomain = contSubDomains.nextElement();
            contSubDomain.removeEquation(contVar);
            if (contSubDomain instanceof MembraneSubDomain) {
                ((MembraneSubDomain) contSubDomain).removeJumpCondition(contVar);
            }
        }
        // 
        // remove all continuous variables for speciesContextSpec parameters (e.g. initial conditions, diffusion rates, boundary conditions, velocities)
        // 
        SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(stochSpeciesContext);
        Parameter[] scsParameters = scs.getParameters();
        for (Parameter parameter : scsParameters) {
            Variable continuousScsParmVariable = mathMapping.getMathSymbolMapping().getVariable(parameter);
            allContinuousVars.remove(continuousScsParmVariable);
        }
        // 
        // copy ParticleJumpProcess and ParticleProperties to the continuous math
        // 
        SubDomain contSubDomain = contMathDesc.getSubDomain(contVar.getDomain().getName());
        SubDomain stochSubDomain = mathDesc.getSubDomain(stochVar.getDomain().getName());
        ParticleProperties particleProperties = stochSubDomain.getParticleProperties(stochVar);
        contSubDomain.addParticleProperties(particleProperties);
    }
    // 
    // add all ParticleJumpProcesses to the continuous model
    // 
    Enumeration<SubDomain> enumStochSubdomains = mathDesc.getSubDomains();
    while (enumStochSubdomains.hasMoreElements()) {
        SubDomain stochSubdomain = enumStochSubdomains.nextElement();
        SubDomain contSubdomain = contMathDesc.getSubDomain(stochSubdomain.getName());
        for (ParticleJumpProcess particleJumpProcess : stochSubdomain.getParticleJumpProcesses()) {
            // 
            // modify "selection list" (particleVariables), probability rate, and actions if referenced particleVariable is to be "forced continuous"
            // 
            ParticleVariable[] selectedParticles = particleJumpProcess.getParticleVariables();
            for (ParticleVariable particleVariable : selectedParticles) {
                if (continuousSpeciesParticleVars.contains(particleVariable)) {
                    particleJumpProcess.remove(particleVariable);
                    JumpProcessRateDefinition jumpProcessRateDefinition = particleJumpProcess.getParticleRateDefinition();
                    if (jumpProcessRateDefinition instanceof MacroscopicRateConstant) {
                        MacroscopicRateConstant macroscopicRateConstant = (MacroscopicRateConstant) jumpProcessRateDefinition;
                        macroscopicRateConstant.setExpression(Expression.mult(macroscopicRateConstant.getExpression(), new Expression(particleVariable, null)));
                    } else if (jumpProcessRateDefinition instanceof InteractionRadius) {
                        throw new MappingException("cannot adjust interaction radius for reaction process " + particleJumpProcess.getName() + ", particle " + particleVariable.getName() + " is continuous");
                    } else {
                        throw new MappingException("rate definition type " + jumpProcessRateDefinition.getClass().getSimpleName() + " not yet implemented for hybrid PDE/Particle math generation");
                    }
                }
                Iterator<Action> iterAction = particleJumpProcess.getActions().iterator();
                while (iterAction.hasNext()) {
                    Action action = iterAction.next();
                    if (continuousSpeciesParticleVars.contains(action.getVar())) {
                        iterAction.remove();
                    }
                }
            }
            if (!particleJumpProcess.getActions().isEmpty()) {
                contSubdomain.addParticleJumpProcess(particleJumpProcess);
            }
        }
    }
    // 
    for (MathMappingParameter mathMappingParameter : fieldMathMappingParameters) {
        if (mathMappingParameter instanceof UnitFactorParameter) {
            String name = mathMappingParameter.getName();
            if (!allContinuousVars.containsKey(name)) {
                allContinuousVars.put(name, newFunctionOrConstant(name, mathMappingParameter.getExpression(), null));
            }
        }
    }
    // 
    // add constants and functions from the particle math that aren't already defined in the continuous math
    // 
    Enumeration<Variable> enumVars = mathDesc.getVariables();
    while (enumVars.hasMoreElements()) {
        Variable var = enumVars.nextElement();
        if (var instanceof Constant || var instanceof Function) {
            String name = var.getName();
            if (!allContinuousVars.containsKey(name)) {
                allContinuousVars.put(name, var);
            }
        }
    }
    contMathDesc.setAllVariables(allContinuousVars.values().toArray(new Variable[0]));
    mathDesc = contMathDesc;
    // 
    for (int i = 0; i < fieldMathMappingParameters.length; i++) {
        if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
            GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
            Variable variable = newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass());
            if (mathDesc.getVariable(variable.getName()) == null) {
                mathDesc.addVariable(variable);
            }
        }
    }
    if (!mathDesc.isValid()) {
        System.out.println(mathDesc.getVCML_database());
        throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
    }
    System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string begin ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
    System.out.println(mathDesc.getVCML());
    System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string end ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
}

Example 39 with SubDomain

use of cbit.vcell.math.SubDomain in project vcell by virtualcell.

the class FiniteVolumeFileWriter method writeModelDescription.

/**
 *# Model description: FEATURE name handle priority boundary_conditions
 *MODEL_BEGIN
 *FEATURE Nucleus 0 value value value value
 *FEATURE Cytosol 1 flux flux value value
 *FEATURE ExtraCellular 2 flux flux value value
 *MEMBRANE Nucleus_Cytosol_membrane Nucleus Cytosol value value value value
 *MEMBRANE Cytosol_ExtraCellular_membrane Cytosol ExtraCellular flux flux value value
 *MODEL_END
 * @throws MathException
 */
private void writeModelDescription() throws MathException {
    Simulation simulation = simTask.getSimulation();
    printWriter.println("# Model description: FEATURE name handle boundary_conditions");
    printWriter.println(FVInputFileKeyword.MODEL_BEGIN);
    MathDescription mathDesc = simulation.getMathDescription();
    Enumeration<SubDomain> enum1 = mathDesc.getSubDomains();
    while (enum1.hasMoreElements()) {
        SubDomain sd = enum1.nextElement();
        if (sd instanceof CompartmentSubDomain) {
            CompartmentSubDomain csd = (CompartmentSubDomain) sd;
            printWriter.print(FVInputFileKeyword.FEATURE + " " + csd.getName() + " " + mathDesc.getHandle(csd) + " ");
            writeFeature_boundaryConditions(csd);
        } else if (sd instanceof MembraneSubDomain) {
            MembraneSubDomain msd = (MembraneSubDomain) sd;
            printWriter.print(FVInputFileKeyword.MEMBRANE + " " + msd.getName() + " " + msd.getInsideCompartment().getName() + " " + msd.getOutsideCompartment().getName() + " ");
            writeMembrane_boundaryConditions(msd);
        }
    }
    printWriter.println(FVInputFileKeyword.MODEL_END);
    printWriter.println();
}

Example 40 with SubDomain

use of cbit.vcell.math.SubDomain in project vcell by virtualcell.

the class FiniteVolumeFileWriter method writeVariables.

/**
 *# Variables : type name unit time_dependent_flag advection_flag solve_whole_mesh_flag solve_regions
 *VARIABLE_BEGIN
 *VOLUME_ODE rB uM
 *VOLUME_PDE rf uM false false
 *VOLUME_PDE r uM false false
 *VOLUME_ODE rfB uM
 *VARIABLE_END
 * @throws MathException
 * @throws ExpressionException
 * @throws IOException
 */
private void writeVariables() throws MathException, ExpressionException, IOException {
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    printWriter.println("# Variables : type name domain time_dependent_flag advection_flag grad_flag solve_whole_mesh_flag solve_regions");
    printWriter.println(FVInputFileKeyword.VARIABLE_BEGIN);
    MathDescription mathDesc = simSymbolTable.getSimulation().getMathDescription();
    Variable[] vars = simSymbolTable.getVariables();
    ArrayList<RandomVariable> rvList = new ArrayList<RandomVariable>();
    for (int i = 0; i < vars.length; i++) {
        String varName = vars[i].getName();
        String domainName = vars[i].getDomain() == null ? null : vars[i].getDomain().getName();
        if (vars[i] instanceof VolumeRandomVariable || vars[i] instanceof MembraneRandomVariable) {
            rvList.add((RandomVariable) vars[i]);
        } else if (vars[i] instanceof VolVariable) {
            if (bChomboSolver && domainName == null) {
                throw new MathException(simTask.getSimulation().getSolverTaskDescription().getSolverDescription().getDisplayLabel() + " requires that every variable is defined in a single domain");
            }
            VolVariable volVar = (VolVariable) vars[i];
            if (mathDesc.isPDE(volVar)) {
                boolean hasTimeVaryingDiffusionOrAdvection = simSymbolTable.hasTimeVaryingDiffusionOrAdvection(volVar);
                final boolean hasVelocity = mathDesc.hasVelocity(volVar);
                final boolean hasGradient = mathDesc.hasGradient(volVar);
                if (mathDesc.isPdeSteady(volVar)) {
                    printWriter.print("VOLUME_PDE_STEADY ");
                } else {
                    printWriter.print("VOLUME_PDE ");
                }
                printWriter.print(varName + " " + domainName + " " + hasTimeVaryingDiffusionOrAdvection + " " + hasVelocity + " " + hasGradient);
            } else {
                printWriter.print("VOLUME_ODE " + varName + " " + domainName);
            }
            if (domainName == null) {
                Vector<SubDomain> listOfSubDomains = new Vector<SubDomain>();
                int totalNumCompartments = 0;
                Enumeration<SubDomain> subDomainEnum = mathDesc.getSubDomains();
                while (subDomainEnum.hasMoreElements()) {
                    SubDomain subDomain = subDomainEnum.nextElement();
                    if (subDomain instanceof CompartmentSubDomain) {
                        CompartmentSubDomain compartmentSubDomain = (CompartmentSubDomain) subDomain;
                        totalNumCompartments++;
                        Equation varEquation = subDomain.getEquation(vars[i]);
                        if (varEquation != null) {
                            if (!(varEquation instanceof PdeEquation) || !((PdeEquation) varEquation).isDummy(simSymbolTable, compartmentSubDomain)) {
                                listOfSubDomains.add(compartmentSubDomain);
                            }
                        }
                    }
                }
                if ((totalNumCompartments == listOfSubDomains.size()) || (listOfSubDomains.size() == 0 && simTask.getSimulation().getSolverTaskDescription().getSolverDescription().equals(SolverDescription.SundialsPDE))) {
                    printWriter.print(" true");
                } else {
                    printWriter.print(" false");
                    for (int j = 0; j < listOfSubDomains.size(); j++) {
                        CompartmentSubDomain compartmentSubDomain = (CompartmentSubDomain) listOfSubDomains.elementAt(j);
                        printWriter.print(" " + compartmentSubDomain.getName());
                    }
                }
                printWriter.println();
            } else {
                printWriter.println(" false " + domainName);
            }
        } else if (vars[i] instanceof VolumeParticleVariable) {
            printWriter.println(FVInputFileKeyword.VOLUME_PARTICLE + " " + varName + " " + domainName);
        } else if (vars[i] instanceof MembraneParticleVariable) {
            printWriter.println(FVInputFileKeyword.MEMBRANE_PARTICLE + " " + varName + " " + domainName);
        } else if (vars[i] instanceof VolumeRegionVariable) {
            printWriter.println("VOLUME_REGION " + varName + " " + domainName);
        } else if (vars[i] instanceof MemVariable) {
            if (bChomboSolver && domainName == null) {
                throw new MathException(simTask.getSimulation().getSolverTaskDescription().getSolverDescription().getDisplayLabel() + " requires that every variable is defined in a single domain");
            }
            MemVariable memVar = (MemVariable) vars[i];
            if (mathDesc.isPDE(memVar)) {
                printWriter.println("MEMBRANE_PDE " + varName + " " + domainName + " " + simSymbolTable.hasTimeVaryingDiffusionOrAdvection(memVar));
            } else {
                printWriter.println("MEMBRANE_ODE " + varName + " " + domainName);
            }
        } else if (vars[i] instanceof MembraneRegionVariable) {
            printWriter.println("MEMBRANE_REGION " + varName + " " + domainName);
        } else if (vars[i] instanceof FilamentVariable) {
            throw new RuntimeException("Filament application not supported yet");
        }
    }
    int numRandomVariables = rvList.size();
    if (numRandomVariables > 0) {
        ISize samplingSize = simTask.getSimulation().getMeshSpecification().getSamplingSize();
        String[] varNameArr = new String[numRandomVariables];
        VariableType[] varTypeArr = new VariableType[numRandomVariables];
        double[][] dataArr = new double[numRandomVariables][];
        for (int i = 0; i < numRandomVariables; i++) {
            RandomVariable rv = rvList.get(i);
            varNameArr[i] = rv.getName();
            int numRandomNumbers = 0;
            if (rv instanceof VolumeRandomVariable) {
                printWriter.print("VOLUME_RANDOM");
                varTypeArr[i] = VariableType.VOLUME;
                numRandomNumbers = samplingSize.getXYZ();
            } else if (rv instanceof MembraneRandomVariable) {
                printWriter.print("MEMBRANE_RANDOM");
                varTypeArr[i] = VariableType.MEMBRANE;
                numRandomNumbers = resampledGeometry.getGeometrySurfaceDescription().getSurfaceCollection().getTotalPolygonCount();
            } else {
                throw new RuntimeException("Unknown RandomVariable type");
            }
            printWriter.println(" " + varNameArr[i]);
            dataArr[i] = generateRandomNumbers(rv, numRandomNumbers);
        }
        File rvFile = new File(workingDirectory, simTask.getSimulationJobID() + RANDOM_VARIABLE_FILE_EXTENSION);
        DataSet.writeNew(rvFile, varNameArr, varTypeArr, samplingSize, dataArr);
    }
    printWriter.println(FVInputFileKeyword.VARIABLE_END);
    printWriter.println();
}