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

Example 41 with Variable

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

the class MathTestingUtilities method constructOdesForSensitivity.

// 
// This method is used to solve for sensitivity of variables to a given parameter.
// The mathDescription and the sensitivity parameter are passed as arguments.
// New variables and ODEs are constructed according to the rule listed below and are added to the mathDescription.
// The method returns the modified mathDescription.
// 
public static MathDescription constructOdesForSensitivity(MathDescription mathDesc, Constant sensParam) throws ExpressionException, MathException, MappingException {
    // 
    // For each ODE :
    // 
    // dX/dt = F(X, P)
    // 
    // (where P is the sensitivity parameter)
    // 
    // we create two other ODEs :
    // 
    // dX_1/dt = F(X_1, P_1)    and
    // 
    // dX_2/dt = F(X_2, P_2)
    // 
    // where P_1 = P + epsilon, and
    // P_2 = P - epsilon.
    // 
    // We keep the initial conditions for both the new ODEs to be the same,
    // i.e., X_1_init = X_2_init.
    // 
    // Then, solving for X_1 & X_2, sensitivity of X wrt P can be computed as :
    // 
    // dX = (X_1 - X_2)
    // --	 -----------   .
    // dP	 (P_1 - P_2)
    // 
    // 
    // REMOVE PRINTS AFTER CHECKING !!!
    System.out.println(" \n\n------------  Old Math Description -----------------");
    System.out.println(mathDesc.getVCML_database());
    if (mathDesc.getGeometry().getDimension() > 0) {
        throw new RuntimeException("Suppport for Spatial systems not yet implemented.");
    }
    VariableHash varHash = new VariableHash();
    Enumeration<Variable> enumVar = mathDesc.getVariables();
    while (enumVar.hasMoreElements()) {
        varHash.addVariable(enumVar.nextElement());
    }
    // 
    // Get 2 values of senstivity parameter (P + epsilon) & (P - epsilon)
    // 
    Constant epsilon = new Constant("epsilon", new Expression(sensParam.getConstantValue() * 1e-3));
    Constant sensParam1 = new Constant(sensParam.getName() + "_1", new Expression(sensParam.getConstantValue() + epsilon.getConstantValue()));
    Constant sensParam2 = new Constant(sensParam.getName() + "_2", new Expression(sensParam.getConstantValue() - epsilon.getConstantValue()));
    // 
    // Iterate through each subdomain (only 1 in compartmental case), and each equation in the subdomain
    // 
    Enumeration<SubDomain> subDomainEnum = mathDesc.getSubDomains();
    // 
    // Create a vector of equations to store the 2 equations for each ODE variable in the subdomain.
    // Later, add it to the equations list in the subdomain.
    // 
    Vector<Equation> equnsVector = new Vector<Equation>();
    Vector<Variable> varsVector = new Vector<Variable>();
    Vector<Variable> var1s = new Vector<Variable>();
    Vector<Variable> var2s = new Vector<Variable>();
    while (subDomainEnum.hasMoreElements()) {
        SubDomain subDomain = subDomainEnum.nextElement();
        Enumeration<Equation> equationEnum = subDomain.getEquations();
        Domain domain = new Domain(subDomain);
        while (equationEnum.hasMoreElements()) {
            Equation equation = equationEnum.nextElement();
            if (equation instanceof OdeEquation) {
                OdeEquation odeEquation = (OdeEquation) equation;
                // Similar to substituteWithExactSolutions, to bind and substitute functions in the ODE
                Expression substitutedRateExp = substituteFunctions(odeEquation.getRateExpression(), mathDesc);
                String varName = odeEquation.getVariable().getName();
                VolVariable var = new VolVariable(varName, domain);
                varsVector.addElement(var);
                // 
                // Create the variable var1, and get the initExpr and rateExpr from the original ODE.
                // Substitute the new vars (var1 and param1) in the old initExpr and rateExpr and create a new ODE
                // 
                String varName1 = new String("__" + varName + "_1");
                Expression initExpr1 = odeEquation.getInitialExpression();
                Expression rateExpr1 = new Expression(substitutedRateExp);
                rateExpr1.substituteInPlace(new Expression(varName), new Expression(varName1));
                rateExpr1.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam1.getName()));
                VolVariable var1 = new VolVariable(varName1, domain);
                var1s.addElement(var1);
                OdeEquation odeEqun1 = new OdeEquation(var1, initExpr1, rateExpr1);
                equnsVector.addElement(odeEqun1);
                // 
                // Create the variable var2, and get the initExpr and rateExpr from the original ODE.
                // Substitute the new vars (var2 and param2) in the old initExpr and rateExpr and create a new ODE
                // 
                String varName2 = new String("__" + varName + "_2");
                Expression initExpr2 = odeEquation.getInitialExpression();
                Expression rateExpr2 = new Expression(substitutedRateExp);
                rateExpr2.substituteInPlace(new Expression(varName), new Expression(varName2));
                rateExpr2.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam2.getName()));
                VolVariable var2 = new VolVariable(varName2, domain);
                var2s.addElement(var2);
                OdeEquation odeEqun2 = new OdeEquation(var2, initExpr2, rateExpr2);
                equnsVector.addElement(odeEqun2);
                // 
                // Create a function for the sensitivity function expression (X1-X2)/(P1-P2), and save in varHash
                // 
                Expression diffVar = Expression.add(new Expression(var1.getName()), Expression.negate(new Expression(var2.getName())));
                Expression diffParam = Expression.add(new Expression(sensParam1.getName()), Expression.negate(new Expression(sensParam2.getName())));
                Expression sensitivityExpr = Expression.mult(diffVar, Expression.invert(diffParam));
                Function sens_Func = new Function("__sens" + varName + "_wrt_" + sensParam.getName(), sensitivityExpr, domain);
                varHash.addVariable(epsilon);
                varHash.addVariable(sensParam1);
                varHash.addVariable(sensParam2);
                varHash.addVariable(var1);
                varHash.addVariable(var2);
                varHash.addVariable(sens_Func);
            } else {
                // sensitivity not implemented for PDEs or other equation types.
                throw new RuntimeException("SolverTest.constructedExactMath(): equation type " + equation.getClass().getName() + " not yet implemented");
            }
        }
        // 
        // Need to substitute the new variables in the new ODEs.
        // i.e., if Rate Expr for ODE_1 for variable X_1 contains variable Y, variable Z, etc.
        // Rate Expr is already substituted with X_1, but it also needs substitute Y with Y_1, Z with Z_1, etc.
        // So get the volume variables, from the vectors for vars, var_1s and var_2s
        // Substitute the rate expressions for the newly added ODEs in equnsVector.
        // 
        Variable[] vars = (Variable[]) BeanUtils.getArray(varsVector, Variable.class);
        Variable[] var_1s = (Variable[]) BeanUtils.getArray(var1s, Variable.class);
        Variable[] var_2s = (Variable[]) BeanUtils.getArray(var2s, Variable.class);
        Vector<Equation> newEqunsVector = new Vector<Equation>();
        for (int i = 0; i < equnsVector.size(); i++) {
            Equation equn = equnsVector.elementAt(i);
            Expression initEx = equn.getInitialExpression();
            Expression rateEx = equn.getRateExpression();
            for (int j = 0; j < vars.length; j++) {
                if (equn.getVariable().getName().endsWith("_1")) {
                    rateEx.substituteInPlace(new Expression(vars[j].getName()), new Expression(var_1s[j].getName()));
                } else if (equn.getVariable().getName().endsWith("_2")) {
                    rateEx.substituteInPlace(new Expression(vars[j].getName()), new Expression(var_2s[j].getName()));
                }
            }
            OdeEquation odeEqun = new OdeEquation(equn.getVariable(), initEx, rateEx);
            newEqunsVector.addElement(odeEqun);
        }
        // 
        for (int i = 0; i < newEqunsVector.size(); i++) {
            mathDesc.getSubDomain(subDomain.getName()).addEquation((Equation) newEqunsVector.elementAt(i));
        }
        // 
        // FAST SYSTEM
        // If the subdomain has a fast system, create a new fast system by substituting the high-low variables/parameters
        // in the expressions for the fastInvariants and fastRates and adding them to the fast system.
        // 
        Vector<FastInvariant> invarsVector = new Vector<FastInvariant>();
        Vector<FastRate> ratesVector = new Vector<FastRate>();
        Enumeration<FastInvariant> fastInvarsEnum = null;
        Enumeration<FastRate> fastRatesEnum = null;
        // Get the fast invariants and fast rates in the system.
        FastSystem fastSystem = subDomain.getFastSystem();
        if (fastSystem != null) {
            fastInvarsEnum = fastSystem.getFastInvariants();
            fastRatesEnum = fastSystem.getFastRates();
            // 
            while (fastInvarsEnum.hasMoreElements()) {
                FastInvariant fastInvar = fastInvarsEnum.nextElement();
                Expression fastInvarExpr = fastInvar.getFunction();
                fastInvarExpr = MathUtilities.substituteFunctions(fastInvarExpr, mathDesc);
                Expression fastInvarExpr1 = new Expression(fastInvarExpr);
                Expression fastInvarExpr2 = new Expression(fastInvarExpr);
                for (int i = 0; i < vars.length; i++) {
                    fastInvarExpr1.substituteInPlace(new Expression(vars[i].getName()), new Expression(var_1s[i].getName()));
                    fastInvarExpr2.substituteInPlace(new Expression(vars[i].getName()), new Expression(var_2s[i].getName()));
                }
                fastInvarExpr1.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam1.getName()));
                FastInvariant fastInvar1 = new FastInvariant(fastInvarExpr1);
                invarsVector.addElement(fastInvar1);
                fastInvarExpr2.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam2.getName()));
                FastInvariant fastInvar2 = new FastInvariant(fastInvarExpr2);
                invarsVector.addElement(fastInvar2);
            }
            // Add the newly created fast invariants to the existing list of fast invariants in the fast system.
            for (int i = 0; i < invarsVector.size(); i++) {
                FastInvariant inVar = (FastInvariant) invarsVector.elementAt(i);
                fastSystem.addFastInvariant(inVar);
            }
            // 
            while (fastRatesEnum.hasMoreElements()) {
                FastRate fastRate = fastRatesEnum.nextElement();
                Expression fastRateExpr = fastRate.getFunction();
                fastRateExpr = MathUtilities.substituteFunctions(fastRateExpr, mathDesc);
                Expression fastRateExpr1 = new Expression(fastRateExpr);
                Expression fastRateExpr2 = new Expression(fastRateExpr);
                for (int i = 0; i < vars.length; i++) {
                    fastRateExpr1.substituteInPlace(new Expression(vars[i].getName()), new Expression(var_1s[i].getName()));
                    fastRateExpr2.substituteInPlace(new Expression(vars[i].getName()), new Expression(var_2s[i].getName()));
                }
                fastRateExpr1.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam1.getName()));
                FastRate fastRate1 = new FastRate(fastRateExpr1);
                ratesVector.addElement(fastRate1);
                fastRateExpr2.substituteInPlace(new Expression(sensParam.getName()), new Expression(sensParam2.getName()));
                FastRate fastRate2 = new FastRate(fastRateExpr2);
                ratesVector.addElement(fastRate2);
            }
            // Add the newly created fast rates to the existing list of fast rates in the fast system.
            for (int i = 0; i < ratesVector.size(); i++) {
                FastRate rate = (FastRate) ratesVector.elementAt(i);
                fastSystem.addFastRate(rate);
            }
        }
    }
    // Reset all variables in mathDesc.
    mathDesc.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
    // REMOVE PRINTS AFTER CHECKING
    System.out.println(" \n\n------------  New Math Description -----------------");
    System.out.println(mathDesc.getVCML_database());
    return mathDesc;
}

Example 42 with Variable

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

the class MathTestingUtilities method getOutwardNormal.

/**
 * Insert the method's description here.
 * Creation date: (1/23/2003 10:30:23 PM)
 * @return cbit.vcell.parser.Expression
 * @param analyticSubDomainExp cbit.vcell.parser.Expression
 */
public static Function[] getOutwardNormal(Expression analyticSubVolume, String baseName) throws ExpressionException, MappingException, MathException {
    VariableHash varHash = new VariableHash();
    Expression[] insideOutsideFunctions = getInsideOutsideFunctions(analyticSubVolume);
    StringBuffer normalBufferX = new StringBuffer("0.0");
    StringBuffer normalBufferY = new StringBuffer("0.0");
    StringBuffer normalBufferZ = new StringBuffer("0.0");
    for (int i = 0; i < insideOutsideFunctions.length; i++) {
        // 
        // each one gets a turn being the minimum
        // 
        Function[] functions = getOutwardNormalFromInsideOutsideFunction(insideOutsideFunctions[i], baseName + i);
        for (int j = 0; j < functions.length; j++) {
            varHash.addVariable(functions[j]);
        }
        String closestName = baseName + i + "_closest";
        StringBuffer closestBuffer = new StringBuffer("1.0");
        for (int j = 0; j < insideOutsideFunctions.length; j++) {
            if (i != j) {
                closestBuffer.append(" && (" + baseName + i + "_distance < " + baseName + j + "_distance)");
            }
        }
        Expression closest = new Expression(closestBuffer.toString());
        varHash.addVariable(new Function(closestName, closest, null));
        normalBufferX.append(" + (" + baseName + i + "_closest * " + baseName + i + "_Nx)");
        normalBufferY.append(" + (" + baseName + i + "_closest * " + baseName + i + "_Ny)");
        normalBufferZ.append(" + (" + baseName + i + "_closest * " + baseName + i + "_Nz)");
    }
    varHash.addVariable(new Function(baseName + "_Nx", new Expression(normalBufferX.toString()), null));
    varHash.addVariable(new Function(baseName + "_Ny", new Expression(normalBufferY.toString()), null));
    varHash.addVariable(new Function(baseName + "_Nz", new Expression(normalBufferZ.toString()), null));
    Variable[] vars = varHash.getAlphabeticallyOrderedVariables();
    java.util.Vector<Variable> varList = new java.util.Vector<Variable>(java.util.Arrays.asList(vars));
    return (Function[]) BeanUtils.getArray(varList, Function.class);
}

Example 43 with Variable

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

the class MathTestingUtilities method constructExactMath.

/**
 * constructExactMath()
 *
 * take an equation of the form:
 *
 *     d A
 *     --- = F(A,t)
 *     d t
 *
 * and create a new equation with a known exact solution 'A_exact' by adding a forcing function R(t)
 *
 *     d A
 *     --- = F(A,t) + R(t)
 *     d t
 *
 * where:
 *
 *             d A_exact
 *      R(t) = --------- - F(A_exact,t)
 *                d t
 *
 * solving for R(t) is done analytically.
 *
 * Creation date: (1/21/2003 10:47:54 AM)
 * @return cbit.vcell.math.MathDescription
 * @param mathDesc cbit.vcell.math.MathDescription
 */
public static MathDescription constructExactMath(MathDescription mathDesc, java.util.Random random, ConstructedSolutionTemplate constructedSolutionTemplate) throws ExpressionException, MathException, MappingException {
    if (mathDesc.hasFastSystems()) {
        throw new RuntimeException("SolverTest.constructExactMath() suppport for fastSystems not yet implemented.");
    }
    MathDescription exactMath = null;
    try {
        exactMath = (MathDescription) BeanUtils.cloneSerializable(mathDesc);
        exactMath.setDescription("constructed exact solution from MathDescription (" + mathDesc.getName() + ")");
        exactMath.setName("exact from " + mathDesc.getName());
    } catch (Throwable e) {
        e.printStackTrace(System.out);
        throw new RuntimeException("error cloning MathDescription: " + e.getMessage());
    }
    // 
    // preload the VariableHash with existing Variables (and Constants,Functions,etc) and then sort all at once.
    // 
    VariableHash varHash = new VariableHash();
    Enumeration<Variable> enumVar = exactMath.getVariables();
    while (enumVar.hasMoreElements()) {
        varHash.addVariable(enumVar.nextElement());
    }
    java.util.Enumeration<SubDomain> subDomainEnum = exactMath.getSubDomains();
    while (subDomainEnum.hasMoreElements()) {
        SubDomain subDomain = subDomainEnum.nextElement();
        Domain domain = new Domain(subDomain);
        java.util.Enumeration<Equation> equationEnum = subDomain.getEquations();
        if (subDomain instanceof MembraneSubDomain) {
            MembraneSubDomain memSubDomain = (MembraneSubDomain) subDomain;
            AnalyticSubVolume insideAnalyticSubVolume = (AnalyticSubVolume) exactMath.getGeometry().getGeometrySpec().getSubVolume(memSubDomain.getInsideCompartment().getName());
            Function[] outwardNormalFunctions = getOutwardNormal(insideAnalyticSubVolume.getExpression(), "_" + insideAnalyticSubVolume.getName());
            for (int i = 0; i < outwardNormalFunctions.length; i++) {
                varHash.addVariable(outwardNormalFunctions[i]);
            }
        }
        while (equationEnum.hasMoreElements()) {
            Equation equation = equationEnum.nextElement();
            if (equation.getExactSolution() != null) {
                throw new RuntimeException("exact solution already exists");
            }
            Enumeration<Constant> origMathConstants = mathDesc.getConstants();
            if (equation instanceof OdeEquation) {
                OdeEquation odeEquation = (OdeEquation) equation;
                Expression substitutedRateExp = substituteWithExactSolution(odeEquation.getRateExpression(), (CompartmentSubDomain) subDomain, exactMath);
                SolutionTemplate solutionTemplate = constructedSolutionTemplate.getSolutionTemplate(equation.getVariable().getName(), subDomain.getName());
                String varName = odeEquation.getVariable().getName();
                String initName = null;
                while (origMathConstants.hasMoreElements()) {
                    Constant constant = origMathConstants.nextElement();
                    if (constant.getName().startsWith(varName + "_" + subDomain.getName() + DiffEquMathMapping.MATH_FUNC_SUFFIX_SPECIES_INIT_CONC_UNIT_PREFIX)) {
                        initName = constant.getName();
                    }
                }
                String exactName = varName + "_" + subDomain.getName() + "_exact";
                String errorName = varName + "_" + subDomain.getName() + "_error";
                String origRateName = "_" + varName + "_" + subDomain.getName() + "_origRate";
                String substitutedRateName = "_" + varName + "_" + subDomain.getName() + "_substitutedRate";
                String exactTimeDerivativeName = "_" + varName + "_" + subDomain.getName() + "_exact_dt";
                Expression exactExp = solutionTemplate.getTemplateExpression();
                Expression errorExp = new Expression(exactName + " - " + varName);
                Expression origRateExp = new Expression(odeEquation.getRateExpression());
                Expression exactTimeDerivativeExp = exactExp.differentiate("t").flatten();
                Expression newRate = new Expression(origRateName + " - " + substitutedRateName + " + " + exactTimeDerivativeName);
                Constant[] constants = solutionTemplate.getConstants();
                for (int i = 0; i < constants.length; i++) {
                    varHash.addVariable(constants[i]);
                }
                Expression initExp = new Expression(exactExp);
                initExp.substituteInPlace(new Expression("t"), new Expression(0.0));
                varHash.addVariable(new Function(initName, initExp.flatten(), domain));
                varHash.addVariable(new Function(exactName, exactExp, domain));
                varHash.addVariable(new Function(errorName, errorExp, domain));
                varHash.addVariable(new Function(exactTimeDerivativeName, exactTimeDerivativeExp, domain));
                varHash.addVariable(new Function(origRateName, origRateExp, domain));
                varHash.addVariable(new Function(substitutedRateName, substitutedRateExp, domain));
                odeEquation.setRateExpression(newRate);
                odeEquation.setInitialExpression(new Expression(initName));
                odeEquation.setExactSolution(new Expression(exactName));
            } else if (equation instanceof PdeEquation) {
                PdeEquation pdeEquation = (PdeEquation) equation;
                Expression substitutedRateExp = substituteWithExactSolution(pdeEquation.getRateExpression(), (CompartmentSubDomain) subDomain, exactMath);
                SolutionTemplate solutionTemplate = constructedSolutionTemplate.getSolutionTemplate(equation.getVariable().getName(), subDomain.getName());
                String varName = pdeEquation.getVariable().getName();
                String initName = null;
                while (origMathConstants.hasMoreElements()) {
                    Constant constant = origMathConstants.nextElement();
                    if (constant.getName().startsWith(varName + "_" + subDomain.getName() + DiffEquMathMapping.MATH_FUNC_SUFFIX_SPECIES_INIT_CONC_UNIT_PREFIX)) {
                        initName = constant.getName();
                    }
                }
                String diffusionRateName = "_" + varName + "_" + subDomain.getName() + "_diffusionRate";
                String exactName = varName + "_" + subDomain.getName() + "_exact";
                String errorName = varName + "_" + subDomain.getName() + "_error";
                String origRateName = "_" + varName + "_" + subDomain.getName() + "_origRate";
                String substitutedRateName = "_" + varName + "_" + subDomain.getName() + "_substitutedRate";
                String exactTimeDerivativeName = "_" + varName + "_" + subDomain.getName() + "_exact_dt";
                String exactDxName = "_" + varName + "_" + subDomain.getName() + "_exact_dx";
                String exactDyName = "_" + varName + "_" + subDomain.getName() + "_exact_dy";
                String exactDzName = "_" + varName + "_" + subDomain.getName() + "_exact_dz";
                String exactDx2Name = "_" + varName + "_" + subDomain.getName() + "_exact_dx2";
                String exactDy2Name = "_" + varName + "_" + subDomain.getName() + "_exact_dy2";
                String exactDz2Name = "_" + varName + "_" + subDomain.getName() + "_exact_dz2";
                String exactLaplacianName = "_" + varName + "_" + subDomain.getName() + "_exact_laplacian";
                Expression exactExp = solutionTemplate.getTemplateExpression();
                Expression errorExp = new Expression(exactName + " - " + varName);
                Expression origRateExp = new Expression(pdeEquation.getRateExpression());
                Expression initExp = new Expression(exactExp);
                initExp.substituteInPlace(new Expression("t"), new Expression(0.0));
                initExp = initExp.flatten();
                Expression exactTimeDerivativeExp = exactExp.differentiate("t").flatten();
                Expression exactDxExp = exactExp.differentiate("x").flatten();
                Expression exactDx2Exp = exactDxExp.differentiate("x").flatten();
                Expression exactDyExp = exactExp.differentiate("y").flatten();
                Expression exactDy2Exp = exactDxExp.differentiate("y").flatten();
                Expression exactDzExp = exactExp.differentiate("z").flatten();
                Expression exactDz2Exp = exactDxExp.differentiate("z").flatten();
                Expression exactLaplacianExp = Expression.add(Expression.add(exactDx2Exp, exactDy2Exp), exactDz2Exp).flatten();
                Expression newRate = new Expression(origRateName + " - " + substitutedRateName + " - ((" + diffusionRateName + ")*" + exactLaplacianName + ")" + " + " + exactTimeDerivativeName);
                Constant[] constants = solutionTemplate.getConstants();
                for (int i = 0; i < constants.length; i++) {
                    varHash.addVariable(constants[i]);
                }
                varHash.addVariable(new Function(initName, initExp, domain));
                varHash.addVariable(new Function(diffusionRateName, new Expression(pdeEquation.getDiffusionExpression()), domain));
                varHash.addVariable(new Function(exactName, exactExp, domain));
                varHash.addVariable(new Function(errorName, errorExp, domain));
                varHash.addVariable(new Function(exactTimeDerivativeName, exactTimeDerivativeExp, domain));
                varHash.addVariable(new Function(origRateName, origRateExp, domain));
                varHash.addVariable(new Function(substitutedRateName, substitutedRateExp, domain));
                varHash.addVariable(new Function(exactDxName, exactDxExp, domain));
                varHash.addVariable(new Function(exactDyName, exactDyExp, domain));
                varHash.addVariable(new Function(exactDzName, exactDzExp, domain));
                varHash.addVariable(new Function(exactDx2Name, exactDx2Exp, domain));
                varHash.addVariable(new Function(exactDy2Name, exactDy2Exp, domain));
                varHash.addVariable(new Function(exactDz2Name, exactDz2Exp, domain));
                varHash.addVariable(new Function(exactLaplacianName, exactLaplacianExp, domain));
                pdeEquation.setRateExpression(newRate);
                pdeEquation.setInitialExpression(new Expression(initName));
                pdeEquation.setDiffusionExpression(new Expression(diffusionRateName));
                pdeEquation.setExactSolution(new Expression(exactName));
                CompartmentSubDomain compartmentSubDomain = (CompartmentSubDomain) subDomain;
                if (compartmentSubDomain.getBoundaryConditionXm().isDIRICHLET()) {
                    Expression origExp = pdeEquation.getBoundaryXm();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryXm(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
                    } else {
                        pdeEquation.setBoundaryXm(exactExp);
                    }
                } else if (compartmentSubDomain.getBoundaryConditionXm().isNEUMANN()) {
                    Expression origExp = pdeEquation.getBoundaryXm();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryXm(new Expression(origExp + "-" + substitutedExp + "-" + diffusionRateName + "*" + exactDxName));
                    } else {
                        pdeEquation.setBoundaryXm(new Expression("-" + diffusionRateName + "*" + exactDxName));
                    }
                } else {
                    throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionXm());
                }
                if (compartmentSubDomain.getBoundaryConditionXp().isDIRICHLET()) {
                    Expression origExp = pdeEquation.getBoundaryXp();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryXp(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
                    } else {
                        pdeEquation.setBoundaryXp(exactExp);
                    }
                } else if (compartmentSubDomain.getBoundaryConditionXp().isNEUMANN()) {
                    Expression origExp = pdeEquation.getBoundaryXp();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryXp(new Expression(origExp + "-" + substitutedExp + "+" + diffusionRateName + "*" + exactDxName));
                    } else {
                        pdeEquation.setBoundaryXp(new Expression(diffusionRateName + "*" + exactDxName));
                    }
                } else {
                    throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionXp());
                }
                if (compartmentSubDomain.getBoundaryConditionYm().isDIRICHLET()) {
                    Expression origExp = pdeEquation.getBoundaryYm();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryYm(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
                    } else {
                        pdeEquation.setBoundaryYm(exactExp);
                    }
                } else if (compartmentSubDomain.getBoundaryConditionYm().isNEUMANN()) {
                    Expression origExp = pdeEquation.getBoundaryYm();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryYm(new Expression(origExp + "-" + substitutedExp + "-" + diffusionRateName + "*" + exactDyName));
                    } else {
                        pdeEquation.setBoundaryYm(new Expression("-" + diffusionRateName + "*" + exactDyName));
                    }
                } else {
                    throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionYm());
                }
                if (compartmentSubDomain.getBoundaryConditionYp().isDIRICHLET()) {
                    Expression origExp = pdeEquation.getBoundaryYp();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryYp(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
                    } else {
                        pdeEquation.setBoundaryYp(exactExp);
                    }
                } else if (compartmentSubDomain.getBoundaryConditionYp().isNEUMANN()) {
                    Expression origExp = pdeEquation.getBoundaryYp();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryYp(new Expression(origExp + "-" + substitutedExp + "+" + diffusionRateName + "*" + exactDyName));
                    } else {
                        pdeEquation.setBoundaryYp(new Expression(diffusionRateName + "*" + exactDyName));
                    }
                } else {
                    throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionYp());
                }
                if (compartmentSubDomain.getBoundaryConditionZm().isDIRICHLET()) {
                    Expression origExp = pdeEquation.getBoundaryZm();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryZm(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
                    } else {
                        pdeEquation.setBoundaryZm(exactExp);
                    }
                } else if (compartmentSubDomain.getBoundaryConditionZm().isNEUMANN()) {
                    Expression origExp = pdeEquation.getBoundaryZm();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryZm(new Expression(origExp + "-" + substitutedExp + "-" + diffusionRateName + "*" + exactDzName));
                    } else {
                        pdeEquation.setBoundaryZm(new Expression("-" + diffusionRateName + "*" + exactDzName));
                    }
                } else {
                    throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionXm());
                }
                if (compartmentSubDomain.getBoundaryConditionZp().isDIRICHLET()) {
                    Expression origExp = pdeEquation.getBoundaryZp();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryZp(new Expression(origExp + "-" + substitutedExp + "+" + exactExp));
                    } else {
                        pdeEquation.setBoundaryZp(exactExp);
                    }
                } else if (compartmentSubDomain.getBoundaryConditionZp().isNEUMANN()) {
                    Expression origExp = pdeEquation.getBoundaryZp();
                    if (origExp != null) {
                        Expression substitutedExp = substituteWithExactSolution(origExp, compartmentSubDomain, exactMath);
                        pdeEquation.setBoundaryZp(new Expression(origExp + "-" + substitutedExp + "+" + diffusionRateName + "*" + exactDzName));
                    } else {
                        pdeEquation.setBoundaryZp(new Expression(diffusionRateName + "*" + exactDzName));
                    }
                } else {
                    throw new RuntimeException("unsupported boundary condition type " + compartmentSubDomain.getBoundaryConditionZp());
                }
            } else {
                throw new RuntimeException("SolverTest.constructedExactMath(): equation type " + equation.getClass().getName() + " not yet implemented");
            }
        }
        if (subDomain instanceof MembraneSubDomain) {
            MembraneSubDomain membraneSubDomain = (MembraneSubDomain) subDomain;
            Enumeration<JumpCondition> enumJumpConditions = membraneSubDomain.getJumpConditions();
            while (enumJumpConditions.hasMoreElements()) {
                JumpCondition jumpCondition = enumJumpConditions.nextElement();
                Expression origInfluxExp = jumpCondition.getInFluxExpression();
                Expression origOutfluxExp = jumpCondition.getOutFluxExpression();
                Expression substitutedInfluxExp = substituteWithExactSolution(origInfluxExp, membraneSubDomain, exactMath);
                Expression substitutedOutfluxExp = substituteWithExactSolution(origOutfluxExp, membraneSubDomain, exactMath);
                String varName = jumpCondition.getVariable().getName();
                String origInfluxName = "_" + varName + "_" + subDomain.getName() + "_origInflux";
                String origOutfluxName = "_" + varName + "_" + subDomain.getName() + "_origOutflux";
                String substitutedInfluxName = "_" + varName + "_" + subDomain.getName() + "_substitutedInflux";
                String substitutedOutfluxName = "_" + varName + "_" + subDomain.getName() + "_substitutedOutflux";
                String diffusionRateInsideName = "_" + varName + "_" + membraneSubDomain.getInsideCompartment().getName() + "_diffusionRate";
                String diffusionRateOutsideName = "_" + varName + "_" + membraneSubDomain.getOutsideCompartment().getName() + "_diffusionRate";
                String exactInsideDxName = "_" + varName + "_" + membraneSubDomain.getInsideCompartment().getName() + "_exact_dx";
                String exactInsideDyName = "_" + varName + "_" + membraneSubDomain.getInsideCompartment().getName() + "_exact_dy";
                String exactInsideDzName = "_" + varName + "_" + membraneSubDomain.getInsideCompartment().getName() + "_exact_dz";
                String exactOutsideDxName = "_" + varName + "_" + membraneSubDomain.getOutsideCompartment().getName() + "_exact_dx";
                String exactOutsideDyName = "_" + varName + "_" + membraneSubDomain.getOutsideCompartment().getName() + "_exact_dy";
                String exactOutsideDzName = "_" + varName + "_" + membraneSubDomain.getOutsideCompartment().getName() + "_exact_dz";
                String outwardNormalXName = "_" + membraneSubDomain.getInsideCompartment().getName() + "_Nx";
                String outwardNormalYName = "_" + membraneSubDomain.getInsideCompartment().getName() + "_Ny";
                String outwardNormalZName = "_" + membraneSubDomain.getInsideCompartment().getName() + "_Nz";
                String exactInfluxName = "_" + varName + "_" + membraneSubDomain.getName() + "_exactInflux";
                String exactOutfluxName = "_" + varName + "_" + membraneSubDomain.getName() + "_exactOutflux";
                Expression exactInfluxExp = new Expression(diffusionRateInsideName + " * (" + outwardNormalXName + "*" + exactInsideDxName + " + " + outwardNormalYName + "*" + exactInsideDyName + " + " + outwardNormalZName + "*" + exactInsideDzName + ")");
                Expression exactOutfluxExp = new Expression("-" + diffusionRateOutsideName + " * (" + outwardNormalXName + "*" + exactOutsideDxName + " + " + outwardNormalYName + "*" + exactOutsideDyName + " + " + outwardNormalZName + "*" + exactOutsideDzName + ")");
                Expression newInfluxExp = new Expression(origInfluxName + " - " + substitutedInfluxName + " + " + exactInfluxName);
                Expression newOutfluxExp = new Expression(origOutfluxName + " - " + substitutedOutfluxName + " + " + exactOutfluxName);
                varHash.addVariable(new Function(origInfluxName, origInfluxExp, domain));
                varHash.addVariable(new Function(origOutfluxName, origOutfluxExp, domain));
                varHash.addVariable(new Function(exactInfluxName, exactInfluxExp, domain));
                varHash.addVariable(new Function(exactOutfluxName, exactOutfluxExp, domain));
                varHash.addVariable(new Function(substitutedInfluxName, substitutedInfluxExp, domain));
                varHash.addVariable(new Function(substitutedOutfluxName, substitutedOutfluxExp, domain));
                jumpCondition.setInFlux(newInfluxExp);
                jumpCondition.setOutFlux(newOutfluxExp);
            }
        }
    }
    exactMath.setAllVariables(varHash.getAlphabeticallyOrderedVariables());
    if (!exactMath.isValid()) {
        throw new RuntimeException("generated Math is not valid: " + exactMath.getWarning());
    }
    return exactMath;
}

Example 44 with Variable

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

the class FVSolverStandalone method writeVCGAndResampleFieldData.

protected void writeVCGAndResampleFieldData() throws SolverException {
    fireSolverStarting(SimulationMessage.MESSAGE_SOLVEREVENT_STARTING_PROC_GEOM);
    try {
        // write subdomains file
        SubdomainInfo.write(new File(getSaveDirectory(), baseName + SimDataConstants.SUBDOMAINS_FILE_SUFFIX), simTask.getSimulation().getMathDescription());
        PrintWriter pw = new PrintWriter(new FileWriter(new File(getSaveDirectory(), baseName + SimDataConstants.VCG_FILE_EXTENSION)));
        GeometryFileWriter.write(pw, getResampledGeometry());
        pw.close();
        FieldDataIdentifierSpec[] argFieldDataIDSpecs = simTask.getSimulationJob().getFieldDataIdentifierSpecs();
        if (argFieldDataIDSpecs != null && argFieldDataIDSpecs.length > 0) {
            fireSolverStarting(SimulationMessage.MESSAGE_SOLVEREVENT_STARTING_RESAMPLE_FD);
            FieldFunctionArguments psfFieldFunc = null;
            Variable var = simTask.getSimulationJob().getSimulationSymbolTable().getVariable(Simulation.PSF_FUNCTION_NAME);
            if (var != null) {
                FieldFunctionArguments[] ffas = FieldUtilities.getFieldFunctionArguments(var.getExpression());
                if (ffas == null || ffas.length == 0) {
                    throw new DataAccessException("Point Spread Function " + Simulation.PSF_FUNCTION_NAME + " can only be a single field function.");
                } else {
                    Expression newexp;
                    try {
                        newexp = new Expression(ffas[0].infix());
                        if (!var.getExpression().compareEqual(newexp)) {
                            throw new DataAccessException("Point Spread Function " + Simulation.PSF_FUNCTION_NAME + " can only be a single field function.");
                        }
                        psfFieldFunc = ffas[0];
                    } catch (ExpressionException e) {
                        e.printStackTrace();
                        throw new DataAccessException(e.getMessage());
                    }
                }
            }
            boolean[] bResample = new boolean[argFieldDataIDSpecs.length];
            Arrays.fill(bResample, true);
            for (int i = 0; i < argFieldDataIDSpecs.length; i++) {
                argFieldDataIDSpecs[i].getFieldFuncArgs().getTime().bindExpression(simTask.getSimulationJob().getSimulationSymbolTable());
                if (argFieldDataIDSpecs[i].getFieldFuncArgs().equals(psfFieldFunc)) {
                    bResample[i] = false;
                }
            }
            int numMembraneElements = getResampledGeometry().getGeometrySurfaceDescription().getSurfaceCollection().getTotalPolygonCount();
            CartesianMesh simpleMesh = CartesianMesh.createSimpleCartesianMesh(getResampledGeometry().getOrigin(), getResampledGeometry().getExtent(), simTask.getSimulation().getMeshSpecification().getSamplingSize(), getResampledGeometry().getGeometrySurfaceDescription().getRegionImage());
            String secondarySimDataDir = PropertyLoader.getProperty(PropertyLoader.secondarySimDataDirInternalProperty, null);
            DataSetControllerImpl dsci = new DataSetControllerImpl(null, getSaveDirectory().getParentFile(), secondarySimDataDir == null ? null : new File(secondarySimDataDir));
            dsci.writeFieldFunctionData(null, argFieldDataIDSpecs, bResample, simpleMesh, simResampleInfoProvider, numMembraneElements, HESM_OVERWRITE_AND_CONTINUE);
        }
    } catch (Exception e) {
        throw new SolverException(e.getMessage());
    }
}

Example 45 with Variable

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

the class NetCDFWriter method getVariableSymbols.

private String[] getVariableSymbols(String[] symbols) {
    Simulation simulation = simTask.getSimulation();
    Vector<String> vars = new Vector<String>();
    if (symbols != null) {
        for (int i = 0; i < symbols.length; i++) {
            Variable v = simulation.getMathDescription().getVariable(symbols[i]);
            if ((v != null) && (v instanceof StochVolVariable)) {
                vars.add(symbols[i]);
            }
        }
        return vars.toArray(new String[vars.size()]);
    } else
        return vars.toArray(new String[0]);
}