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

Example 26 with MathDescription

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

the class Simulation method getRequiredFeatures.

// public void gatherIssues(IssueContext issueContext, List<Issue> issueList) {
// 
// getMathOverrides().gatherIssues(issueContext, issueList);
// 
// //
// // Check if the math corresponding to this simulation has fast systems and if the solverTaskDescription contains a non-null sensitivity parameter.
// // If so, the simulation is invalid.
// //
// if (fieldMathDescription != null && getSolverTaskDescription() != null) {
// if (getMathDescription().hasFastSystems() && (getSolverTaskDescription().getSensitivityParameter() != null)) {
// Issue issue = new Issue(this, issueContext, IssueCategory.Simulation_SensAnal_And_FastSystem,
// VCellErrorMessages.getErrorMessage(VCellErrorMessages.SIMULATION_SENSANAL_FASTSYSTEM,getName()),
// Issue.SEVERITY_ERROR);
// issueList.add(issue);
// }
// }
// if (fieldMathDescription==null || !fieldMathDescription.isValid()){
// Issue issue = new Issue(this, issueContext, IssueCategory.MathDescription_MathException,fieldMathDescription.getWarning(),Issue.SEVERITY_ERROR);
// issueList.add(issue);
// }
// 
// Set<SolverFeature> supportedFeatures = getSolverTaskDescription().getSolverDescription().getSupportedFeatures();
// Set<SolverFeature> missingFeatures = getRequiredFeatures();
// missingFeatures.removeAll(supportedFeatures);
// 
// String text = "The selected Solver does not support the following required features: \n";
// for (SolverFeature sf : missingFeatures) {
// text += sf.getName() + "\n";
// }
// 
// if (!missingFeatures.isEmpty()) {
// System.out.println(this.getKey());
// String tooltip = "The selected Solver " + getSolverTaskDescription().getSolverDescription().getDisplayLabel() +
// " does not support the following required features: <br>";
// for (SolverFeature sf : missingFeatures) {
// tooltip += "&nbsp;&nbsp;&nbsp;" + sf.getName() + "<br>";
// }
// Collection<SolverDescription >goodSolvers = SolverDescription.getSolverDescriptions(getRequiredFeatures());
// assert goodSolvers != null;
// if (!goodSolvers.isEmpty()) {
// tooltip += "Please choose one of the solvers : <br>";
// for (SolverDescription sd : goodSolvers) {
// tooltip += "&nbsp;&nbsp;&nbsp;" + sd.getDisplayLabel() + "<br>";
// }
// }
// Issue issue = new Issue(this,issueContext, IssueCategory.MathDescription_MathException, text, tooltip, Issue.SEVERITY_ERROR);
// issueList.add(issue);
// }
// }
public Set<SolverFeature> getRequiredFeatures() {
    Set<SolverFeature> requiredFeatures = new HashSet<SolverFeature>();
    final MathDescription md = getMathDescription();
    if (isSpatial()) {
        requiredFeatures.add(SolverFeature.Feature_Spatial);
    } else {
        requiredFeatures.add(SolverFeature.Feature_NonSpatial);
    }
    final boolean hybrid = md.isSpatialHybrid();
    final boolean stoch = md.isNonSpatialStoch() || md.isSpatialStoch();
    final boolean ruleBased = md.isRuleBased();
    if (hybrid) {
        requiredFeatures.add(SolverFeature.Feature_Hybrid);
    }
    if (stoch && !hybrid) {
        requiredFeatures.add(SolverFeature.Feature_Stochastic);
    }
    if (!stoch && !ruleBased) {
        requiredFeatures.add(SolverFeature.Feature_Deterministic);
    }
    if (md.hasFastSystems()) {
        requiredFeatures.add(SolverFeature.Feature_FastSystem);
    }
    if (md.hasPeriodicBoundaryCondition()) {
        requiredFeatures.add(SolverFeature.Feature_PeriodicBoundaryCondition);
    }
    if (md.hasEvents()) {
        requiredFeatures.add(SolverFeature.Feature_Events);
    }
    if (md.hasRandomVariables()) {
        requiredFeatures.add(SolverFeature.Feature_RandomVariables);
    }
    if (getSolverTaskDescription().getStopAtSpatiallyUniformErrorTolerance() != null) {
        requiredFeatures.add(SolverFeature.Feature_StopAtSpatiallyUniform);
    }
    if (getDataProcessingInstructions() != null) {
        requiredFeatures.add(SolverFeature.Feature_DataProcessingInstructions);
    }
    if (md.getVariable(PSF_FUNCTION_NAME) != null) {
        requiredFeatures.add(SolverFeature.Feature_PSF);
    }
    if (isSerialParameterScan()) {
        requiredFeatures.add(SolverFeature.Feature_SerialParameterScans);
    }
    if (md.hasVolumeRegionEquations()) {
        requiredFeatures.add(SolverFeature.Feature_VolumeRegionEquations);
    }
    if (md.hasRegionSizeFunctions()) {
        requiredFeatures.add(SolverFeature.Feature_RegionSizeFunctions);
    }
    if (md.hasGradient()) {
        requiredFeatures.add(SolverFeature.Feature_GradientSourceTerm);
    }
    if (md.getPostProcessingBlock().getNumDataGenerators() > 0) {
        requiredFeatures.add(SolverFeature.Feature_PostProcessingBlock);
    }
    if (md.isRuleBased()) {
        requiredFeatures.add(SolverFeature.Feature_Rulebased);
    }
    return requiredFeatures;
}

Example 27 with MathDescription

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

the class SimulationWarning method analyzeDiffusion.

/**
 * make sure diffusion expressions are constants, store for later use
 * @throws ExpressionException
 */
private static Map<MembraneSubDomain, List<DiffusionValue>> analyzeDiffusion(Simulation simulation, double timeStep, IssueContext issueContext, List<Issue> issueList) throws ExpressionException {
    Map<MembraneSubDomain, List<DiffusionValue>> diffusionValuesMap = new IdentityHashMap<>();
    diffusionValuesMap.clear();
    MutableDouble value = new MutableDouble();
    MathDescription cm = simulation.getMathDescription();
    Objects.requireNonNull(cm);
    MathDescription localMath = new MathDescription(cm);
    SimulationSymbolTable symTable = new SimulationSymbolTable(simulation, 0);
    Map<MembraneSubDomain, List<DiffusionValue>> dvMap = new HashMap<>();
    double maxDiffValue = Double.MIN_VALUE;
    List<DiffusionValue> diffusionList = new ArrayList<>();
    for (SubDomain sd : localMath.getSubDomainCollection()) {
        final boolean isMembrane = sd instanceof MembraneSubDomain;
        diffusionList.clear();
        for (ParticleProperties pp : sd.getParticleProperties()) {
            String name = pp.getVariable().getName();
            Expression diffExp = pp.getDiffusion();
            Expression flattened = MathUtilities.substituteFunctions(diffExp, symTable).flatten();
            if (isConstant(flattened, value)) {
                if (isMembrane) {
                    DiffusionValue dv = new DiffusionValue(name, value.doubleValue());
                    maxDiffValue = Math.max(maxDiffValue, dv.value);
                    diffusionList.add(dv);
                }
            } else {
                String s = "Smoldyn only supports constant diffusion, " + name + " is variable";
                Issue i = new Issue(simulation, issueContext, IssueCategory.SMOLYDN_DIFFUSION, s, s, Severity.ERROR);
                issueList.add(i);
            }
        }
        if (isMembrane && !diffusionList.isEmpty()) {
            dvMap.put((MembraneSubDomain) sd, diffusionList);
        }
    }
    diffusionValuesMap.putAll(dvMap);
    MeshSpecification ms = simulation.getMeshSpecification();
    Geometry g = ms.getGeometry();
    int dim = g.getDimension();
    double minDelta = Double.MAX_VALUE;
    switch(dim) {
        case 3:
            minDelta = Math.min(minDelta, ms.getDz(true));
        // fall-through
        case 2:
            minDelta = Math.min(minDelta, ms.getDy(true));
        // fall-through
        case 1:
            minDelta = Math.min(minDelta, ms.getDx(true));
            break;
        default:
            throw new RuntimeException("Invalid dimension " + dim + " for smoldyn solver");
    }
    double minArea = minDelta * minDelta / 2;
    double limit = PRECHECK_LIMIT_ADJUST * minArea / maxDiffValue;
    boolean warn = (timeStep > limit);
    if (lg.isDebugEnabled()) {
        lg.debug("Min delta " + minDelta + ", min area " + minArea + " time limit " + limit + " timeStep " + timeStep + " -> warn = " + warn);
    }
    if (warn) {
        String s = "Time step " + timeStep + " may be too large, performing further analysis ...";
        Issue i = new Issue(simulation, issueContext, IssueCategory.SMOLYDN_DIFFUSION, s, s, Severity.WARNING);
        issueList.add(i);
    }
    lg.debug("end of diffusion analysis");
    return diffusionValuesMap;
}

Example 28 with MathDescription

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

the class IDAFileWriter method writeEquations.

/**
 * Insert the method's description here.
 * Creation date: (3/8/00 10:31:52 PM)
 */
protected String writeEquations(HashMap<Discontinuity, String> discontinuityNameMap) throws MathException, ExpressionException {
    Simulation simulation = simTask.getSimulation();
    StringBuffer sb = new StringBuffer();
    MathDescription mathDescription = simulation.getMathDescription();
    if (mathDescription.hasFastSystems()) {
        // 
        // define vector of original variables
        // 
        SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
        CompartmentSubDomain subDomain = (CompartmentSubDomain) mathDescription.getSubDomains().nextElement();
        FastSystem fastSystem = subDomain.getFastSystem();
        FastSystemAnalyzer fs_Analyzer = new FastSystemAnalyzer(fastSystem, simSymbolTable);
        int numIndependent = fs_Analyzer.getNumIndependentVariables();
        int systemDim = mathDescription.getStateVariableNames().size();
        int numDependent = systemDim - numIndependent;
        // 
        // get all variables from fast system (dependent and independent)
        // 
        HashSet<String> fastSystemVarHash = new HashSet<String>();
        Enumeration<Variable> dependentfastSystemVarEnum = fs_Analyzer.getDependentVariables();
        while (dependentfastSystemVarEnum.hasMoreElements()) {
            fastSystemVarHash.add(dependentfastSystemVarEnum.nextElement().getName());
        }
        Enumeration<Variable> independentfastSystemVarEnum = fs_Analyzer.getIndependentVariables();
        while (independentfastSystemVarEnum.hasMoreElements()) {
            fastSystemVarHash.add(independentfastSystemVarEnum.nextElement().getName());
        }
        // 
        // get all equations including for variables that are not in the fastSystem (ode equations for "slow system")
        // 
        RationalExpMatrix origInitVector = new RationalExpMatrix(systemDim, 1);
        RationalExpMatrix origSlowRateVector = new RationalExpMatrix(systemDim, 1);
        RationalExpMatrix origVarColumnVector = new RationalExpMatrix(systemDim, 1);
        Enumeration<Equation> enumEquations = subDomain.getEquations();
        int varIndex = 0;
        while (enumEquations.hasMoreElements()) {
            Equation equation = enumEquations.nextElement();
            origVarColumnVector.set_elem(varIndex, 0, new RationalExp(equation.getVariable().getName()));
            Expression rateExpr = equation.getRateExpression();
            rateExpr.bindExpression(varsSymbolTable);
            rateExpr = MathUtilities.substituteFunctions(rateExpr, varsSymbolTable);
            origSlowRateVector.set_elem(varIndex, 0, new RationalExp("(" + rateExpr.flatten().infix() + ")"));
            Expression initExpr = new Expression(equation.getInitialExpression());
            initExpr.substituteInPlace(new Expression("t"), new Expression(0.0));
            initExpr = MathUtilities.substituteFunctions(initExpr, varsSymbolTable).flatten();
            origInitVector.set_elem(varIndex, 0, new RationalExp("(" + initExpr.flatten().infix() + ")"));
            varIndex++;
        }
        // 
        // make symbolic matrix for fast invariants (from FastSystem's fastInvariants as well as a new fast invariant for each variable not included in the fast system.
        // 
        RationalExpMatrix fastInvarianceMatrix = new RationalExpMatrix(numDependent, systemDim);
        int row = 0;
        for (int i = 0; i < origVarColumnVector.getNumRows(); i++) {
            // 
            if (!fastSystemVarHash.contains(origVarColumnVector.get(i, 0).infixString())) {
                fastInvarianceMatrix.set_elem(row, i, RationalExp.ONE);
                row++;
            }
        }
        Enumeration<FastInvariant> enumFastInvariants = fastSystem.getFastInvariants();
        while (enumFastInvariants.hasMoreElements()) {
            FastInvariant fastInvariant = enumFastInvariants.nextElement();
            Expression fastInvariantExpression = fastInvariant.getFunction();
            for (int col = 0; col < systemDim; col++) {
                Expression coeff = fastInvariantExpression.differentiate(origVarColumnVector.get(col, 0).infixString()).flatten();
                coeff = simSymbolTable.substituteFunctions(coeff);
                fastInvarianceMatrix.set_elem(row, col, RationalExpUtils.getRationalExp(coeff));
            }
            row++;
        }
        for (int i = 0; i < systemDim; i++) {
            sb.append("VAR " + origVarColumnVector.get(i, 0).infixString() + " INIT " + origInitVector.get(i, 0).infixString() + ";\n");
        }
        RationalExpMatrix fullMatrix = null;
        RationalExpMatrix inverseFullMatrix = null;
        RationalExpMatrix newSlowRateVector = null;
        try {
            RationalExpMatrix fastMat = ((RationalExpMatrix) fastInvarianceMatrix.transpose().findNullSpace());
            fullMatrix = new RationalExpMatrix(systemDim, systemDim);
            row = 0;
            for (int i = 0; i < fastInvarianceMatrix.getNumRows(); i++) {
                for (int col = 0; col < systemDim; col++) {
                    fullMatrix.set_elem(row, col, fastInvarianceMatrix.get(i, col));
                }
                row++;
            }
            for (int i = 0; i < fastMat.getNumRows(); i++) {
                for (int col = 0; col < systemDim; col++) {
                    fullMatrix.set_elem(row, col, fastMat.get(i, col));
                }
                row++;
            }
            inverseFullMatrix = new RationalExpMatrix(systemDim, systemDim);
            inverseFullMatrix.identity();
            RationalExpMatrix copyOfFullMatrix = new RationalExpMatrix(fullMatrix);
            copyOfFullMatrix.gaussianElimination(inverseFullMatrix);
            newSlowRateVector = new RationalExpMatrix(numDependent, 1);
            newSlowRateVector.matmul(fastInvarianceMatrix, origSlowRateVector);
        } catch (MatrixException ex) {
            ex.printStackTrace();
            throw new MathException(ex.getMessage());
        }
        sb.append("TRANSFORM\n");
        for (row = 0; row < systemDim; row++) {
            for (int col = 0; col < systemDim; col++) {
                sb.append(fullMatrix.get(row, col).getConstant().doubleValue() + " ");
            }
            sb.append("\n");
        }
        sb.append("INVERSETRANSFORM\n");
        for (row = 0; row < systemDim; row++) {
            for (int col = 0; col < systemDim; col++) {
                sb.append(inverseFullMatrix.get(row, col).getConstant().doubleValue() + " ");
            }
            sb.append("\n");
        }
        int numDifferential = numDependent;
        int numAlgebraic = numIndependent;
        sb.append("RHS DIFFERENTIAL " + numDifferential + " ALGEBRAIC " + numAlgebraic + "\n");
        int equationIndex = 0;
        while (equationIndex < numDependent) {
            // print row of mass matrix followed by slow rate corresponding to fast invariant
            Expression slowRateExp = new Expression(newSlowRateVector.get(equationIndex, 0).infixString()).flatten();
            slowRateExp.bindExpression(simSymbolTable);
            slowRateExp = MathUtilities.substituteFunctions(slowRateExp, varsSymbolTable).flatten();
            Vector<Discontinuity> v = slowRateExp.getDiscontinuities();
            for (Discontinuity od : v) {
                od = getSubsitutedAndFlattened(od, varsSymbolTable);
                String dname = discontinuityNameMap.get(od);
                if (dname == null) {
                    dname = ROOT_VARIABLE_PREFIX + discontinuityNameMap.size();
                    discontinuityNameMap.put(od, dname);
                }
                slowRateExp.substituteInPlace(od.getDiscontinuityExp(), new Expression("(" + dname + "==1)"));
            }
            sb.append(slowRateExp.infix() + ";\n");
            equationIndex++;
        }
        Enumeration<FastRate> enumFastRates = fastSystem.getFastRates();
        while (enumFastRates.hasMoreElements()) {
            // print the fastRate for this row
            Expression fastRateExp = new Expression(enumFastRates.nextElement().getFunction());
            fastRateExp = MathUtilities.substituteFunctions(fastRateExp, varsSymbolTable).flatten();
            Vector<Discontinuity> v = fastRateExp.getDiscontinuities();
            for (Discontinuity od : v) {
                od = getSubsitutedAndFlattened(od, varsSymbolTable);
                String dname = discontinuityNameMap.get(od);
                if (dname == null) {
                    dname = ROOT_VARIABLE_PREFIX + discontinuityNameMap.size();
                    discontinuityNameMap.put(od, dname);
                }
                fastRateExp.substituteInPlace(od.getDiscontinuityExp(), new Expression("(" + dname + "==1)"));
            }
            sb.append(fastRateExp.flatten().infix() + ";\n");
            equationIndex++;
        }
    } else {
        for (int i = 0; i < getStateVariableCount(); i++) {
            StateVariable stateVar = getStateVariable(i);
            Expression initExpr = new Expression(stateVar.getInitialRateExpression());
            initExpr = MathUtilities.substituteFunctions(initExpr, varsSymbolTable);
            initExpr.substituteInPlace(new Expression("t"), new Expression(0.0));
            sb.append("VAR " + stateVar.getVariable().getName() + " INIT " + initExpr.flatten().infix() + ";\n");
        }
        sb.append("TRANSFORM\n");
        for (int row = 0; row < getStateVariableCount(); row++) {
            for (int col = 0; col < getStateVariableCount(); col++) {
                sb.append((row == col ? 1 : 0) + " ");
            }
            sb.append("\n");
        }
        sb.append("INVERSETRANSFORM\n");
        for (int row = 0; row < getStateVariableCount(); row++) {
            for (int col = 0; col < getStateVariableCount(); col++) {
                sb.append((row == col ? 1 : 0) + " ");
            }
            sb.append("\n");
        }
        sb.append("RHS DIFFERENTIAL " + getStateVariableCount() + " ALGEBRAIC 0\n");
        for (int i = 0; i < getStateVariableCount(); i++) {
            StateVariable stateVar = getStateVariable(i);
            Expression rateExpr = new Expression(stateVar.getRateExpression());
            rateExpr = MathUtilities.substituteFunctions(rateExpr, varsSymbolTable).flatten();
            Vector<Discontinuity> v = rateExpr.getDiscontinuities();
            for (Discontinuity od : v) {
                od = getSubsitutedAndFlattened(od, varsSymbolTable);
                String dname = discontinuityNameMap.get(od);
                if (dname == null) {
                    dname = ROOT_VARIABLE_PREFIX + discontinuityNameMap.size();
                    discontinuityNameMap.put(od, dname);
                }
                rateExpr.substituteInPlace(od.getDiscontinuityExp(), new Expression("(" + dname + "==1)"));
            }
            sb.append(rateExpr.infix() + ";\n");
        }
    }
    return sb.toString();
}

Example 29 with MathDescription

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

the class SundialsSolver method createStateVariables.

/*
	This method was created in Visual Age
*/
private StateVariable[] createStateVariables() throws MathException, ExpressionException {
    Vector<StateVariable> stateVariables = new Vector<StateVariable>();
    // get Ode's from MathDescription and create ODEStateVariables
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    MathDescription mathDescription = simSymbolTable.getSimulation().getMathDescription();
    Enumeration<Equation> enum1 = ((SubDomain) mathDescription.getSubDomains().nextElement()).getEquations();
    while (enum1.hasMoreElements()) {
        Equation equation = (Equation) enum1.nextElement();
        if (equation instanceof OdeEquation) {
            stateVariables.addElement(new ODEStateVariable((OdeEquation) equation, simSymbolTable));
        } else {
            throw new MathException("encountered non-ode equation, unsupported");
        }
    }
    // Get sensitivity variables
    Variable[] variables = simSymbolTable.getVariables();
    Vector<SensVariable> sensVariables = new Vector<SensVariable>();
    if (getSensitivityParameter() != null) {
        for (int i = 0; i < variables.length; i++) {
            if (variables[i] instanceof VolVariable) {
                VolVariable volVariable = (VolVariable) variables[i];
                SensVariable sv = new SensVariable(volVariable, getSensitivityParameter());
                sensVariables.addElement(sv);
            }
        }
    }
    if (rateSensitivity == null) {
        rateSensitivity = new RateSensitivity(mathDescription, mathDescription.getSubDomains().nextElement());
    }
    if (jacobian == null) {
        jacobian = new Jacobian(mathDescription, mathDescription.getSubDomains().nextElement());
    }
    for (int v = 0; v < sensVariables.size(); v++) {
        stateVariables.addElement(new SensStateVariable((SensVariable) sensVariables.elementAt(v), rateSensitivity, jacobian, sensVariables, simSymbolTable));
    }
    if (stateVariables.size() == 0) {
        throw new MathException("there are no equations defined");
    }
    StateVariable[] stateVars = (StateVariable[]) BeanUtils.getArray(stateVariables, StateVariable.class);
    return (stateVars);
}

Example 30 with MathDescription

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

the class OdeFileWriter method createStateVariables.

private void createStateVariables() throws Exception {
    Simulation simulation = simTask.getSimulation();
    MathDescription mathDescription = simulation.getMathDescription();
    SolverTaskDescription solverTaskDescription = simulation.getSolverTaskDescription();
    // get Ode's from MathDescription and create ODEStateVariables
    Enumeration<Equation> enum1 = mathDescription.getSubDomains().nextElement().getEquations();
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    while (enum1.hasMoreElements()) {
        Equation equation = enum1.nextElement();
        if (equation instanceof OdeEquation) {
            fieldStateVariables.addElement(new ODEStateVariable((OdeEquation) equation, simSymbolTable));
        } else {
            throw new MathException("encountered non-ode equation, unsupported");
        }
    }
    // Get sensitivity variables
    Variable[] variables = simSymbolTable.getVariables();
    Vector<SensVariable> sensVariables = new Vector<SensVariable>();
    Constant sensitivityParameter = solverTaskDescription.getSensitivityParameter();
    if (sensitivityParameter != null) {
        Constant origSensParam = sensitivityParameter;
        Constant overriddenSensParam = (Constant) simSymbolTable.getVariable(origSensParam.getName());
        for (int i = 0; i < variables.length; i++) {
            if (variables[i] instanceof VolVariable) {
                VolVariable volVariable = (VolVariable) variables[i];
                SensVariable sv = new SensVariable(volVariable, overriddenSensParam);
                sensVariables.addElement(sv);
            }
        }
    }
    if (rateSensitivity == null) {
        rateSensitivity = new RateSensitivity(mathDescription, mathDescription.getSubDomains().nextElement());
    }
    if (jacobian == null) {
        jacobian = new Jacobian(mathDescription, mathDescription.getSubDomains().nextElement());
    }
    // get Jacobian and RateSensitivities from MathDescription and create SensStateVariables
    for (int v = 0; v < sensVariables.size(); v++) {
        fieldStateVariables.addElement(new SensStateVariable(sensVariables.elementAt(v), rateSensitivity, jacobian, sensVariables, simSymbolTable));
    }
}