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

Example 61 with MathException

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

the class CopasiOptimizationSolver method getOdeSolverResultSet.

private static ODESolverResultSet getOdeSolverResultSet(RowColumnResultSet rcResultSet, SimulationSymbolTable simSymbolTable, String[] parameterNames, double[] parameterValues) {
    // 
    // get simulation results - copy from RowColumnResultSet into OdeSolverResultSet
    // 
    ODESolverResultSet odeSolverResultSet = new ODESolverResultSet();
    for (int i = 0; i < rcResultSet.getDataColumnCount(); i++) {
        odeSolverResultSet.addDataColumn(new ODESolverResultSetColumnDescription(rcResultSet.getColumnDescriptions(i).getName()));
    }
    for (int i = 0; i < rcResultSet.getRowCount(); i++) {
        odeSolverResultSet.addRow(rcResultSet.getRow(i));
    }
    // 
    // add appropriate Function columns to result set
    // 
    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).flatten();
                // 
                for (int j = 0; parameterNames != null && j < parameterNames.length; j++) {
                    exp1.substituteInPlace(new Expression(parameterNames[j]), new Expression(parameterValues[j]));
                }
            } 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);
            }
        }
    }
    return odeSolverResultSet;
}

Example 62 with MathException

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

the class GibsonSolver method getStochSolverResultSet.

/**
 * Get data columns and function columns to be ready for plot.
 * Creation date: (8/15/2006 11:36:43 AM)
 * @return cbit.vcell.solver.stoch.StochSolverResultSet
 */
public ODESolverResultSet getStochSolverResultSet() {
    // read .stoch file, this funciton here equals to getODESolverRestultSet()+getStateVariableResultSet()  in ODE.
    ODESolverResultSet stSolverResultSet = new ODESolverResultSet();
    FileInputStream inputStream = null;
    try {
        inputStream = new FileInputStream(getBaseName() + IDA_DATA_EXTENSION);
        InputStreamReader inputStreamReader = new InputStreamReader(inputStream);
        BufferedReader bufferedReader = new BufferedReader(inputStreamReader);
        // Read header
        String line = bufferedReader.readLine();
        if (line == null) {
            // throw exception
            System.out.println("There is no data in output file!");
            return null;
        }
        while (line.indexOf(':') > 0) {
            String name = line.substring(0, line.indexOf(':'));
            stSolverResultSet.addDataColumn(new ODESolverResultSetColumnDescription(name));
            line = line.substring(line.indexOf(':') + 1);
        }
        // Read data
        while ((line = bufferedReader.readLine()) != null) {
            line = line + "\t";
            double[] values = new double[stSolverResultSet.getDataColumnCount()];
            boolean bCompleteRow = true;
            for (int i = 0; i < stSolverResultSet.getDataColumnCount(); i++) {
                if (line.indexOf('\t') == -1) {
                    bCompleteRow = false;
                    break;
                } else {
                    String value = line.substring(0, line.indexOf('\t')).trim();
                    values[i] = Double.valueOf(value).doubleValue();
                    line = line.substring(line.indexOf('\t') + 1);
                }
            }
            if (bCompleteRow) {
                stSolverResultSet.addRow(values);
            } else {
                break;
            }
        }
    // 
    } catch (Exception e) {
        e.printStackTrace(System.out);
        return null;
    } finally {
        try {
            if (inputStream != null) {
                inputStream.close();
            }
        } catch (Exception ex) {
            lg.error(ex.getMessage(), ex);
        }
    }
    /*
	Add appropriate Function columns to result set if the stochastic simulation is to display the trajectory.
	No function columns for the results of multiple stochastic trials
	*/
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    if (simSymbolTable.getSimulation().getSolverTaskDescription().getStochOpt().getNumOfTrials() == 1) {
        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);
                    stSolverResultSet.addFunctionColumn(cd);
                } catch (ExpressionException e) {
                    e.printStackTrace(System.out);
                }
            }
        }
    }
    return stSolverResultSet;
}

Example 63 with MathException

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

the class HybridSolver method getHybridSolverResultSet.

/**
 * Get data columns and function columns to be ready for plot.
 * Creation date: (8/15/2006 11:36:43 AM)
 * @return cbit.vcell.solver.stoch.StochSolverResultSet
 */
private ODESolverResultSet getHybridSolverResultSet() {
    // read .stoch file, this funciton here equals to getODESolverRestultSet()+getStateVariableResultSet()  in ODE.
    ODESolverResultSet stSolverResultSet = new ODESolverResultSet();
    try {
        String filename = getBaseName() + NETCDF_DATA_EXTENSION;
        NetCDFEvaluator ncEva = new NetCDFEvaluator();
        NetCDFReader ncReader = null;
        try {
            ncEva.setNetCDFTarget(filename);
            ncReader = ncEva.getNetCDFReader();
        } catch (Exception e) {
            e.printStackTrace(System.err);
            throw new RuntimeException("Cannot open simulation result file: " + filename + "!");
        }
        // Read result according to trial number
        if (ncReader.getNumTrials() == 1) {
            // Read header
            String[] varNames = ncReader.getSpeciesNames_val();
            // first column will be time t.
            stSolverResultSet.addDataColumn(new ODESolverResultSetColumnDescription("t"));
            // following columns are stoch variables
            for (int i = 0; i < varNames.length; i++) {
                stSolverResultSet.addDataColumn(new ODESolverResultSetColumnDescription(varNames[i]));
            }
            // Read data
            // data only, no time points
            ArrayDouble data = (ArrayDouble) ncEva.getTimeSeriesData(1);
            double[] timePoints = ncReader.getTimePoints();
            System.out.println("time points length is " + timePoints.length);
            // shape[0]:num of timepoints, shape[1]: num of species
            int[] shape = data.getShape();
            if (// one species
            shape.length == 1) {
                ArrayDouble.D1 temData = (ArrayDouble.D1) data;
                System.out.println("one species in time series data and size is " + temData.getSize());
                for (// rows
                int k = 0; // rows
                k < timePoints.length; // rows
                k++) {
                    double[] values = new double[stSolverResultSet.getDataColumnCount()];
                    values[0] = timePoints[k];
                    for (int i = 1; i < stSolverResultSet.getDataColumnCount(); i++) {
                        values[i] = temData.get(k);
                    }
                    stSolverResultSet.addRow(values);
                }
            }
            if (// more than one species
            shape.length == 2) {
                ArrayDouble.D2 temData = (ArrayDouble.D2) data;
                System.out.println("multiple species in time series, the length of time series is :" + data.getShape()[0] + ", and the total number of speceis is: " + data.getShape()[1]);
                for (// rows
                int k = 0; // rows
                k < timePoints.length; // rows
                k++) {
                    double[] values = new double[stSolverResultSet.getDataColumnCount()];
                    values[0] = timePoints[k];
                    for (int i = 1; i < stSolverResultSet.getDataColumnCount(); i++) {
                        values[i] = temData.get(k, i - 1);
                    }
                    stSolverResultSet.addRow(values);
                }
            }
        } else if (ncReader.getNumTrials() > 1) {
            // Read header
            String[] varNames = ncReader.getSpeciesNames_val();
            // first column will be time t.
            stSolverResultSet.addDataColumn(new ODESolverResultSetColumnDescription("TrialNo"));
            // following columns are stoch variables
            for (int i = 0; i < varNames.length; i++) {
                stSolverResultSet.addDataColumn(new ODESolverResultSetColumnDescription(varNames[i]));
            }
            // Read data
            // data only, no trial numbers
            ArrayDouble data = (ArrayDouble) ncEva.getDataOverTrials(ncReader.getTimePoints().length - 1);
            int[] trialNum = ncEva.getNetCDFReader().getTrialNumbers();
            // System.out.println("total trials are "+trialNum.length);
            // shape[0]:number of trials, shape[1]: num of species
            int[] shape = data.getShape();
            if (// one species
            shape.length == 1) {
                ArrayDouble.D1 temData = (ArrayDouble.D1) data;
                // System.out.println("one species over trials, size is: "+temData.getSize());
                for (// rows
                int k = 0; // rows
                k < trialNum.length; // rows
                k++) {
                    double[] values = new double[stSolverResultSet.getDataColumnCount()];
                    values[0] = trialNum[k];
                    for (int i = 1; i < stSolverResultSet.getDataColumnCount(); i++) {
                        values[i] = temData.get(k);
                    }
                    stSolverResultSet.addRow(values);
                }
            }
            if (// more than one species
            shape.length == 2) {
                ArrayDouble.D2 temData = (ArrayDouble.D2) data;
                // System.out.println("multiple species in multiple trials, the length of trials is :"+data.getShape()[0]+", and the total number of speceis is: "+data.getShape()[1]);
                for (// rows
                int k = 0; // rows
                k < trialNum.length; // rows
                k++) {
                    double[] values = new double[stSolverResultSet.getDataColumnCount()];
                    values[0] = trialNum[k];
                    for (int i = 1; i < stSolverResultSet.getDataColumnCount(); i++) {
                        values[i] = temData.get(k, i - 1);
                    }
                    stSolverResultSet.addRow(values);
                }
            }
        } else {
            throw new RuntimeException("Number of trials should be a countable positive value, from 1 to N.");
        }
    } catch (Exception e) {
        e.printStackTrace(System.err);
        throw new RuntimeException("Problem encountered in parsing hybrid simulation results.\n" + e.getMessage());
    }
    /*
	 *Add appropriate Function columns to result set if the stochastic simulation is to display the trajectory.
	 *No function columns for the results of multiple stochastic trials.
	 *In stochastic simulation the functions include probability functions and clamped variable.
	 */
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    if (simSymbolTable.getSimulation().getSolverTaskDescription().getStochOpt().getNumOfTrials() == 1) {
        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);
                    stSolverResultSet.addFunctionColumn(cd);
                } catch (ExpressionException e) {
                    e.printStackTrace(System.out);
                }
            }
        }
    }
    return stSolverResultSet;
}

Example 64 with MathException

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

the class PropensitySolver method solvePropensity.

public static PropensityFunction solvePropensity(Expression original_expression, String[] speciesNames, int maxOrder) throws ExpressionException, MathException {
    boolean[][] rootExists = new boolean[speciesNames.length][maxOrder];
    ArrayList<Binomial> binomialList = new ArrayList<Binomial>();
    int[] order = new int[speciesNames.length];
    int totalOrder = 0;
    // Actually, if order is 3, the roots can be either 0 (implies var^3) or 0,1,2 (implies var*(var-1)*(var-2))
    for (int i = 0; i < speciesNames.length; i++) {
        order[i] = -1;
        Expression diff_i = new Expression(original_expression);
        for (int j = 0; j < rootExists[0].length; j++) {
            Expression exp_subX_j = original_expression.getSubstitutedExpression(new Expression(speciesNames[i]), new Expression(j));
            if (ExpressionUtils.functionallyEquivalent(exp_subX_j, new Expression(0.0), false, 1e-8, 1e-8)) {
                rootExists[i][j] = true;
            } else {
                rootExists[i][j] = false;
            }
            diff_i = diff_i.differentiate(speciesNames[i]).flatten();
            if (!ExpressionUtils.functionallyEquivalent(diff_i, new Expression(0.0), false, 1e-8, 1e-8)) {
                order[i] = j + 1;
                if (order[i] > maxOrder) {
                    throw new MathException("order of " + speciesNames[i] + " in expression " + original_expression.infix() + " is larger than expected for a propensity function");
                }
            }
        }
    }
    // check if the roots are valid for mass actions (roots should start with 0 and increase consecutively)
    for (int i = 0; i < speciesNames.length; i++) {
        int ord = order[i];
        totalOrder = totalOrder + order[i];
        int expectedRootCount = 0;
        int totalRootCount = 0;
        // expectedRootCount should be 1 or order
        for (int j = 0; j < ord; j++) {
            if (rootExists[i][j]) {
                expectedRootCount++;
            }
        }
        // total root count
        for (int j = 0; j < rootExists[i].length; j++) {
            if (rootExists[i][j]) {
                totalRootCount++;
            }
        }
        // either root is 0, or number of roots equals to order.
        if (!((totalRootCount == 1 && expectedRootCount == totalRootCount && rootExists[i][0]) || (totalRootCount > 1 && expectedRootCount == totalRootCount))) {
            throw new MathException("\n\nSpecies \'" + speciesNames[i] + "\' has wrong form in propensity function. Hybrid solvers require strict propensity functions which should only contain the reactants in the function. \n\nIntegers in '" + speciesNames[i] + "' binomials should start with 0 and increase consecutively." + "e.g. species1*(species1-1)*(species1-2)*species2*(speceis2-1)...");
        }
    }
    // get all the constraints(for expponets vars)
    ArrayList<String> constraintSymbols = new ArrayList<String>();
    ArrayList<Expression> constraints = new ArrayList<Expression>();
    ArrayList<ArrayList<String>> exponents = new ArrayList<ArrayList<String>>();
    Expression canonical_expression = new Expression(1.0);
    for (int i = 0; i < rootExists.length; i++) {
        exponents.add(new ArrayList<String>());
        for (int j = 0; j < rootExists[i].length; j++) {
            if (rootExists[i][j]) {
                String exponentVar = "n_" + speciesNames[i] + "_" + j;
                constraintSymbols.add(exponentVar);
                exponents.get(i).add(exponentVar);
                canonical_expression = Expression.mult(canonical_expression, new Expression("(" + speciesNames[i] + " - " + j + ")^" + exponentVar));
                Binomial binomial = new Binomial();
                binomial.varName = speciesNames[i];
                binomial.root = j;
                binomial.powerVarName = exponentVar;
                binomialList.add(binomial);
            }
        }
        Expression expConstraint = new Expression(0.0);
        for (int j = 0; j < exponents.get(i).size(); j++) {
            expConstraint = Expression.add(expConstraint, new Expression(exponents.get(i).get(j)));
            constraints.add(new Expression(exponents.get(i).get(j) + " >= 1"));
            constraints.add(new Expression(exponents.get(i).get(j) + " <= " + order[i]));
        }
        constraints.add(new Expression(expConstraint.flatten().infix() + " == " + order[i]));
    }
    SimpleSymbolTable constraintSymbolTable = new SimpleSymbolTable(constraintSymbols.toArray(new String[constraintSymbols.size()]));
    for (int i = 0; i < constraints.size(); i++) {
        constraints.get(i).bindExpression(constraintSymbolTable);
    }
    // the array list to save different conbination values of orders of binomials. one element saves one conbination
    ArrayList<int[]> orderList = new ArrayList<int[]>();
    try {
        PropensitySolver.RootOrderIterator iter = new PropensitySolver.RootOrderIterator(totalOrder, constraintSymbols.size());
        while (iter.nextOrder() != null) {
            orderList.add(iter.getCurrentOrders().clone());
        }
    } catch (Exception e) {
        e.printStackTrace();
        throw new MathException("Did not recognize propensity function from \"" + original_expression.infix() + " product of species shoulde have order less than or equal to 3 and the form like 'rateConstant*species1*(species1-1)*species2'.");
    }
    Expression expression_for_K = Expression.mult(original_expression.flatten(), Expression.invert(canonical_expression.flatten()));
    Random rand = new Random(0);
    for (int[] values : orderList) {
        // 
        // for each valid combination of orders (via constraints), choose the one that always gives the same answer for K
        // 
        Expression order_substituted_canonical_expression = new Expression(canonical_expression);
        Expression order_substituted_expression_for_K = new Expression(expression_for_K);
        for (int j = 0; j < values.length; j++) {
            order_substituted_canonical_expression.substituteInPlace(new Expression(constraintSymbols.get(j)), new Expression(values[j]));
            order_substituted_expression_for_K.substituteInPlace(new Expression(constraintSymbols.get(j)), new Expression(values[j]));
        }
        order_substituted_canonical_expression = order_substituted_canonical_expression.flatten();
        order_substituted_expression_for_K = order_substituted_expression_for_K.flatten();
        Hashtable<String, Integer> k_hash = new Hashtable<String, Integer>();
        boolean bFailed = false;
        // for the specific conbination of orders, replace species names with random numbers, therefore to get rate constant value.
        for (int i = 0; !bFailed && i < 40; i++) {
            Expression substituted_expression_for_K = new Expression(order_substituted_expression_for_K);
            // substitute variables with random number and save the results in a hashtable.
            for (int j = 0; j < speciesNames.length; j++) {
                int randInt = rand.nextInt(20) - 10;
                while (randInt < rootExists[j].length && randInt >= 0 && rootExists[j][randInt]) {
                    // make sure it is not a root of the system
                    randInt = rand.nextInt(20) - 10;
                }
                substituted_expression_for_K.substituteInPlace(new Expression(speciesNames[j]), new Expression(randInt));
                substituted_expression_for_K = substituted_expression_for_K.flatten();
            }
            String k_key = substituted_expression_for_K.infix();
            if (k_hash.containsKey(k_key)) {
                int numOccur = k_hash.get(k_key);
                k_hash.put(k_key, numOccur + 1);
            } else {
                k_hash.put(k_key, 1);
            }
        }
        // to see if the results in the hashtable are the same all the time. if so, the conbination is correct.
        Set<String> keySet = k_hash.keySet();
        String[] keys = keySet.toArray(new String[k_hash.size()]);
        Expression k_exp_0 = new Expression(keys[0]);
        Expression k_most_popular = k_exp_0;
        Integer mostOccurances = k_hash.get(keys[0]);
        // why not check if there is only one key, the conbination is correct??
        for (int i = 1; !bFailed && i < keys.length; i++) {
            Expression k_exp_i = new Expression(keys[i]);
            Integer occurances = k_hash.get(keys[i]);
            if (occurances > mostOccurances) {
                mostOccurances = occurances;
                k_most_popular = k_exp_i;
            }
            if (!ExpressionUtils.functionallyEquivalent(k_exp_0, k_exp_i, false, 1e-8, 1e-8)) {
                bFailed = true;
            }
        }
        if (!bFailed) {
            PropensityFunction propensityFunction = new PropensityFunction();
            propensityFunction.speciesNames = speciesNames.clone();
            propensityFunction.speciesOrders = order.clone();
            propensityFunction.canonicalExpression = Expression.mult(k_most_popular, order_substituted_canonical_expression);
            propensityFunction.k_expression = k_most_popular;
            for (int i = 0; i < binomialList.size(); i++) {
                binomialList.get(i).power = (int) values[i];
            }
            propensityFunction.binomials = binomialList.toArray(new Binomial[binomialList.size()]);
            return propensityFunction;
        }
    }
    throw new MathException("Did not recognize propensity function from \"" + original_expression.infix() + "\", looking for product of species (e.g. \"rateConstant*species1*(species1-1)*species2\")");
}

Example 65 with MathException

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

the class XmlReader method getVarIniPoissonExpectedCount.

private VarIniCondition getVarIniPoissonExpectedCount(Element param, MathDescription md) throws XmlParseException, MathException, ExpressionException {
    // retrieve values
    Expression exp = unMangleExpression(param.getText());
    String name = unMangle(param.getAttributeValue(XMLTags.NameAttrTag));
    Variable var = md.getVariable(name);
    if (var == null) {
        throw new MathFormatException("variable " + name + " not defined");
    }
    if (!(var instanceof StochVolVariable)) {
        throw new MathFormatException("variable " + name + " not a Stochastic Volume Variable");
    }
    try {
        VarIniCondition varIni = new VarIniPoissonExpectedCount(var, exp);
        return varIni;
    } catch (Exception e) {
        e.printStackTrace();
    }
    return null;
}