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

Example 31 with MathException

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

the class DefaultODESolver method createODESolverResultSet.

/**
 */
private ODESolverResultSet createODESolverResultSet() throws ExpressionException {
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    // 
    // create symbol table for binding expression
    // 
    String[] symbols = new String[fieldIdentifiers.size()];
    for (int i = 0; i < symbols.length; i++) {
        symbols[i] = ((Variable) fieldIdentifiers.elementAt(i)).getName();
    }
    // Initialize the ResultSet...
    ODESolverResultSet odeSolverResultSet = new ODESolverResultSet();
    odeSolverResultSet.addDataColumn(new ODESolverResultSetColumnDescription(ReservedVariable.TIME.getName()));
    for (int i = 0; i < getStateVariableCount(); i++) {
        StateVariable stateVariable = getStateVariable(i);
        if (stateVariable instanceof SensStateVariable) {
            SensStateVariable sensStateVariable = (SensStateVariable) stateVariable;
            odeSolverResultSet.addDataColumn(new ODESolverResultSetColumnDescription(sensStateVariable.getVariable().getName(), sensStateVariable.getParameter().getName(), sensStateVariable.getVariable().getName()));
        } else {
            odeSolverResultSet.addDataColumn(new ODESolverResultSetColumnDescription(stateVariable.getVariable().getName()));
        }
    }
    Variable[] variables = simSymbolTable.getVariables();
    for (int i = 0; i < variables.length; i++) {
        if (variables[i] instanceof Function && SimulationSymbolTable.isFunctionSaved((Function) variables[i])) {
            Function function = (Function) variables[i];
            Expression exp1 = new Expression(function.getExpression());
            try {
                exp1 = simSymbolTable.substituteFunctions(exp1);
            } catch (MathException e) {
                e.printStackTrace(System.out);
                throw new RuntimeException("Substitute function failed on function " + function.getName() + " " + e.getMessage());
            }
            odeSolverResultSet.addFunctionColumn(new FunctionColumnDescription(exp1.flatten(), function.getName(), null, function.getName(), false));
        }
    }
    // 
    if (getSensitivityParameter() != null) {
        if (odeSolverResultSet.findColumn(getSensitivityParameter().getName()) == -1) {
            FunctionColumnDescription fcd = new FunctionColumnDescription(new Expression(getSensitivityParameter().getConstantValue()), getSensitivityParameter().getName(), null, getSensitivityParameter().getName(), false);
            odeSolverResultSet.addFunctionColumn(fcd);
        }
        StateVariable[] stateVars = (StateVariable[]) org.vcell.util.BeanUtils.getArray(fieldStateVariables, StateVariable.class);
        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());
                try {
                    depSensFnExpr = simSymbolTable.substituteFunctions(depSensFnExpr);
                } catch (MathException e) {
                    e.printStackTrace(System.out);
                    throw new RuntimeException("Substitute function failed on function " + depSensFunction.getName() + " " + e.getMessage());
                }
                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);
                }
            }
        }
    }
    return (odeSolverResultSet);
}

Example 32 with MathException

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

the class DefaultODESolver method initialize.

/**
 * This method was created by a SmartGuide.
 * @exception SolverException The exception description.
 */
protected void initialize() throws SolverException {
    SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
    Simulation sim = simSymbolTable.getSimulation();
    try {
        // create a fast system if necessary
        fieldFastAlgebraicSystem = null;
        if (getSubDomain().getFastSystem() != null) {
            fieldFastAlgebraicSystem = new FastAlgebraicSystem(new FastSystemAnalyzer(getSubDomain().getFastSystem(), simSymbolTable));
        }
        // refreshIdentifiers();
        fieldIdentifiers = createIdentifiers();
        fieldSensVariables = createSensitivityVariables();
        // refreshStateVariables();
        fieldStateVariables = createStateVariables();
        // 
        // allocate ValueVectors object
        fieldValueVectors = new ValueVectors(getValueVectorCount(), fieldIdentifiers.size());
        // initialize indexes of variables
        fieldVariableIndexes = new int[getStateVariableCount()];
        for (int i = 0; i < getStateVariableCount(); i++) {
            fieldVariableIndexes[i] = getStateVariable(i).getVariable().getIndex();
        }
        // initialize constants
        double[] initialValues = getValueVector(0);
        for (int i = 0; i < fieldIdentifiers.size(); i++) {
            if (fieldIdentifiers.elementAt(i) instanceof Constant) {
                Constant constant = (Constant) fieldIdentifiers.elementAt(i);
                constant.bind(simSymbolTable);
                if (constant.isConstant()) {
                    // constant.getValue();
                    initialValues[constant.getIndex()] = constant.getExpression().evaluateConstant();
                } else {
                    throw new SolverException("cannot evaluate constant '" + constant.getName() + "' = " + constant.getExpression());
                }
            }
        }
        // initialize variables
        for (int i = 0; i < getStateVariableCount(); i++) {
            initialValues[getVariableIndex(i)] = getStateVariable(i).evaluateIC(initialValues);
        }
        fieldODESolverResultSet = createODESolverResultSet();
        // reset - in the ** default ** solvers we don't pick up from where we left off, we can override that behaviour in integrate() if ever necessary
        fieldCurrentTime = sim.getSolverTaskDescription().getTimeBounds().getStartingTime();
    } catch (ExpressionException expressionException) {
        expressionException.printStackTrace(System.out);
        throw new SolverException(expressionException.getMessage());
    } catch (MathException mathException) {
        mathException.printStackTrace(System.out);
        throw new SolverException(mathException.getMessage());
    }
}

Example 33 with MathException

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

the class IDAFileWriter method writeEquations.

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

Example 34 with MathException

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

the class Jacobian method parseMathDesc.

/**
 * This method was created by a SmartGuide.
 * @exception java.lang.Exception The exception description.
 */
private void parseMathDesc() throws MathException {
    Vector equationList = new Vector();
    Enumeration enum1 = subDomain.getEquations();
    while (enum1.hasMoreElements()) {
        Equation equ = (Equation) enum1.nextElement();
        if (equ instanceof OdeEquation) {
            equationList.addElement(equ);
        } else {
            throw new MathException("encountered non-ode equation, unsupported");
        }
    }
    Vector variableList = new Vector();
    enum1 = mathDesc.getVariables();
    while (enum1.hasMoreElements()) {
        Variable var = (Variable) enum1.nextElement();
        if (var instanceof VolVariable) {
            variableList.addElement(var);
        }
    }
    if (equationList.size() != variableList.size()) {
        throw new MathException("there are " + equationList.size() + " equations and " + variableList.size() + " variables");
    }
    numVariables = variableList.size();
    rates = new Expression[numVariables];
    vars = new VolVariable[numVariables];
    for (int i = 0; i < numVariables; i++) {
        OdeEquation odeEqu = (OdeEquation) equationList.elementAt(i);
        rates[i] = odeEqu.getRateExpression();
        vars[i] = (VolVariable) odeEqu.getVariable();
    }
}

Example 35 with MathException

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

the class SimulationSymbolTable method createAnnotatedFunctionsList.

public Vector<AnnotatedFunction> createAnnotatedFunctionsList(MathDescription mathDescription) throws InconsistentDomainException {
    // Get the list of (volVariables) in the simulation. Needed to determine 'type' of  functions
    boolean bSpatial = getSimulation().isSpatial();
    String[] variableNames = null;
    VariableType[] variableTypes = null;
    if (bSpatial) {
        Variable[] allVariables = getVariables();
        Vector<Variable> varVector = new Vector<Variable>();
        for (int i = 0; i < allVariables.length; i++) {
            if ((allVariables[i] instanceof VolVariable) || (allVariables[i] instanceof VolumeRegionVariable) || (allVariables[i] instanceof MemVariable) || (allVariables[i] instanceof MembraneRegionVariable) || (allVariables[i] instanceof FilamentVariable) || (allVariables[i] instanceof FilamentRegionVariable) || (allVariables[i] instanceof PointVariable) || (allVariables[i] instanceof ParticleVariable) || (allVariables[i] instanceof InsideVariable) || (allVariables[i] instanceof OutsideVariable)) {
                varVector.addElement(allVariables[i]);
            } else if (allVariables[i] instanceof Constant || (allVariables[i] instanceof Function)) {
            } else {
                System.err.println("SimulationSymbolTable.createAnnotatedFunctionsList() found unexpected variable type " + allVariables[i].getClass().getSimpleName() + " in spatial simulation");
            }
        }
        variableNames = new String[varVector.size()];
        for (int i = 0; i < variableNames.length; i++) {
            variableNames[i] = varVector.get(i).getName();
        }
        // Lookup table for variableType for each variable in 'variables' array.
        variableTypes = new VariableType[variableNames.length];
        for (int i = 0; i < variableNames.length; i++) {
            variableTypes[i] = VariableType.getVariableType(varVector.get(i));
        }
    }
    // 
    // Bind and substitute functions to simulation before storing them in the '.functions' file
    // 
    Function[] functions = getFunctions();
    Vector<AnnotatedFunction> annotatedFunctionVector = new Vector<AnnotatedFunction>();
    for (int i = 0; i < functions.length; i++) {
        if (isFunctionSaved(functions[i])) {
            String errString = "";
            VariableType funcType = null;
            try {
                Expression substitutedExp = substituteFunctions(functions[i].getExpression());
                substitutedExp.bindExpression(this);
                functions[i].setExpression(substitutedExp.flatten());
            } catch (MathException e) {
                e.printStackTrace(System.out);
                errString = errString + ", " + e.getMessage();
            // throw new RuntimeException(e.getMessage());
            } catch (ExpressionException e) {
                e.printStackTrace(System.out);
                errString = errString + ", " + e.getMessage();
            // throw new RuntimeException(e.getMessage());
            }
            // 
            // get function's data type from the types of it's identifiers
            // 
            funcType = bSpatial ? getFunctionVariableType(functions[i], mathDescription, variableNames, variableTypes, bSpatial) : VariableType.NONSPATIAL;
            AnnotatedFunction annotatedFunc = new AnnotatedFunction(functions[i].getName(), functions[i].getExpression(), functions[i].getDomain(), errString, funcType, FunctionCategory.PREDEFINED);
            annotatedFunctionVector.addElement(annotatedFunc);
        }
    }
    return annotatedFunctionVector;
}