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

Example 21 with VolVariable

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

the class AbstractSolver method getFunctionSensitivity.

/**
 * This method was created by a SmartGuide.
 * @return double[]
 * @param identifier java.lang.String
 */
public Expression getFunctionSensitivity(Expression funcExpr, Constant constant, StateVariable[] stateVariables) throws ExpressionException {
    if (stateVariables == null || stateVariables.length == 0) {
        return null;
    }
    // this uses the chain rule
    // d F(x)   del F(x)        del F(x)   del xi
    // ------ = -------- + Sum (-------- * ------)
    // d k      del k      i    del xi    del k
    // explicit     via state variables
    // dependence
    // 
    // collect the explicit term
    // 
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    Expression sensFuncExp = funcExpr.differentiate(constant.getName());
    sensFuncExp.bindExpression(simSymbolTable);
    sensFuncExp = sensFuncExp.flatten();
    for (int i = 0; i < stateVariables.length; i++) {
        if (stateVariables[i] instanceof ODEStateVariable) {
            ODEStateVariable sv = (ODEStateVariable) stateVariables[i];
            VolVariable volVar = (VolVariable) sv.getVariable();
            // 
            // get corresponding sensitivity state variable
            // 
            SensStateVariable sensStateVariable = null;
            for (int j = 0; j < stateVariables.length; j++) {
                if (stateVariables[j] instanceof SensStateVariable && ((SensVariable) ((SensStateVariable) stateVariables[j]).getVariable()).getVolVariable().getName().equals(volVar.getName()) && ((SensStateVariable) stateVariables[j]).getParameter().getName().equals(constant.getName())) {
                    sensStateVariable = (SensStateVariable) stateVariables[j];
                }
            }
            if (sensStateVariable == null) {
                System.out.println("sens of " + volVar.getName() + " wrt " + constant.getName() + " is not found");
                return null;
            }
            // 
            // get coefficient of proportionality   (e.g.  A = total + b*B + c*C ... gives dA/dK = b*dB/dK + c*dC/dK)
            // 
            Expression coeffExpression = funcExpr.differentiate(volVar.getName());
            coeffExpression.bindExpression(simSymbolTable);
            coeffExpression = MathUtilities.substituteFunctions(coeffExpression, simSymbolTable);
            coeffExpression.bindExpression(simSymbolTable);
            coeffExpression = coeffExpression.flatten();
            sensFuncExp = Expression.add(sensFuncExp, Expression.mult(coeffExpression, new Expression(sensStateVariable.getVariable().getName())));
        }
    }
    return new Expression(sensFuncExp);
}

Example 22 with VolVariable

use of cbit.vcell.math.VolVariable 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 23 with VolVariable

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

the class DisplayImageOp method displayImage.

public void displayImage(final Image image, String title, WindowListener listener) {
    final ImagePlaneManagerPanel imagePanel = new ImagePlaneManagerPanel();
    double[] doublePixels = image.getDoublePixels();
    double minPixel = Double.MAX_VALUE;
    double maxPixel = -Double.MAX_VALUE;
    for (int i = 0; i < doublePixels.length; i++) {
        double pixel = doublePixels[i];
        doublePixels[i] = pixel;
        minPixel = Math.min(minPixel, pixel);
        maxPixel = Math.max(maxPixel, pixel);
    }
    Range newRange = new Range(minPixel, maxPixel);
    SourceDataInfo source = new SourceDataInfo(SourceDataInfo.RAW_VALUE_TYPE, doublePixels, image.getExtent(), image.getOrigin(), newRange, 0, image.getNumX(), 1, image.getNumY(), image.getNumX(), image.getNumZ(), image.getNumX() * image.getNumY());
    imagePanel.setDisplayAdapterServicePanelVisible(true);
    imagePanel.setCurveValueProvider(new CurveValueProvider() {

        @Override
        public void curveAdded(Curve curve) {
            System.out.println("called curveAdded(" + curve + "), do nothing for now");
        }

        @Override
        public void curveRemoved(Curve curve) {
            System.out.println("called curveRemoved(" + curve + ")");
        }

        @Override
        public String getCurveValue(CurveSelectionInfo csi) {
            System.out.println("called getCurveValue(CurveSelectionInfo " + csi);
            return null;
        }

        @Override
        public CurveSelectionInfo getInitalCurveSelection(int tool, Coordinate wc) {
            System.out.println("called getInitialCurveSelection(tool=" + tool + ", coord=" + wc + ")");
            return null;
        }

        @Override
        public boolean isAddControlPointOK(int tool, Coordinate wc, Curve addedToThisCurve) {
            System.out.println("called isAddControlPointOK");
            return true;
        }

        @Override
        public boolean providesInitalCurve(int tool, Coordinate wc) {
            System.out.println("called providesInitialCurve(tool=" + tool + " (TOOL_LINE=" + CurveEditorTool.TOOL_LINE + "), coord=" + wc);
            return false;
        }

        @Override
        public void setDescription(Curve curve) {
            System.out.println("called setDescription(" + curve + ")");
            curve.setDescription(CurveValueProvider.DESCRIPTION_VOLUME);
        }

        @Override
        public CurveSelectionInfo findChomboCurveSelectionInfoForPoint(CoordinateIndex ci) {
            System.out.println("called find ChomboCurveSelectionInfoForPoint(coord=" + ci + ")");
            return null;
        }
    });
    DisplayAdapterService das = imagePanel.getDisplayAdapterServicePanel().getDisplayAdapterService();
    das.setValueDomain(null);
    das.addColorModelForValues(DisplayAdapterService.createGrayColorModel(), DisplayAdapterService.createGraySpecialColors(), DisplayAdapterService.GRAY);
    das.addColorModelForValues(DisplayAdapterService.createBlueRedColorModel(), DisplayAdapterService.createBlueRedSpecialColors(), DisplayAdapterService.BLUERED);
    das.setActiveColorModelID(DisplayAdapterService.BLUERED);
    final JFrame jframe = new JFrame();
    jframe.setTitle(title);
    JPanel panel = new JPanel();
    panel.setLayout(new GridBagLayout());
    GridBagConstraints imageConstraints = new GridBagConstraints();
    imageConstraints.gridx = 0;
    imageConstraints.gridy = 0;
    imageConstraints.weightx = 1.0;
    imageConstraints.weighty = 1.0;
    imageConstraints.fill = GridBagConstraints.BOTH;
    panel.add(imagePanel, imageConstraints);
    JButton plotButton = new JButton("plot");
    plotButton.addActionListener(new ActionListener() {

        @Override
        public void actionPerformed(ActionEvent e) {
            try {
                Curve curve = imagePanel.getCurveRenderer().getSelection().getCurve();
                VariableType variableType = VariableType.VOLUME;
                Curve samplerCurve = curve.getSampledCurve();
                samplerCurve.setDescription(curve.getDescription());
                VCImage vcImage = new VCImageUncompressed(null, new byte[image.getISize().getXYZ()], image.getExtent(), image.getISize().getX(), image.getISize().getY(), image.getISize().getZ());
                int dimension = 1 + (image.getISize().getY() > 0 ? 1 : 0) + (image.getISize().getZ() > 0 ? 1 : 0);
                RegionImage regionImage = new RegionImage(vcImage, dimension, image.getExtent(), image.getOrigin(), RegionImage.NO_SMOOTHING);
                CartesianMesh mesh = CartesianMesh.createSimpleCartesianMesh(image.getOrigin(), image.getExtent(), image.getISize(), regionImage);
                SpatialSelectionVolume ssVolume = new SpatialSelectionVolume(new CurveSelectionInfo(samplerCurve), variableType, mesh);
                String varName = "var";
                SymbolTableEntry[] symbolTableEntries = new SymbolTableEntry[] { new VolVariable(varName, null) };
                PlotData plotData = getLineScan(ssVolume, image, mesh);
                PlotPane plotPane = new PlotPane();
                DataSymbolMetadataResolver resolver = null;
                Plot2D plot2D = new Plot2D(symbolTableEntries, resolver, new String[] { varName }, new PlotData[] { plotData }, new String[] { "Values along curve", "Distance (\u00b5m)", "[" + varName + "]" });
                plotPane.setPlot2D(plot2D);
                DialogUtils.showComponentCloseDialog(jframe, plotPane, "plot");
            } catch (ImageException | IOException | DataAccessException | MathException e1) {
                e1.printStackTrace();
            }
        }
    });
    GridBagConstraints plotButtonConstraints = new GridBagConstraints();
    plotButtonConstraints.gridx = 0;
    plotButtonConstraints.gridy = 1;
    panel.add(plotButton, plotButtonConstraints);
    jframe.getContentPane().add(panel);
    jframe.setSize(500, 500);
    jframe.addWindowListener(listener);
    jframe.setVisible(true);
    imagePanel.setSourceDataInfo(source);
}

Example 24 with VolVariable

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

the class SimulationWorkspace method getExpectedSizeBytes.

/**
 * Insert the method's description here.
 * Creation date: (6/12/2001 10:09:25 AM)
 * @return boolean
 * @param simulation cbit.vcell.solver.Simulation
 */
private static long getExpectedSizeBytes(SimulationSymbolTable simSymbolTable) {
    Simulation simulation = simSymbolTable.getSimulation();
    long numTimepoints;
    if (simulation.getMathDescription().isNonSpatialStoch()) {
        numTimepoints = getEstimatedNumTimePointsForStoch(simSymbolTable);
    } else {
        numTimepoints = getExpectedNumTimePoints(simulation);
    }
    int x, y, z;
    int numVariables = 0;
    if (simulation.isSpatial()) {
        x = simulation.getMeshSpecification().getSamplingSize().getX();
        y = simulation.getMeshSpecification().getSamplingSize().getY();
        z = simulation.getMeshSpecification().getSamplingSize().getZ();
        // 
        // compute number of volume variables only (they are multiplied by x*y*z)
        // 
        numVariables = 0;
        Enumeration<Variable> variables = simulation.getMathDescription().getVariables();
        while (variables.hasMoreElements()) {
            Variable var = variables.nextElement();
            if (var instanceof VolVariable) {
                numVariables++;
            }
        }
    } else {
        x = 1;
        y = 1;
        z = 1;
        numVariables = 0;
        Enumeration<Variable> variables = simulation.getMathDescription().getVariables();
        while (variables.hasMoreElements()) {
            Variable var = variables.nextElement();
            if ((var instanceof VolVariable) || (var instanceof Function)) {
                numVariables++;
            }
        }
    }
    // approximate, don't compute header size since it's negligible whenever we approach quota size anyway...
    // values are actually written as longs
    long expectedSize = numTimepoints * numVariables * x * y * z * 8;
    return expectedSize;
}

Example 25 with VolVariable

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

the class MatlabOdeFileCoder method write_V6_MFile.

/**
 * Insert the method's description here.
 * Creation date: (3/8/00 10:31:52 PM)
 */
public void write_V6_MFile(java.io.PrintWriter pw, String functionName) throws MathException, ExpressionException {
    MathDescription mathDesc = simulation.getMathDescription();
    if (!mathDesc.isValid()) {
        throw new MathException("invalid math description\n" + mathDesc.getWarning());
    }
    if (mathDesc.isSpatial()) {
        throw new MathException("spatial math description, cannot create ode file");
    }
    if (mathDesc.hasFastSystems()) {
        throw new MathException("math description contains algebraic constraints, cannot create .m file");
    }
    // 
    // print function declaration
    // 
    pw.println("function [T,Y,yinit,param, allNames, allValues] = " + functionName + "(argTimeSpan,argYinit,argParam)");
    pw.println("% [T,Y,yinit,param] = " + functionName + "(argTimeSpan,argYinit,argParam)");
    pw.println("%");
    pw.println("% input:");
    pw.println("%     argTimeSpan is a vector of start and stop times (e.g. timeSpan = [0 10.0])");
    pw.println("%     argYinit is a vector of initial conditions for the state variables (optional)");
    pw.println("%     argParam is a vector of values for the parameters (optional)");
    pw.println("%");
    pw.println("% output:");
    pw.println("%     T is the vector of times");
    pw.println("%     Y is the vector of state variables");
    pw.println("%     yinit is the initial conditions that were used");
    pw.println("%     param is the parameter vector that was used");
    pw.println("%     allNames is the output solution variable names");
    pw.println("%     allValues is the output solution variable values corresponding to the names");
    pw.println("%");
    pw.println("%     example of running this file: [T,Y,yinit,param,allNames,allValues] = myMatlabFunc; <-(your main function name)");
    pw.println("%");
    VariableHash varHash = new VariableHash();
    for (Variable var : simulationSymbolTable.getVariables()) {
        varHash.addVariable(var);
    }
    Variable[] variables = varHash.getTopologicallyReorderedVariables();
    CompartmentSubDomain subDomain = (CompartmentSubDomain) mathDesc.getSubDomains().nextElement();
    // 
    // collect "true" constants (Constants without identifiers)
    // 
    // 
    // collect "variables" (VolVariables only)
    // 
    // 
    // collect "functions" (Functions and Constants with identifiers)
    // 
    Vector<Constant> constantList = new Vector<Constant>();
    Vector<VolVariable> volVarList = new Vector<VolVariable>();
    Vector<Variable> functionList = new Vector<Variable>();
    for (int i = 0; i < variables.length; i++) {
        if (variables[i] instanceof Constant) {
            Constant constant = (Constant) variables[i];
            String[] symbols = constant.getExpression().getSymbols();
            if (symbols == null || symbols.length == 0) {
                constantList.addElement(constant);
            } else {
                functionList.add(constant);
            }
        } else if (variables[i] instanceof VolVariable) {
            volVarList.addElement((VolVariable) variables[i]);
        } else if (variables[i] instanceof Function) {
            functionList.addElement(variables[i]);
        }
    }
    Constant[] constants = (Constant[]) BeanUtils.getArray(constantList, Constant.class);
    VolVariable[] volVars = (VolVariable[]) BeanUtils.getArray(volVarList, VolVariable.class);
    Variable[] functions = (Variable[]) BeanUtils.getArray(functionList, Variable.class);
    int numVars = volVarList.size() + functionList.size();
    String varNamesForStringArray = "";
    String varNamesForValueArray = "";
    for (Variable var : volVarList) {
        varNamesForStringArray = varNamesForStringArray + "'" + var.getName() + "';";
        varNamesForValueArray = varNamesForValueArray + var.getName() + " ";
    }
    for (Variable func : functionList) {
        varNamesForStringArray = varNamesForStringArray + "'" + func.getName() + "';";
        varNamesForValueArray = varNamesForValueArray + func.getName() + " ";
    }
    pw.println("");
    pw.println("%");
    pw.println("% Default time span");
    pw.println("%");
    double beginTime = 0.0;
    double endTime = simulation.getSolverTaskDescription().getTimeBounds().getEndingTime();
    pw.println("timeSpan = [" + beginTime + " " + endTime + "];");
    pw.println("");
    pw.println("% output variable lengh and names");
    pw.println("numVars = " + numVars + ";");
    pw.println("allNames = {" + varNamesForStringArray + "};");
    pw.println("");
    pw.println("if nargin >= 1");
    pw.println("\tif length(argTimeSpan) > 0");
    pw.println("\t\t%");
    pw.println("\t\t% TimeSpan overridden by function arguments");
    pw.println("\t\t%");
    pw.println("\t\ttimeSpan = argTimeSpan;");
    pw.println("\tend");
    pw.println("end");
    pw.println("%");
    pw.println("% Default Initial Conditions");
    pw.println("%");
    pw.println("yinit = [");
    for (int j = 0; j < volVars.length; j++) {
        Expression initial = subDomain.getEquation(volVars[j]).getInitialExpression();
        double defaultInitialCondition = 0;
        try {
            initial.bindExpression(mathDesc);
            defaultInitialCondition = initial.evaluateConstant();
            pw.println("\t" + defaultInitialCondition + ";\t\t% yinit(" + (j + 1) + ") is the initial condition for '" + cbit.vcell.parser.SymbolUtils.getEscapedTokenMatlab(volVars[j].getName()) + "'");
        } catch (ExpressionException e) {
            e.printStackTrace(System.out);
            pw.println("\t" + initial.infix_Matlab() + ";\t\t% yinit(" + (j + 1) + ") is the initial condition for '" + cbit.vcell.parser.SymbolUtils.getEscapedTokenMatlab(volVars[j].getName()) + "'");
        // throw new RuntimeException("error evaluating initial condition for variable "+volVars[j].getName());
        }
    }
    pw.println("];");
    pw.println("if nargin >= 2");
    pw.println("\tif length(argYinit) > 0");
    pw.println("\t\t%");
    pw.println("\t\t% initial conditions overridden by function arguments");
    pw.println("\t\t%");
    pw.println("\t\tyinit = argYinit;");
    pw.println("\tend");
    pw.println("end");
    pw.println("%");
    pw.println("% Default Parameters");
    pw.println("%   constants are only those \"Constants\" from the Math Description that are just floating point numbers (no identifiers)");
    pw.println("%   note: constants of the form \"A_init\" are really initial conditions and are treated in \"yinit\"");
    pw.println("%");
    pw.println("param = [");
    int paramIndex = 0;
    for (int i = 0; i < constants.length; i++) {
        pw.println("\t" + constants[i].getExpression().infix_Matlab() + ";\t\t% param(" + (paramIndex + 1) + ") is '" + cbit.vcell.parser.SymbolUtils.getEscapedTokenMatlab(constants[i].getName()) + "'");
        paramIndex++;
    }
    pw.println("];");
    pw.println("if nargin >= 3");
    pw.println("\tif length(argParam) > 0");
    pw.println("\t\t%");
    pw.println("\t\t% parameter values overridden by function arguments");
    pw.println("\t\t%");
    pw.println("\t\tparam = argParam;");
    pw.println("\tend");
    pw.println("end");
    pw.println("%");
    pw.println("% invoke the integrator");
    pw.println("%");
    pw.println("[T,Y] = ode15s(@f,timeSpan,yinit,odeset('OutputFcn',@odeplot),param,yinit);");
    pw.println("");
    pw.println("% get the solution");
    pw.println("all = zeros(size(T), numVars);");
    pw.println("for i = 1:size(T)");
    pw.println("\tall(i,:) = getRow(T(i), Y(i,:), yinit, param);");
    pw.println("end");
    pw.println("");
    pw.println("allValues = all;");
    pw.println("end");
    // get row data for solution
    pw.println("");
    pw.println("% -------------------------------------------------------");
    pw.println("% get row data");
    pw.println("function rowValue = getRow(t,y,y0,p)");
    // 
    // print volVariables (in order and assign to var vector)
    // 
    pw.println("\t% State Variables");
    for (int i = 0; i < volVars.length; i++) {
        pw.println("\t" + cbit.vcell.parser.SymbolUtils.getEscapedTokenMatlab(volVars[i].getName()) + " = y(" + (i + 1) + ");");
    }
    // 
    // print constants
    // 
    pw.println("\t% Constants");
    paramIndex = 0;
    for (int i = 0; i < constants.length; i++) {
        pw.println("\t" + cbit.vcell.parser.SymbolUtils.getEscapedTokenMatlab(constants[i].getName()) + " = p(" + (paramIndex + 1) + ");");
        paramIndex++;
    }
    // 
    // print variables
    // 
    pw.println("\t% Functions");
    for (int i = 0; i < functions.length; i++) {
        pw.println("\t" + cbit.vcell.parser.SymbolUtils.getEscapedTokenMatlab(functions[i].getName()) + " = " + functions[i].getExpression().infix_Matlab() + ";");
    }
    pw.println("");
    pw.println("\trowValue = [" + varNamesForValueArray + "];");
    pw.println("end");
    // 
    // print ode-rate
    // 
    pw.println("");
    pw.println("% -------------------------------------------------------");
    pw.println("% ode rate");
    pw.println("function dydt = f(t,y,p,y0)");
    // 
    // print volVariables (in order and assign to var vector)
    // 
    pw.println("\t% State Variables");
    for (int i = 0; i < volVars.length; i++) {
        pw.println("\t" + cbit.vcell.parser.SymbolUtils.getEscapedTokenMatlab(volVars[i].getName()) + " = y(" + (i + 1) + ");");
    }
    // 
    // print constants
    // 
    pw.println("\t% Constants");
    paramIndex = 0;
    for (int i = 0; i < constants.length; i++) {
        pw.println("\t" + cbit.vcell.parser.SymbolUtils.getEscapedTokenMatlab(constants[i].getName()) + " = p(" + (paramIndex + 1) + ");");
        paramIndex++;
    }
    // 
    // print variables
    // 
    pw.println("\t% Functions");
    for (int i = 0; i < functions.length; i++) {
        pw.println("\t" + cbit.vcell.parser.SymbolUtils.getEscapedTokenMatlab(functions[i].getName()) + " = " + functions[i].getExpression().infix_Matlab() + ";");
    }
    pw.println("\t% Rates");
    pw.println("\tdydt = [");
    for (int i = 0; i < volVars.length; i++) {
        pw.println("\t\t" + subDomain.getEquation(volVars[i]).getRateExpression().infix_Matlab() + ";    % rate for " + cbit.vcell.parser.SymbolUtils.getEscapedTokenMatlab(volVars[i].getName()));
    }
    pw.println("\t];");
    pw.println("end");
}