Examples with DefaultOutputTimeSpec cbit.vcell.solver.DefaultOutputTimeSpec used on opensource projects

Example 21 with DefaultOutputTimeSpec

use of cbit.vcell.solver.DefaultOutputTimeSpec in project vcell by virtualcell.

the class DefaultODESolver method integrate.

/**
 * This method was created by a SmartGuide.
 */
protected void integrate() throws SolverException, UserStopException, IOException {
    try {
        SolverTaskDescription taskDescription = simTask.getSimulation().getSolverTaskDescription();
        double timeStep = taskDescription.getTimeStep().getDefaultTimeStep();
        fieldCurrentTime = taskDescription.getTimeBounds().getStartingTime();
        // before computation begins, settle fast equilibrium
        if (getFastAlgebraicSystem() != null) {
            fieldValueVectors.copyValues(0, 1);
            getFastAlgebraicSystem().initVars(getValueVector(0), getValueVector(1));
            getFastAlgebraicSystem().solveSystem(getValueVector(0), getValueVector(1));
            fieldValueVectors.copyValues(1, 0);
        }
        // check for failure
        check(getValueVector(0));
        updateResultSet();
        // 
        int iteration = 0;
        while (fieldCurrentTime < taskDescription.getTimeBounds().getEndingTime()) {
            checkForUserStop();
            step(fieldCurrentTime, timeStep);
            // update (old = new)
            fieldValueVectors.copyValuesDown();
            // compute fast system
            if (getFastAlgebraicSystem() != null) {
                fieldValueVectors.copyValues(0, 1);
                getFastAlgebraicSystem().initVars(getValueVector(0), getValueVector(1));
                getFastAlgebraicSystem().solveSystem(getValueVector(0), getValueVector(1));
                fieldValueVectors.copyValues(1, 0);
            }
            // check for failure
            check(getValueVector(0));
            // fieldCurrentTime += timeStep;
            iteration++;
            fieldCurrentTime = taskDescription.getTimeBounds().getStartingTime() + iteration * timeStep;
            // Print results if it coincides with a save interval...
            if (taskDescription.getOutputTimeSpec().isDefault()) {
                int keepEvery = ((DefaultOutputTimeSpec) taskDescription.getOutputTimeSpec()).getKeepEvery();
                if ((iteration % keepEvery) == 0) {
                    updateResultSet();
                }
            }
        }
        // store last time point
        if (taskDescription.getOutputTimeSpec().isDefault()) {
            int keepEvery = ((DefaultOutputTimeSpec) taskDescription.getOutputTimeSpec()).getKeepEvery();
            if ((iteration % keepEvery) != 0)
                updateResultSet();
        }
    } catch (ExpressionException | MathException e) {
        throw new SolverException("Solver failed: " + e.getMessage(), e);
    }
}

Example 22 with DefaultOutputTimeSpec

use of cbit.vcell.solver.DefaultOutputTimeSpec in project vcell by virtualcell.

the class OdeFileWriter method write.

/**
 * Insert the method's description here.
 * Creation date: (3/8/00 10:31:52 PM)
 */
public void write(String[] parameterNames) throws Exception {
    createStateVariables();
    createSymbolTable();
    Simulation simulation = simTask.getSimulation();
    if (simulation.getSolverTaskDescription().getUseSymbolicJacobian()) {
        throw new RuntimeException("symbolic jacobian option not yet supported in interpreted Stiff solver");
    }
    writeJMSParamters();
    SolverTaskDescription solverTaskDescription = simulation.getSolverTaskDescription();
    TimeBounds timeBounds = solverTaskDescription.getTimeBounds();
    ErrorTolerance errorTolerance = solverTaskDescription.getErrorTolerance();
    printWriter.println("SOLVER " + getSolverName());
    printWriter.println("STARTING_TIME " + timeBounds.getStartingTime());
    printWriter.println("ENDING_TIME " + timeBounds.getEndingTime());
    printWriter.println("RELATIVE_TOLERANCE " + errorTolerance.getRelativeErrorTolerance());
    printWriter.println("ABSOLUTE_TOLERANCE " + errorTolerance.getAbsoluteErrorTolerance());
    printWriter.println("MAX_TIME_STEP " + simulation.getSolverTaskDescription().getTimeStep().getMaximumTimeStep());
    OutputTimeSpec ots = simulation.getSolverTaskDescription().getOutputTimeSpec();
    if (ots.isDefault()) {
        printWriter.println("KEEP_EVERY " + ((DefaultOutputTimeSpec) ots).getKeepEvery());
    } else if (ots.isUniform()) {
        printWriter.println("OUTPUT_TIME_STEP " + ((UniformOutputTimeSpec) ots).getOutputTimeStep());
    } else if (ots.isExplicit()) {
        printWriter.println("OUTPUT_TIMES " + ((ExplicitOutputTimeSpec) ots).getNumTimePoints());
        printWriter.println(((ExplicitOutputTimeSpec) ots).toSpaceSeperatedMultiLinesOfString());
    }
    if (parameterNames != null && parameterNames.length != 0) {
        printWriter.println("NUM_PARAMETERS " + parameterNames.length);
        for (int i = 0; i < parameterNames.length; i++) {
            printWriter.println(parameterNames[i]);
        }
    }
    HashMap<Discontinuity, String> discontinuityNameMap = new HashMap<Discontinuity, String>();
    String eventString = null;
    if (simulation.getMathDescription().hasEvents()) {
        eventString = writeEvents(discontinuityNameMap);
    }
    String equationString = writeEquations(discontinuityNameMap);
    if (discontinuityNameMap.size() > 0) {
        printWriter.println("DISCONTINUITIES " + discontinuityNameMap.size());
        for (Discontinuity od : discontinuityNameMap.keySet()) {
            printWriter.println(discontinuityNameMap.get(od) + " " + od.getDiscontinuityExp().flatten().infix() + "; " + od.getRootFindingExp().flatten().infix() + ";");
        }
    }
    if (eventString != null) {
        printWriter.print(eventString);
    }
    printWriter.println("NUM_EQUATIONS " + getStateVariableCount());
    printWriter.println(equationString);
}

Example 23 with DefaultOutputTimeSpec

use of cbit.vcell.solver.DefaultOutputTimeSpec in project vcell by virtualcell.

the class RungeKuttaFehlbergSolver method integrate.

protected void integrate() throws cbit.vcell.solver.SolverException, UserStopException, IOException {
    try {
        // Get machine epsilon
        final double DBL_EPSILON = 1.0E-12;
        final double epsilon = DBL_EPSILON;
        final double twentySixEpsilon = 26 * epsilon;
        // 
        SolverTaskDescription taskDescription = simTask.getSimulation().getSolverTaskDescription();
        double startingTime = taskDescription.getTimeBounds().getStartingTime();
        double endingTime = taskDescription.getTimeBounds().getEndingTime();
        double relativeErrorTolerance = taskDescription.getErrorTolerance().getRelativeErrorTolerance();
        double absoluteErrorTolerance = taskDescription.getErrorTolerance().getAbsoluteErrorTolerance();
        fieldCurrentTime = taskDescription.getTimeBounds().getStartingTime();
        // 
        double[] oldValues = getValueVector(0);
        double[] newValues = getValueVector(1);
        // The minimumRelativeError is the minimum acceptable value of
        // relativeError.  Attempts to obtain higher accuracy with this
        // subroutine are usually very expensive and often unsuccessful.
        final double minimumRelativeError = 1e-12;
        // Check input parameters...
        if (relativeErrorTolerance < 0.0 || absoluteErrorTolerance < 0.0 || startingTime >= endingTime) {
            throw new SolverException("Invalid parameters");
        }
        // difficulties arising from impossible accuracy requests
        if (relativeErrorTolerance < 2 * epsilon + minimumRelativeError) {
            // Or just set relativeError = 2 * epsilon + minimumRelativeError???
            throw new SolverException("Relative error too small");
        }
        // before computation begins, settle fast equilibrium
        if (getFastAlgebraicSystem() != null) {
            fieldValueVectors.copyValues(0, 1);
            getFastAlgebraicSystem().initVars(getValueVector(0), getValueVector(1));
            getFastAlgebraicSystem().solveSystem(getValueVector(0), getValueVector(1));
            fieldValueVectors.copyValues(1, 0);
        }
        // check for failure
        check(getValueVector(0));
        // 
        double timeRemaining = endingTime - fieldCurrentTime;
        double tolerance = 0.0;
        double maximumTolerance = 0.0;
        final double maximumTimeStep = taskDescription.getTimeStep().getMaximumTimeStep();
        double h = Math.min(Math.abs(timeRemaining), taskDescription.getTimeStep().getMaximumTimeStep());
        for (int i = 0; i < getStateVariableCount(); i++) {
            int I = getVariableIndex(i);
            tolerance = relativeErrorTolerance * Math.abs(oldValues[I]) + absoluteErrorTolerance;
            if (tolerance > 0.0) {
                double yp = Math.abs(evaluate(oldValues, i));
                if (yp * Math.pow(h, 5.0) > tolerance) {
                    h = Math.pow(tolerance / yp, 0.2);
                }
            }
            maximumTolerance = Math.max(maximumTolerance, tolerance);
        }
        if (maximumTolerance <= 0.0)
            h = 0.0;
        h = Math.max(h, twentySixEpsilon * Math.max(timeRemaining, Math.abs(fieldCurrentTime)));
        // To avoid premature underflow in the error
        // tolerance function,  scale the tolerances...
        final double scale = 2.0 / relativeErrorTolerance;
        final double scaledAbsoluteError = scale * absoluteErrorTolerance;
        boolean previousStepFailed = false;
        // Check for failure...
        check(getValueVector(0));
        updateResultSet();
        // Integrate...
        int iteration = 0;
        while (fieldCurrentTime < endingTime) {
            checkForUserStop();
            // Set smallest allowable stepsize...
            final double minimumTimeStep = Math.max(twentySixEpsilon * Math.abs(fieldCurrentTime), taskDescription.getTimeStep().getMinimumTimeStep());
            h = Math.min(h, taskDescription.getTimeStep().getMaximumTimeStep());
            // Adjust stepsize if necessary to hit the output point.
            // Look ahead two steps to avoid drastic changes in the stepsize and
            // thus lessen the impact of output points on the code.
            timeRemaining = endingTime - fieldCurrentTime;
            if (timeRemaining < 2 * h) {
                if (timeRemaining <= h) {
                    h = timeRemaining;
                } else {
                    h = 0.5 * timeRemaining;
                }
            }
            // If too close to output point, extrapolate and return
            if (timeRemaining <= twentySixEpsilon * Math.abs(fieldCurrentTime)) {
                for (int i = 0; i < getStateVariableCount(); i++) {
                    int I = getVariableIndex(i);
                    newValues[I] = oldValues[I] + timeRemaining * evaluate(oldValues, i);
                }
                // update (old = new)
                fieldValueVectors.copyValuesDown();
                // compute fast system
                if (getFastAlgebraicSystem() != null) {
                    fieldValueVectors.copyValues(0, 1);
                    getFastAlgebraicSystem().initVars(getValueVector(0), getValueVector(1));
                    getFastAlgebraicSystem().solveSystem(getValueVector(0), getValueVector(1));
                    fieldValueVectors.copyValues(1, 0);
                }
                // check for failure
                check(getValueVector(0));
                fieldCurrentTime = endingTime;
                return;
            } else {
                // Advance an approximate solution over one step of length h
                // RungeKuttaStep(Model,y,t,h,ss);
                step(fieldCurrentTime, h);
                // compute fast system
                if (getFastAlgebraicSystem() != null) {
                    fieldValueVectors.copyValues(1, 2);
                    getFastAlgebraicSystem().initVars(getValueVector(1), getValueVector(2));
                    getFastAlgebraicSystem().solveSystem(getValueVector(1), getValueVector(2));
                    fieldValueVectors.copyValues(2, 1);
                }
                // check for failure
                check(getValueVector(1));
                // Compute and test allowable tolerances versus local error estimates
                // and remove scaling of tolerances. Note that relative error is
                // measured with respect to the average of the magnitudes of the
                // solution at the beginning and end of the step.
                double estimatedErrorOverErrorTolerance = 0.0;
                for (int i = 0; i < getStateVariableCount(); i++) {
                    int I = getVariableIndex(i);
                    double errorTolerance = Math.abs(oldValues[I]) + Math.abs(newValues[I]) + scaledAbsoluteError;
                    // Inappropriate error tolerance
                    if (errorTolerance <= 0.0) {
                        throw new SolverException("Error tolerance too small");
                    }
                    double estimatedError = Math.abs(calculateErrorTerm(i));
                    estimatedErrorOverErrorTolerance = Math.max(estimatedErrorOverErrorTolerance, estimatedError / errorTolerance);
                }
                double estimatedTolerance = h * estimatedErrorOverErrorTolerance * scale;
                if (estimatedTolerance > 1.0) {
                    // Unsuccessful step.  Reduce the stepsize and try again.
                    // The decrease is limited to a factor of 1/10...
                    previousStepFailed = true;
                    double s = 0.1;
                    if (estimatedTolerance < 59049.0) {
                        // s = 0.1  @  estimatedTolerance = 59049
                        s = 0.9 / Math.pow(estimatedTolerance, 0.2);
                    }
                    h *= s;
                    if (h < minimumTimeStep) {
                        throw new SolverException("Requested error unattainable at smallest allowable stepsize");
                    }
                } else {
                    // Successful step.  Store solution at t+h and evaluate derivatives there.
                    fieldValueVectors.copyValuesDown();
                    fieldCurrentTime += h;
                    iteration++;
                    if (taskDescription.getOutputTimeSpec().isDefault()) {
                        int keepEvery = ((DefaultOutputTimeSpec) taskDescription.getOutputTimeSpec()).getKeepEvery();
                        if ((iteration % keepEvery) == 0)
                            updateResultSet();
                    }
                    // 
                    // Choose next stepsize.  The increase is limited to a factor of 5. If
                    // step failure has just occurred, next stepsize is not allowed to increase.
                    // 
                    double s = 5;
                    if (estimatedTolerance > 0.0001889568) {
                        // s = 5  @  estimatedTolerance = 0.0001889568
                        s = 0.9 / Math.pow(estimatedTolerance, 0.2);
                    }
                    if (previousStepFailed) {
                        s = Math.min(1.0, s);
                    }
                    h = Math.min(Math.max(minimumTimeStep, s * h), maximumTimeStep);
                    previousStepFailed = false;
                }
            }
        }
        // store last time point
        if (taskDescription.getOutputTimeSpec().isDefault()) {
            int keepEvery = ((DefaultOutputTimeSpec) taskDescription.getOutputTimeSpec()).getKeepEvery();
            if ((iteration % keepEvery) != 0)
                updateResultSet();
        }
    } 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 24 with DefaultOutputTimeSpec

use of cbit.vcell.solver.DefaultOutputTimeSpec in project vcell by virtualcell.

the class NFSimSolver method getMathExecutableCommand.

@Override
protected String[] getMathExecutableCommand() {
    String executableName = null;
    try {
        executableName = SolverUtilities.getExes(SolverDescription.NFSim)[0].getAbsolutePath();
    } catch (IOException e) {
        throw new RuntimeException("failed to get executable for solver " + SolverDescription.NFSim.getDisplayLabel() + ": " + e.getMessage(), e);
    }
    String inputFilename = getInputFilename();
    String outputFilename = getOutputFilename();
    String speciesOutputFilename = getSpeciesOutputFilename();
    NFsimSimulationOptions nfsso = simTask.getSimulation().getSolverTaskDescription().getNFSimSimulationOptions();
    ArrayList<String> adv = new ArrayList<String>();
    boolean observableComputationOff = nfsso.getObservableComputationOff();
    if (observableComputationOff == true) {
        // false is by default, no need to specify
        adv.add("-notf");
    }
    Integer moleculeDistance = nfsso.getMoleculeDistance();
    if (moleculeDistance != null) {
        adv.add("-utl");
        adv.add(moleculeDistance + "");
    }
    boolean aggregateBookkeeping = nfsso.getAggregateBookkeeping();
    if (aggregateBookkeeping == true || simTask.getSimulation().getMathDescription().hasSpeciesObservable()) {
        // false is by default, no need to specify
        adv.add("-cb");
    }
    Integer maxMoleculesPerType = nfsso.getMaxMoleculesPerType();
    if (maxMoleculesPerType != null) {
        adv.add("-gml");
        adv.add(maxMoleculesPerType + "");
    }
    Integer equilibrateTime = nfsso.getEquilibrateTime();
    if (equilibrateTime != null) {
        adv.add("-eq");
        adv.add(equilibrateTime + "");
    }
    boolean preventIntraBonds = nfsso.getPreventIntraBonds();
    if (preventIntraBonds == true) {
        // false is by default, no need to specify
        adv.add("-bscb");
    }
    TimeBounds tb = getSimulationJob().getSimulation().getSolverTaskDescription().getTimeBounds();
    double dtime = tb.getEndingTime() - tb.getStartingTime();
    String timeSpecOption1 = "-oSteps";
    String timeSpecOption2 = "10";
    OutputTimeSpec outputTimeSpec = getSimulationJob().getSimulation().getSolverTaskDescription().getOutputTimeSpec();
    if (outputTimeSpec instanceof DefaultOutputTimeSpec) {
        DefaultOutputTimeSpec dots = (DefaultOutputTimeSpec) outputTimeSpec;
        int steps = dots.getKeepAtMost();
        timeSpecOption1 = "-oSteps";
        timeSpecOption2 = Integer.toString(steps);
    } else if (outputTimeSpec instanceof UniformOutputTimeSpec) {
        UniformOutputTimeSpec dots = (UniformOutputTimeSpec) outputTimeSpec;
        double steps = dtime / dots.getOutputTimeStep();
        timeSpecOption1 = "-oSteps";
        int stepsi = (int) Math.round(steps);
        timeSpecOption2 = Integer.toString(stepsi);
    } else {
        throw new RuntimeException("Unsupported output time spec class");
    }
    String[] baseCommands = { "-xml", inputFilename, "-o", outputFilename, "-sim", Double.toString(dtime), "-ss", speciesOutputFilename };
    ArrayList<String> cmds = new ArrayList<String>();
    cmds.add(executableName);
    Integer seed = nfsso.getRandomSeed();
    if (seed != null) {
        cmds.add("-seed");
        cmds.add(seed.toString());
    } else {
        long randomSeed = System.currentTimeMillis();
        randomSeed = randomSeed + simTask.getSimulationJob().getJobIndex();
        // multiply with a large prime number to spread numbers that are too close and in sequence
        randomSeed = randomSeed * 89611;
        Integer rs = (int) randomSeed;
        String str = rs.toString();
        if (str.startsWith("-")) {
            // NFSim wants a positive integer, for anything else is initializing with 0
            str = str.substring(1);
        }
        cmds.add("-seed");
        cmds.add(str);
    // PrintWriter writer;
    // try {
    // writer = new PrintWriter("c:\\TEMP\\aaa\\" + randomSeed + ".txt", "UTF-8");
    // writer.println(str);
    // writer.close();
    // } catch (FileNotFoundException | UnsupportedEncodingException e) {
    // // Auto-generated catch block
    // e.printStackTrace();
    // }
    }
    cmds.add("-vcell");
    cmds.addAll(new ArrayList<String>(Arrays.asList(baseCommands)));
    cmds.add(timeSpecOption1);
    cmds.add(timeSpecOption2);
    cmds.addAll(adv);
    if (bMessaging) {
        cmds.add("-v");
    }
    return cmds.toArray(new String[cmds.size()]);
}