Search in sources :

Example 21 with Function

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

the class SundialsSolver method getODESolverResultSet.

/**
 * This method was created in VisualAge.
 * @return double[]
 * @param vectorIndex int
 */
public ODESolverResultSet getODESolverResultSet() {
    // 
    // read .ida file
    // 
    ODESolverResultSet odeSolverResultSet = getStateVariableResultSet();
    if (odeSolverResultSet == null) {
        return null;
    }
    // 
    // add appropriate Function columns to result set
    // 
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    Function[] functions = simSymbolTable.getFunctions();
    for (int i = 0; i < functions.length; i++) {
        if (SimulationSymbolTable.isFunctionSaved(functions[i])) {
            Expression exp1 = new Expression(functions[i].getExpression());
            try {
                exp1 = simSymbolTable.substituteFunctions(exp1);
            } catch (MathException e) {
                e.printStackTrace(System.out);
                throw new RuntimeException("Substitute function failed on function " + functions[i].getName() + " " + e.getMessage());
            } catch (ExpressionException e) {
                e.printStackTrace(System.out);
                throw new RuntimeException("Substitute function failed on function " + functions[i].getName() + " " + e.getMessage());
            }
            try {
                FunctionColumnDescription cd = new FunctionColumnDescription(exp1.flatten(), functions[i].getName(), null, functions[i].getName(), false);
                odeSolverResultSet.addFunctionColumn(cd);
            } catch (ExpressionException e) {
                e.printStackTrace(System.out);
            }
        }
    }
    if (getSensitivityParameter() != null) {
        try {
            if (odeSolverResultSet.findColumn(getSensitivityParameter().getName()) == -1) {
                FunctionColumnDescription fcd = new FunctionColumnDescription(new Expression(getSensitivityParameter().getConstantValue()), getSensitivityParameter().getName(), null, getSensitivityParameter().getName(), false);
                odeSolverResultSet.addFunctionColumn(fcd);
            }
            Variable[] variables = simSymbolTable.getVariables();
            StateVariable[] stateVars = createStateVariables();
            for (int i = 0; i < variables.length; i++) {
                if (variables[i] instanceof Function && SimulationSymbolTable.isFunctionSaved((Function) variables[i])) {
                    Function depSensFunction = (Function) variables[i];
                    Expression depSensFnExpr = new Expression(depSensFunction.getExpression());
                    depSensFnExpr = simSymbolTable.substituteFunctions(depSensFnExpr);
                    depSensFnExpr = getFunctionSensitivity(depSensFnExpr, getSensitivityParameter(), stateVars);
                    // depSensFnExpr = depSensFnExpr.flatten(); 	// already bound and flattened in getFunctionSensitivity, no need here.....
                    String depSensFnName = new String("sens_" + depSensFunction.getName() + "_wrt_" + getSensitivityParameter().getName());
                    if (depSensFunction != null) {
                        FunctionColumnDescription cd = new FunctionColumnDescription(depSensFnExpr.flatten(), depSensFnName, getSensitivityParameter().getName(), depSensFnName, false);
                        odeSolverResultSet.addFunctionColumn(cd);
                    }
                }
            }
        } catch (MathException e) {
            e.printStackTrace(System.out);
            throw new RuntimeException("Error adding function to resultSet: " + e.getMessage());
        } catch (ExpressionException e) {
            e.printStackTrace(System.out);
            throw new RuntimeException("Error adding function to resultSet: " + e.getMessage());
        }
    }
    return odeSolverResultSet;
}
Also used : Variable(cbit.vcell.math.Variable) VolVariable(cbit.vcell.math.VolVariable) SimulationSymbolTable(cbit.vcell.solver.SimulationSymbolTable) ExpressionException(cbit.vcell.parser.ExpressionException) Function(cbit.vcell.math.Function) AnnotatedFunction(cbit.vcell.solver.AnnotatedFunction) Expression(cbit.vcell.parser.Expression) MathException(cbit.vcell.math.MathException) FunctionColumnDescription(cbit.vcell.math.FunctionColumnDescription)

Example 22 with Function

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

the class SimulationData method getVarAndFunctionDataIdentifiers.

/**
 * This method was created in VisualAge.
 * @return java.lang.String[]
 */
public synchronized DataIdentifier[] getVarAndFunctionDataIdentifiers(OutputContext outputContext) throws IOException, DataAccessException {
    // Is this zip format?
    boolean bIsChombo = false;
    try {
        bIsChombo = isChombo();
    } catch (FileNotFoundException e) {
        e.printStackTrace(System.out);
    }
    File zipFile1 = getZipFile(bIsChombo, null);
    File zipFile2 = getZipFile(bIsChombo, 0);
    bZipFormat1 = false;
    bZipFormat2 = false;
    if (zipFile1.exists()) {
        bZipFormat1 = true;
    } else if (zipFile2.exists()) {
        bZipFormat2 = true;
    }
    refreshLogFile();
    if (!isComsol()) {
        try {
            refreshMeshFile();
        } catch (MathException e) {
            e.printStackTrace(System.out);
            throw new DataAccessException(e.getMessage());
        }
    }
    if (!isRulesData && !getIsODEData() && !isComsol() && dataFilenames != null) {
        // read variables only when I have never read the file since variables don't change
        if (dataSetIdentifierList.size() == 0) {
            File file = getPDEDataFile(0.0);
            DataSet dataSet = getPDEDataSet(file, 0.0);
            String[] varNames = dataSet.getDataNames();
            int[] varTypeInts = dataSet.getVariableTypeIntegers();
            if (varNames == null) {
                return null;
            }
            dataSetIdentifierList.clear();
            for (int i = 0; i < varNames.length; i++) {
                VariableType varType = null;
                try {
                    varType = VariableType.getVariableTypeFromInteger(varTypeInts[i]);
                } catch (IllegalArgumentException e) {
                    if (LG.isEnabledFor(Level.WARN)) {
                        LG.warn("Exception typing " + varNames[i] + " has unsupported type " + varTypeInts[i] + ": " + e.getMessage());
                    }
                    varType = SimulationData.getVariableTypeFromLength(mesh, dataSet.getDataLength(varNames[i]));
                }
                Domain domain = Variable.getDomainFromCombinedIdentifier(varNames[i]);
                String varName = Variable.getNameFromCombinedIdentifier(varNames[i]);
                dataSetIdentifierList.addElement(new DataSetIdentifier(varName, varType, domain));
            }
            refreshDataProcessingOutputInfo(outputContext);
            if (dataProcessingOutputInfo != null) {
                for (int i = 0; i < dataProcessingOutputInfo.getVariableNames().length; i++) {
                    if (dataProcessingOutputInfo.getPostProcessDataType(dataProcessingOutputInfo.getVariableNames()[i]).equals(DataProcessingOutputInfo.PostProcessDataType.image)) {
                        dataSetIdentifierList.addElement(new DataSetIdentifier(dataProcessingOutputInfo.getVariableNames()[i], VariableType.POSTPROCESSING, null));
                    }
                }
            }
        }
        // always read functions file since functions might change
        getFunctionDataIdentifiers(outputContext);
    }
    if ((isRulesData || getIsODEData()) && dataSetIdentifierList.size() == 0) {
        ODEDataBlock odeDataBlock = getODEDataBlock();
        if (odeDataBlock == null) {
            throw new DataAccessException("Results are not availabe yet. Please try again later.");
        }
        ODESimData odeSimData = odeDataBlock.getODESimData();
        int colCount = odeSimData.getColumnDescriptionsCount();
        // assume index=0 is time "t"
        int DATA_OFFSET = 1;
        dataSetIdentifierList.clear();
        for (int i = 0; i < (colCount - DATA_OFFSET); i++) {
            String varName = odeSimData.getColumnDescriptions(i + DATA_OFFSET).getDisplayName();
            // TODO domain
            Domain domain = null;
            dataSetIdentifierList.addElement(new DataSetIdentifier(varName, VariableType.NONSPATIAL, domain));
        }
    }
    if (isComsol() && dataSetIdentifierList.size() == 0) {
        ComsolSimFiles comsolSimFiles = getComsolSimFiles();
        if (comsolSimFiles.simTaskXMLFile != null) {
            try {
                String xmlString = FileUtils.readFileToString(comsolSimFiles.simTaskXMLFile);
                SimulationTask simTask = XmlHelper.XMLToSimTask(xmlString);
                Enumeration<Variable> variablesEnum = simTask.getSimulation().getMathDescription().getVariables();
                while (variablesEnum.hasMoreElements()) {
                    Variable var = variablesEnum.nextElement();
                    if (var instanceof VolVariable) {
                        dataSetIdentifierList.addElement(new DataSetIdentifier(var.getName(), VariableType.VOLUME, var.getDomain()));
                    } else if (var instanceof MemVariable) {
                        dataSetIdentifierList.addElement(new DataSetIdentifier(var.getName(), VariableType.MEMBRANE, var.getDomain()));
                    } else if (var instanceof Function) {
                        VariableType varType = VariableType.UNKNOWN;
                        if (var.getDomain() != null && var.getDomain().getName() != null) {
                            SubDomain subDomain = simTask.getSimulation().getMathDescription().getSubDomain(var.getDomain().getName());
                            if (subDomain instanceof CompartmentSubDomain) {
                                varType = VariableType.VOLUME;
                            } else if (subDomain instanceof MembraneSubDomain) {
                                varType = VariableType.MEMBRANE;
                            } else if (subDomain instanceof FilamentSubDomain) {
                                throw new RuntimeException("filament subdomains not supported");
                            } else if (subDomain instanceof PointSubDomain) {
                                varType = VariableType.POINT_VARIABLE;
                            }
                        }
                        dataSetIdentifierList.addElement(new DataSetIdentifier(var.getName(), varType, var.getDomain()));
                    } else if (var instanceof Constant) {
                        System.out.println("ignoring Constant " + var.getName());
                    } else if (var instanceof InsideVariable) {
                        System.out.println("ignoring InsideVariable " + var.getName());
                    } else if (var instanceof OutsideVariable) {
                        System.out.println("ignoring OutsideVariable " + var.getName());
                    } else {
                        throw new RuntimeException("unexpected variable " + var.getName() + " of type " + var.getClass().getName());
                    }
                }
            } catch (XmlParseException | ExpressionException e) {
                e.printStackTrace();
                throw new RuntimeException("failed to read sim task file, msg: " + e.getMessage(), e);
            }
        }
    }
    DataIdentifier[] dis = new DataIdentifier[dataSetIdentifierList.size()];
    for (int i = 0; i < dataSetIdentifierList.size(); i++) {
        DataSetIdentifier dsi = (DataSetIdentifier) dataSetIdentifierList.elementAt(i);
        String displayName = dsi.getName();
        if (dsi.isFunction()) {
            AnnotatedFunction f = null;
            for (int j = 0; j < annotatedFunctionList.size(); j++) {
                AnnotatedFunction function = (AnnotatedFunction) annotatedFunctionList.elementAt(j);
                if (function.getName().equals(dsi.getName())) {
                    f = function;
                    break;
                }
            }
            if (f != null) {
                displayName = f.getDisplayName();
            }
        }
        dis[i] = new DataIdentifier(dsi.getName(), dsi.getVariableType(), dsi.getDomain(), dsi.isFunction(), displayName);
    }
    return dis;
}
Also used : MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) SimulationTask(cbit.vcell.messaging.server.SimulationTask) InsideVariable(cbit.vcell.math.InsideVariable) VolVariable(cbit.vcell.math.VolVariable) ReservedVariable(cbit.vcell.math.ReservedVariable) MemVariable(cbit.vcell.math.MemVariable) OutsideVariable(cbit.vcell.math.OutsideVariable) Variable(cbit.vcell.math.Variable) VCDataIdentifier(org.vcell.util.document.VCDataIdentifier) VCSimulationDataIdentifier(cbit.vcell.solver.VCSimulationDataIdentifier) ExternalDataIdentifier(org.vcell.util.document.ExternalDataIdentifier) Constant(cbit.vcell.math.Constant) FileNotFoundException(java.io.FileNotFoundException) PointSubDomain(cbit.vcell.math.PointSubDomain) ODESimData(cbit.vcell.solver.ode.ODESimData) InsideVariable(cbit.vcell.math.InsideVariable) ExpressionException(cbit.vcell.parser.ExpressionException) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) FilamentSubDomain(cbit.vcell.math.FilamentSubDomain) SubDomain(cbit.vcell.math.SubDomain) PointSubDomain(cbit.vcell.math.PointSubDomain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) Function(cbit.vcell.math.Function) AnnotatedFunction(cbit.vcell.solver.AnnotatedFunction) MemVariable(cbit.vcell.math.MemVariable) DataAccessException(org.vcell.util.DataAccessException) AnnotatedFunction(cbit.vcell.solver.AnnotatedFunction) ComsolSimFiles(org.vcell.vis.io.ComsolSimFiles) VariableType(cbit.vcell.math.VariableType) VolVariable(cbit.vcell.math.VolVariable) XmlParseException(cbit.vcell.xml.XmlParseException) FilamentSubDomain(cbit.vcell.math.FilamentSubDomain) MathException(cbit.vcell.math.MathException) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) OutsideVariable(cbit.vcell.math.OutsideVariable) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) FilamentSubDomain(cbit.vcell.math.FilamentSubDomain) SubDomain(cbit.vcell.math.SubDomain) PointSubDomain(cbit.vcell.math.PointSubDomain) Domain(cbit.vcell.math.Variable.Domain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) ZipFile(org.apache.commons.compress.archivers.zip.ZipFile) File(java.io.File)

Example 23 with Function

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

the class OutputFunctionContext method computeFunctionTypeWRTExpression.

// check if the new expression is valid for outputFunction of functionType
public VariableType computeFunctionTypeWRTExpression(AnnotatedFunction outputFunction, Expression exp) throws ExpressionException, InconsistentDomainException {
    MathDescription mathDescription = getSimulationOwner().getMathDescription();
    boolean bSpatial = getSimulationOwner().getGeometry().getDimension() > 0;
    if (!bSpatial) {
        return VariableType.NONSPATIAL;
    }
    Expression newexp = new Expression(exp);
    // making sure that output function is not direct function of constant.
    newexp.bindExpression(this);
    // here use math description as symbol table because we allow
    // new expression itself to be function of constant.
    newexp = MathUtilities.substituteFunctions(newexp, this).flatten();
    String[] symbols = newexp.getSymbols();
    VariableType functionType = outputFunction.getFunctionType();
    String funcName = outputFunction.getName();
    Domain funcDomain = outputFunction.getDomain();
    VariableType[] varTypes = null;
    if (symbols != null && symbols.length > 0) {
        // making sure that new expression is defined in the same domain
        varTypes = new VariableType[symbols.length];
        for (int i = 0; i < symbols.length; i++) {
            if (ReservedMathSymbolEntries.getReservedVariableEntry(symbols[i]) != null) {
                varTypes[i] = functionType;
            } else {
                Variable var = mathDescription.getVariable(symbols[i]);
                if (var == null) {
                    var = mathDescription.getPostProcessingBlock().getDataGenerator(symbols[i]);
                }
                varTypes[i] = VariableType.getVariableType(var);
                if (funcDomain != null) {
                    if (var.getDomain() == null) {
                        // OK
                        continue;
                    }
                    GeometryClass funcGeoClass = simulationOwner.getGeometry().getGeometryClass(funcDomain.getName());
                    GeometryClass varGeoClass = simulationOwner.getGeometry().getGeometryClass(var.getDomain().getName());
                    if (varGeoClass instanceof SubVolume && funcGeoClass instanceof SurfaceClass) {
                        // seems ok if membrane refereces volume
                        if (!((SurfaceClass) funcGeoClass).isAdjacentTo((SubVolume) varGeoClass)) {
                            // but has to be adjacent
                            String errMsg = "'" + funcName + "' defined on Membrane '" + funcDomain.getName() + "' directly or indirectly references " + " variable '" + symbols[i] + "' defined on Volume '" + var.getDomain().getName() + " which is not adjacent to Membrane '" + funcDomain.getName() + "'.";
                            throw new ExpressionException(errMsg);
                        }
                    } else if (!var.getDomain().compareEqual(funcDomain)) {
                        String errMsg = "'" + funcName + "' defined on '" + funcDomain.getName() + "' directly or indirectly references " + " variable '" + symbols[i] + "' defined on '" + var.getDomain().getName() + ".";
                        throw new ExpressionException(errMsg);
                    }
                }
            }
        }
    }
    // if there are no variables (like built in function, vcRegionArea), check with flattened expression to find out the variable type of the new expression
    VariableDomain functionVariableDomain = functionType.getVariableDomain();
    Function flattenedFunction = new Function(funcName, newexp, funcDomain);
    flattenedFunction.bind(this);
    VariableType newVarType = SimulationSymbolTable.getFunctionVariableType(flattenedFunction, getSimulationOwner().getMathDescription(), symbols, varTypes, bSpatial);
    if (!newVarType.getVariableDomain().equals(functionVariableDomain)) {
        String errMsg = "The expression for '" + funcName + "' includes at least one " + newVarType.getVariableDomain().getName() + " variable. Please make sure that only " + functionVariableDomain.getName() + " variables are " + "referenced in " + functionVariableDomain.getName() + " output functions.";
        throw new ExpressionException(errMsg);
    }
    return newVarType;
}
Also used : GeometryClass(cbit.vcell.geometry.GeometryClass) ReservedVariable(cbit.vcell.math.ReservedVariable) InsideVariable(cbit.vcell.math.InsideVariable) OutsideVariable(cbit.vcell.math.OutsideVariable) Variable(cbit.vcell.math.Variable) VariableDomain(cbit.vcell.math.VariableType.VariableDomain) VariableType(cbit.vcell.math.VariableType) MathDescription(cbit.vcell.math.MathDescription) SurfaceClass(cbit.vcell.geometry.SurfaceClass) ExpressionException(cbit.vcell.parser.ExpressionException) Function(cbit.vcell.math.Function) Expression(cbit.vcell.parser.Expression) SubVolume(cbit.vcell.geometry.SubVolume) VariableDomain(cbit.vcell.math.VariableType.VariableDomain) Domain(cbit.vcell.math.Variable.Domain)

Example 24 with Function

use of cbit.vcell.math.Function 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);
}
Also used : MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable) InsideVariable(cbit.vcell.math.InsideVariable) SensVariable(cbit.vcell.solver.ode.SensVariable) FilamentVariable(cbit.vcell.math.FilamentVariable) VolVariable(cbit.vcell.math.VolVariable) VolumeRegionVariable(cbit.vcell.math.VolumeRegionVariable) ReservedVariable(cbit.vcell.math.ReservedVariable) MemVariable(cbit.vcell.math.MemVariable) OutsideVariable(cbit.vcell.math.OutsideVariable) FilamentRegionVariable(cbit.vcell.math.FilamentRegionVariable) Variable(cbit.vcell.math.Variable) VariableHash(cbit.vcell.math.VariableHash) Function(cbit.vcell.math.Function) Expression(cbit.vcell.parser.Expression) Vector(java.util.Vector)

Example 25 with Function

use of cbit.vcell.math.Function 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;
}
Also used : MembraneRegionVariable(cbit.vcell.math.MembraneRegionVariable) InsideVariable(cbit.vcell.math.InsideVariable) SensVariable(cbit.vcell.solver.ode.SensVariable) FilamentVariable(cbit.vcell.math.FilamentVariable) VolVariable(cbit.vcell.math.VolVariable) VolumeRegionVariable(cbit.vcell.math.VolumeRegionVariable) ReservedVariable(cbit.vcell.math.ReservedVariable) MemVariable(cbit.vcell.math.MemVariable) OutsideVariable(cbit.vcell.math.OutsideVariable) FilamentRegionVariable(cbit.vcell.math.FilamentRegionVariable) Variable(cbit.vcell.math.Variable) VariableHash(cbit.vcell.math.VariableHash) Constant(cbit.vcell.math.Constant) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) Function(cbit.vcell.math.Function) Vector(java.util.Vector) VolVariable(cbit.vcell.math.VolVariable) OdeEquation(cbit.vcell.math.OdeEquation) PdeEquation(cbit.vcell.math.PdeEquation) Equation(cbit.vcell.math.Equation) FastRate(cbit.vcell.math.FastRate) FastInvariant(cbit.vcell.math.FastInvariant) OdeEquation(cbit.vcell.math.OdeEquation) Expression(cbit.vcell.parser.Expression) FastSystem(cbit.vcell.math.FastSystem) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) Domain(cbit.vcell.math.Variable.Domain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain)

Aggregations

Function (cbit.vcell.math.Function)52 Expression (cbit.vcell.parser.Expression)42 Variable (cbit.vcell.math.Variable)32 ExpressionException (cbit.vcell.parser.ExpressionException)26 Constant (cbit.vcell.math.Constant)24 VolVariable (cbit.vcell.math.VolVariable)22 CompartmentSubDomain (cbit.vcell.math.CompartmentSubDomain)18 MathException (cbit.vcell.math.MathException)16 MathDescription (cbit.vcell.math.MathDescription)15 InsideVariable (cbit.vcell.math.InsideVariable)13 MemVariable (cbit.vcell.math.MemVariable)13 OutsideVariable (cbit.vcell.math.OutsideVariable)13 ReservedVariable (cbit.vcell.math.ReservedVariable)13 Domain (cbit.vcell.math.Variable.Domain)13 MembraneRegionVariable (cbit.vcell.math.MembraneRegionVariable)12 SubDomain (cbit.vcell.math.SubDomain)12 VolumeRegionVariable (cbit.vcell.math.VolumeRegionVariable)12 Vector (java.util.Vector)12 Element (org.jdom.Element)11 SubVolume (cbit.vcell.geometry.SubVolume)9