C++ OdeSolution::rGetTimes方法代码示例

  // Test a given solver on an ODE which has a stopping event defined
  void MyTestSolverOnOdesWithEvents(AbstractIvpOdeSolver& rSolver)
  {
        // ODE which has solution y0 = cos(t), and stopping event y0<0,
        // ie should stop when t = pi/2;
        OdeSecondOrderWithEvents ode_with_events;

        OdeSolution solutions;
        std::vector<double> state_variables =
            ode_with_events.GetInitialConditions();
        solutions = rSolver.Solve(&ode_with_events, state_variables, 0.0, 2.0,
            0.001, 0.001);

        unsigned num_timesteps = solutions.GetNumberOfTimeSteps();

        // Final time should be around pi/2
        TS_ASSERT_DELTA( solutions.rGetTimes()[num_timesteps], M_PI_2, 0.01);

        // Penultimate y0 should be greater than zero
        TS_ASSERT_LESS_THAN( 0, solutions.rGetSolutions()[num_timesteps-1][0]);

        // Final y0 should be less than zero
        TS_ASSERT_LESS_THAN( solutions.rGetSolutions()[num_timesteps][0], 0);

        // Solver should correctly state the stopping event occurred
        TS_ASSERT_EQUALS(rSolver.StoppingEventOccurred(), true);

        // This is to cover the exception when a stopping event occurs before the first timestep.
        TS_ASSERT_THROWS_ANYTHING(rSolver.Solve(&ode_with_events, state_variables, 2.0, 3.0, 0.001));

        ///////////////////////////////////////////////
        // Repeat with sampling time larger than dt
        ///////////////////////////////////////////////

        state_variables = ode_with_events.GetInitialConditions();
        solutions = rSolver.Solve(&ode_with_events, state_variables, 0.0, 2.0,
            0.001, 0.01);

        num_timesteps = solutions.GetNumberOfTimeSteps();

        // Final time should be around pi/2
        TS_ASSERT_DELTA( solutions.rGetTimes()[num_timesteps], M_PI_2, 0.01);

        // Penultimate y0 should be greater than zero
        TS_ASSERT_LESS_THAN( 0, solutions.rGetSolutions()[num_timesteps-1][0]);

        // Final y0 should be less than zero
        TS_ASSERT_LESS_THAN( solutions.rGetSolutions()[num_timesteps][0], 0);

        // Solver should correctly state the stopping event occurred
        TS_ASSERT_EQUALS(rSolver.StoppingEventOccurred(), true);

        // Cover the check event isn't initially true exception
        std::vector<double> bad_init_cond;
        bad_init_cond.push_back(-1); //y0 < 0 so stopping event true
        bad_init_cond.push_back(0.0);
        TS_ASSERT_THROWS_ANYTHING(rSolver.Solve(&ode_with_events, bad_init_cond, 0.0, 2.0, 0.001, 0.01));
    }

    void TestSolvingOdes() throw(Exception)
    {
        /* First, create an instance of the ODE class to be solved. */
        MyOde my_ode;
        /* Next, create a solver. */
        EulerIvpOdeSolver euler_solver;
        /* We will need to provide an initial condition, which needs to
         * be a {{{std::vector}}}.*/
        std::vector<double> initial_condition;
        initial_condition.push_back(1.0);
        /* Then, just call `Solve`, passing in a pointer to the ODE, the
         * initial condition, the start time, end time, the solving timestep,
         * and sampling timestep (how often we want the solution stored in the returned `OdeSolution` object).
         * Here we solve from 0 to 1, with a timestep of 0.01 but a ''sampling
         * timestep'' of 0.1. The return value is an object of type {{{OdeSolution}}}
         * (which is basically just a list of times and solutions).
         */
        OdeSolution solutions = euler_solver.Solve(&my_ode, initial_condition, 0, 1, 0.01, 0.1);
        /* Let's look at the results, which can be obtained from the {{{OdeSolution}}}
         * object using the methods {{{rGetTimes()}}} and {{{rGetSolutions()}}}, which
         * return a {{{std::vector}}} and a {{{std::vector}}} of {{{std::vector}}}s
         * respectively. */
        for (unsigned i=0; i<solutions.rGetTimes().size(); i++)
        {
            /* The {{{[0]}}} here is because we are getting the zeroth component of y (a 1-dimensional vector). */
            std::cout << solutions.rGetTimes()[i] << " " << solutions.rGetSolutions()[i][0] << "\n";
        }

        /* Alternatively, we can print the solution directly to a file, using the {{{WriteToFile}}}
         * method on the {{{OdeSolution}}} class. */
        solutions.WriteToFile("SolvingOdesTutorial", "my_ode_solution", "sec");
        /* Two files are written
         * * {{{my_ode_solution.dat}}} contains the results (a header line, then one column for time and one column per variable)
         * * {{{my_ode_solution.info}}} contains information for reading the data back, a line about the ODE solver ("{{{ODE SOLVER: EulerIvpOdeSolver}}}") and provenance information.
         */

        /* We can see from the printed out results that y goes above 2.5 somewhere just
         * before 0.6. To solve only up until y=2.5, we can solve the ODE that has the
         * stopping event defined, using the same solver as before. */
        MyOdeWithStoppingEvent my_ode_stopping;

        /* '''Note:''' ''when a {{{std::vector}}} is passed in as an initial condition
         * to a {{{Solve}}} call, it gets updated as the solve takes place''. Therefore, if
         * we want to use the same initial condition again, we have to reset it back to 1.0. */
        initial_condition[0] = 1.0;
        solutions = euler_solver.Solve(&my_ode_stopping, initial_condition, 0, 1, 0.01, 0.1);
        /* We can check with the solver that it stopped because of the stopping event, rather than because
         * it reached to end time. */
        TS_ASSERT(euler_solver.StoppingEventOccurred());
        /* Finally, let's print the time of the stopping event (to the nearest dt or so). */
        std::cout << "Stopping event occurred at t="<<solutions.rGetTimes().back()<<"\n";
    }

OdeSolution AbstractRushLarsenCardiacCell::Compute(double tStart, double tEnd, double tSamp)
{
    // In this method, we iterate over timesteps, doing the following for each:
    //   - update V using a forward Euler step
    //   - do as in ComputeExceptVoltage(t) to update the remaining state variables
    //     using Rush Larsen method or forward Euler as appropriate

    // Check length of time interval
    if (tSamp < mDt)
    {
        tSamp = mDt;
    }
    const unsigned n_steps = (unsigned) floor((tEnd - tStart)/tSamp + 0.5);
    assert(fabs(tStart+n_steps*tSamp - tEnd) < 1e-12);
    const unsigned n_small_steps = (unsigned) floor(tSamp/mDt+0.5);
    assert(fabs(mDt*n_small_steps - tSamp) < 1e-12);

    // Initialise solution store
    OdeSolution solutions;
    solutions.SetNumberOfTimeSteps(n_steps);
    solutions.rGetSolutions().push_back(rGetStateVariables());
    solutions.rGetTimes().push_back(tStart);
    solutions.SetOdeSystemInformation(this->mpSystemInfo);

    std::vector<double> dy(mNumberOfStateVariables, 0);
    std::vector<double> alpha(mNumberOfStateVariables, 0);
    std::vector<double> beta(mNumberOfStateVariables, 0);

    // Loop over time
    for (unsigned i=0; i<n_steps; i++)
    {
        double curr_time = tStart;
        for (unsigned j=0; j<n_small_steps; j++)
        {
            curr_time = tStart + i*tSamp + j*mDt;
            EvaluateEquations(curr_time, dy, alpha, beta);
            UpdateTransmembranePotential(dy);
            ComputeOneStepExceptVoltage(dy, alpha, beta);
            VerifyStateVariables();
        }

        // Update solutions
        solutions.rGetSolutions().push_back(rGetStateVariables());
        solutions.rGetTimes().push_back(curr_time+mDt);
    }

    return solutions;
}

    void TestGarysWntOdeSystemApc2Hit()
    {
#ifdef CHASTE_CVODE
        double wnt_level = 0.5;
        boost::shared_ptr<AbstractCellMutationState> p_apc2(new ApcTwoHitCellMutationState);
        Mirams2010WntOdeSystem wnt_system(wnt_level, p_apc2);

        // Solve system using CVODE solver
        // Matlab's strictest bit uses 0.01 below and relaxes it on flatter bits
        double h_value = 0.01;
        CvodeAdaptor cvode_solver;
        OdeSolution solutions;
        std::vector<double> initial_conditions = wnt_system.GetInitialConditions();

        Timer::Reset();
        solutions = cvode_solver.Solve(&wnt_system, initial_conditions, 0.0, 100.0, h_value, h_value);
        Timer::Print("1. Cvode");

        // Test solutions are OK for a small time increase
        int end = solutions.rGetSolutions().size() - 1;

        // Test the simulation is ending at the right time (going into S phase at 7.8 hours)
        TS_ASSERT_DELTA(solutions.rGetTimes()[end], 100, 1e-2);

        // Check results are correct
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][0], 433.114, 2e-3); // Tolerances relaxed for
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][1], 433.114, 2e-3); // different CVODE versions.
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][2], wnt_level, 1e-4);
#else
        std::cout << "CVODE is not enabled. " << std::endl;
        std::cout << "If required please install and alter your hostconfig settings to switch on chaste support." << std::endl;
#endif //CHASTE_CVODE
    }

    void TestBackwardEulerSystemOf3EquationsWithEvents()
    {
        OdeThirdOrderWithEvents ode_system_with_events;

        double h_value = 0.01;

        // Euler solver solution worked out
        BackwardEulerIvpOdeSolver backward_euler_solver(ode_system_with_events.GetNumberOfStateVariables());
        OdeSolution solutions;

        std::vector<double> state_variables = ode_system_with_events.GetInitialConditions();
        solutions = backward_euler_solver.Solve(&ode_system_with_events, state_variables, 0.0, 2.0, h_value, h_value);
        unsigned last = solutions.GetNumberOfTimeSteps();

        // Final time should be pi/6 (?)
        TS_ASSERT_DELTA( solutions.rGetTimes()[last], 0.5236, 0.01);

        // Penultimate y0 should be greater than -0.5
        TS_ASSERT_LESS_THAN(-0.5,solutions.rGetSolutions()[last-1][0]);

        // Final y0 should be less than -0.5
        TS_ASSERT_LESS_THAN( solutions.rGetSolutions()[last][0], -0.5);

        // Solver should correctly state the stopping event occurred
        TS_ASSERT_EQUALS(backward_euler_solver.StoppingEventOccurred(), true);
    }

    /*
     * === Solving n-dimensional ODEs ===
     *
     * Finally, here's a simple test showing how to solve a 2d ODE using the first method.
     * All that is different is the initial condition has be a vector of length 2, and returned
     * solution is of length 2 at every timestep.
     */
    void TestWith2dOde()
    {
        My2dOde my_2d_ode;
        EulerIvpOdeSolver euler_solver;

        /* Define the initial condition for each state variable. */
        std::vector<double> initial_condition;
        initial_condition.push_back(1.0);
        initial_condition.push_back(0.0);

        /* Solve, and print the solution as [time, y1, y2]. */
        OdeSolution solutions = euler_solver.Solve(&my_2d_ode, initial_condition, 0, 1, 0.01, 0.1);
        for (unsigned i=0; i<solutions.rGetTimes().size(); i++)
        {
            std::cout << solutions.rGetTimes()[i] << " "
                      << solutions.rGetSolutions()[i][0] << " "
                      << solutions.rGetSolutions()[i][1] << "\n";
        }
    }

    void Simulate(const std::string& rOutputDirName,
                  const std::string& rModelName,
                  boost::shared_ptr<AbstractCardiacCellInterface> pCell)
    {
        double end_time = GetAttribute(pCell, "SuggestedCycleLength", 700.0); // ms
        if (pCell->GetSolver() || dynamic_cast<AbstractRushLarsenCardiacCell*>(pCell.get()))
        {
            double dt = GetAttribute(pCell, "SuggestedForwardEulerTimestep", 0.0);
            if (dt > 0.0)
            {
                pCell->SetTimestep(dt);
            }
        }
#ifdef CHASTE_CVODE
        AbstractCvodeSystem* p_cvode_cell = dynamic_cast<AbstractCvodeSystem*>(pCell.get());
        if (p_cvode_cell)
        {
            // Set a larger max internal time steps per sampling interval (CVODE's default is 500)
            p_cvode_cell->SetMaxSteps(1000);
            // Numerical or analytic J for CVODE?
            if (!mUseCvodeJacobian)
            {
                p_cvode_cell->ForceUseOfNumericalJacobian();
            }
        }
#endif
        double sampling_interval = 1.0; // ms; used as max dt for CVODE too
        Timer::Reset();
        OdeSolution solution = pCell->Compute(0.0, end_time, sampling_interval);
        std::stringstream message;
        message << "Model " << rModelName << " writing to " << rOutputDirName << " took";
        Timer::Print(message.str());

        const unsigned output_freq = 10; // Only output every N samples
        solution.WriteToFile(rOutputDirName, rModelName, "ms", output_freq, false);
        // Check an AP was produced
        std::vector<double> voltages = solution.GetVariableAtIndex(pCell->GetVoltageIndex());
        CellProperties props(voltages, solution.rGetTimes());
        props.GetLastActionPotentialDuration(90.0); // Don't catch the exception here if it's thrown
        // Compare against saved results
        CheckResults(rModelName, voltages, solution.rGetTimes(), output_freq);
    }

    /**
     * Test two ODE solvers with this ODE system (correct values calculated using the Matlab solver ode15s).
     *
     */
    void TestAlarcon2004Solver()
    {
        // Set up
        double oxygen_concentration = 1.0;
        Alarcon2004OxygenBasedCellCycleOdeSystem alarcon_system(oxygen_concentration, false);

        // Create ODE solvers
        RungeKutta4IvpOdeSolver rk4_solver;
        RungeKuttaFehlbergIvpOdeSolver rkf_solver;
        BackwardEulerIvpOdeSolver back_solver(6);

        // Set up for solver
        OdeSolution solutions;
        std::vector<double> initial_conditions = alarcon_system.GetInitialConditions();
        double h_value = 1e-4; // maximum tolerance

        // Solve the ODE system using a Runge Kutta fourth order solver
        Timer::Reset();
        solutions = rk4_solver.Solve(&alarcon_system, initial_conditions, 0.0, 10.0, h_value, h_value);
        Timer::Print("1. Runge-Kutta");

        // Reset maximum tolerance for Runge Kutta Fehlber solver
        h_value = 1e-1;

        // Solve the ODE system using a Runge Kutta Fehlber solver
        initial_conditions = alarcon_system.GetInitialConditions();
        Timer::Reset();
        solutions = rkf_solver.Solve(&alarcon_system, initial_conditions, 0.0, 10.0, h_value, 1e-4);
        Timer::Print("1. Runge-Kutta-Fehlberg");

        // Test that solutions are accurate for a small time increase
        int end = solutions.rGetSolutions().size() - 1;

        // Test that the solver stops at the right time
        TS_ASSERT_DELTA(solutions.rGetTimes()[end], 9.286356375, 1e-2);

        // Test solution - note the high tolerances
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][0], 0.004000000000000, 1e-3);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][1], 0.379221366479055, 1e-3);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][2], 0.190488726735972, 1e-3);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][3], 9.962110289977730, 1e-3);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][4], 0.096476600742599, 1e-3);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][5], 1.000000000000000, 1e-3);
    }

    void TestRKFehlbergSystemOf3EquationsWithEvents() throw(Exception)
    {
        OdeThirdOrderWithEvents ode_system_with_events;

        double h_value = 0.1;

        // Euler solver solution worked out
        RungeKuttaFehlbergIvpOdeSolver rkf_solver;
        OdeSolution solutions;

        std::vector<double> state_variables = ode_system_with_events.GetInitialConditions();
        solutions = rkf_solver.Solve(&ode_system_with_events, state_variables, 0.0, 2.0, h_value, 1e-5);
        unsigned last = solutions.GetNumberOfTimeSteps();

//        for (unsigned i=0; i<last+1; i++)
//        {
//            std::cout << "Time = " << solutions.rGetTimes()[i] <<
//                " x = " << solutions.rGetSolutions()[i][0] << "\n" << std::flush;
//        }

        // Final time should be pi/6 (?)
        TS_ASSERT_DELTA( solutions.rGetTimes()[last], M_PI/6.0, h_value);

        // Penultimate y0 should be greater than -0.5
        TS_ASSERT_LESS_THAN(-0.5,solutions.rGetSolutions()[last-1][0]);

        // Final y0 should be less than -0.5
        TS_ASSERT_LESS_THAN( solutions.rGetSolutions()[last][0], -0.5);

        // Solver should correctly state the stopping event occurred
        TS_ASSERT_EQUALS(rkf_solver.StoppingEventOccurred(), true);

        // Coverage of exceptions
        TS_ASSERT_THROWS_THIS(rkf_solver.Solve(&ode_system_with_events, solutions.rGetSolutions()[last], M_PI/6.0, 2.0, 0.1, 1e-5),
                "(Solve with sampling) Stopping event is true for initial condition");
        TS_ASSERT_THROWS_THIS(rkf_solver.Solve(&ode_system_with_events, solutions.rGetSolutions()[last], M_PI/6.0, 2.0, 0.1),
                "(Solve without sampling) Stopping event is true for initial condition");
    }

    void TestGarysWntOdeSystemBetaCatenin1Hit() throw(Exception)
    {
#ifdef CHASTE_CVODE
        double wnt_level = 0.5;
        boost::shared_ptr<AbstractCellMutationState> p_bcat1(new BetaCateninOneHitCellMutationState);
        Mirams2010WntOdeSystem wnt_system(wnt_level, p_bcat1);

        // Solve system using CVODE solver
        // Matlab's strictest bit uses 0.01 below and relaxes it on flatter bits
        double h_value = 0.1;
        CvodeAdaptor cvode_solver;
        OdeSolution solutions;
        std::vector<double> initial_conditions = wnt_system.GetInitialConditions();

        double start_time, end_time, elapsed_time = 0.0;
        start_time = (double) std::clock();
        solutions = cvode_solver.Solve(&wnt_system, initial_conditions, 0.0, 100.0, h_value, h_value);
        end_time = (double) std::clock();
        elapsed_time = (end_time - start_time)/(CLOCKS_PER_SEC);
        std::cout << "1. Cvode Elapsed time = " << elapsed_time << " secs for 100 hours\n";

        // Test solutions are OK for a small time increase
        int end = solutions.rGetSolutions().size() - 1;

        // Tests the simulation is ending at the right time (going into S phase at 7.8 hours)
        TS_ASSERT_DELTA(solutions.rGetTimes()[end], 100, 1e-2);

        // Check results are correct
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][0], 67.5011, 1e-4);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][1], 824.0259, 1e-4);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][2], wnt_level, 1e-4);
#else
        std::cout << "CVODE is not enabled. " << std::endl;
        std::cout << "If required please install and alter your hostconfig settings to switch on chaste support." << std::endl;
#endif //CHASTE_CVODE
    }

    void TestRKFehlbergWithExampleFromBook() throw(Exception)
    {
        /*
         * Book is "Numerical Analysis 6th Edition by R.L. Burden and J. D. Faires
         * This example is on P291 Table 5.9
         */
        RkfTestOde ode;

        double max_step_size = 0.25;
        double start_time = 0.0;
        double end_time = 2.0;
        RungeKuttaFehlbergIvpOdeSolver rkf_solver;

        OdeSolution solutions;

        std::vector<double> state_variables = ode.GetInitialConditions();
        double tolerance = 1e-5;
        solutions = rkf_solver.Solve(&ode, state_variables, start_time, end_time, max_step_size, tolerance);

        // Times (from MatLab Code) to check timstepping is being adapted properly
        TS_ASSERT_DELTA(solutions.rGetTimes()[0], 0, 1e-7);
        TS_ASSERT_DELTA(solutions.rGetTimes()[1], 2.500000000000000e-01, 1e-7);
        TS_ASSERT_DELTA(solutions.rGetTimes()[2], 4.868046415733731e-01, 1e-7);
        TS_ASSERT_DELTA(solutions.rGetTimes()[3], 7.298511818781566e-01, 1e-7);
        TS_ASSERT_DELTA(solutions.rGetTimes()[4], 9.798511818781566e-01, 1e-7);
        TS_ASSERT_DELTA(solutions.rGetTimes()[5], 1.229851181878157e+00, 1e-7);
        TS_ASSERT_DELTA(solutions.rGetTimes()[6], 1.479851181878157e+00, 1e-7);
        TS_ASSERT_DELTA(solutions.rGetTimes()[7], 1.729851181878157e+00, 1e-7);
        TS_ASSERT_DELTA(solutions.rGetTimes()[8], 1.979851181878157e+00, 1e-7);
        TS_ASSERT_DELTA(solutions.rGetTimes()[9], 2.000000000000000e+00, 1e-7);

        TS_ASSERT_EQUALS(solutions.GetNumberOfTimeSteps(), 9u);

        // y values (from analytic result)
        for (unsigned i=0; i<solutions.GetNumberOfTimeSteps(); i++)
        {
            double time = solutions.rGetTimes()[i];
            double y = (time+1.0)*(time+1.0) - 0.5*exp(time);

            // Tolerance set to 1e-5, so 2e-5 to pass here
            TS_ASSERT_DELTA(solutions.rGetSolutions()[i][0], y, 2e-5);
        }
    }

    void TestTysonNovakSolver() throw(Exception)
    {
        TysonNovak2001OdeSystem tyson_novak_system;

        // Solve system using backward Euler solver

        // Matlab's strictest bit uses 0.01 below and relaxes it on flatter bits

        double dt = 0.1/60.0;

        //Euler solver solution worked out
        BackwardEulerIvpOdeSolver backward_euler_solver(6);

        std::vector<double> state_variables = tyson_novak_system.GetInitialConditions();

        Timer::Reset();
        OdeSolution solutions = backward_euler_solver.Solve(&tyson_novak_system, state_variables, 0.0, 75.8350/60.0, dt, dt);
        Timer::Print("1. Tyson Novak Backward Euler");

        // If you run it up to about 75min the ODE will stop, anything less and it will not and this test will fail
        TS_ASSERT_EQUALS(backward_euler_solver.StoppingEventOccurred(), true);

        unsigned end = solutions.rGetSolutions().size() - 1;

        // The following code provides nice output for gnuplot
        // use the command
        // plot "tyson_novak.dat" u 1:2
        // or
        // plot "tyson_novak.dat" u 1:3 etc. for the various proteins...

//        OutputFileHandler handler("");
//        out_stream file=handler.OpenOutputFile("tyson_novak.dat");
//        for (unsigned i=0; i<=end; i++)
//        {
//            (*file) << solutions.rGetTimes()[i]<< "\t" << solutions.rGetSolutions()[i][0] << "\t" << solutions.rGetSolutions()[i][1] << "\t" << solutions.rGetSolutions()[i][2] << "\t" << solutions.rGetSolutions()[i][3] << "\t" << solutions.rGetSolutions()[i][4] << "\t" << solutions.rGetSolutions()[i][5] << "\n" << std::flush;
//        }
//        file->close();

        ColumnDataWriter writer("TysonNovak", "TysonNovak");
        if (PetscTools::AmMaster()) // if master process
        {
            int step_per_row = 1;
            int time_var_id = writer.DefineUnlimitedDimension("Time", "s");

            std::vector<int> var_ids;
            for (unsigned i=0; i<tyson_novak_system.rGetStateVariableNames().size(); i++)
            {
                var_ids.push_back(writer.DefineVariable(tyson_novak_system.rGetStateVariableNames()[i],
                                                        tyson_novak_system.rGetStateVariableUnits()[i]));
            }
            writer.EndDefineMode();

            for (unsigned i = 0; i < solutions.rGetSolutions().size(); i+=step_per_row)
            {
                writer.PutVariable(time_var_id, solutions.rGetTimes()[i]);
                for (unsigned j=0; j<var_ids.size(); j++)
                {
                    writer.PutVariable(var_ids[j], solutions.rGetSolutions()[i][j]);
                }
                writer.AdvanceAlongUnlimitedDimension();
            }
            writer.Close();
        }
        PetscTools::Barrier();

        // Proper values calculated using the Matlab stiff ODE solver ode15s. Note that
        // large tolerances are required for the tests to pass with both chaste solvers
        // and CVODE.
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][0],0.10000000000000, 1e-2);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][1],0.98913684535843, 1e-2);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][2],1.54216806705641, 1e-1);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][3],1.40562614481544, 1e-1);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][4],0.67083371879876, 1e-2);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][5],0.95328206604519, 2e-2);
    }

    void TestInterpolatorTimesAndGenerateReferenceTrace() throw(Exception)
    {
#ifdef CHASTE_CVODE
        OutputFileHandler handler("CvodeCellsWithDataClamp");

        boost::shared_ptr<AbstractIvpOdeSolver> p_empty_solver;
        boost::shared_ptr<AbstractStimulusFunction> p_empty_stimulus;

        // N.B. Because we use the Shannon model as a lot of examples,
        // here it is actually a Shannon->WithModifiers->WithDataClamp->CvodeCell
        // (the WithModifiers doesn't need to be there to use the data clamp!)
        mpModel.reset(new CellShannon2004FromCellMLCvodeDataClamp(p_empty_solver,p_empty_stimulus));

        TS_ASSERT_EQUALS(mpModel->HasParameter("membrane_data_clamp_current_conductance"), true);

        mpModel->SetMaxSteps(5000);
        mpModel->UseCellMLDefaultStimulus();

        // Run a simulation without clamping switched on

        Timer::Reset();
        double end_time = 400.0;
        OdeSolution solution = mpModel->Compute(0, end_time, 0.2);
        Timer::Print("OdeSolution");
        std::vector<double> expt_times = solution.rGetTimes();
        std::vector<double> expt_data = solution.GetAnyVariable("membrane_voltage");
        solution.WriteToFile("CvodeCellsWithDataClamp","shannon_original_no_clamp", "ms", 1, false); // false to clean

        TS_ASSERT_THROWS_THIS(mpModel->TurnOnDataClamp(),
            "Before calling TurnOnDataClamp(), please provide experimental data via the SetExperimentalData() method.");

        // Test the interpolation methods.
        {
            mpModel->SetExperimentalData(expt_times, expt_data);

            // Note - unless the data clamp is switched on the below method just returns DOUBLE_UNSET to save time interpolating.
            double time = 100.0;
            TS_ASSERT_EQUALS(mpModel->GetExperimentalVoltageAtTimeT(time), DOUBLE_UNSET);

            // So now turn on the data clamp
            mpModel->TurnOnDataClamp();

# if CHASTE_SUNDIALS_VERSION >= 20400
            double tol = 5e-3; // mV
#else
            double tol = 0.2; // mV
#endif
            TS_ASSERT_DELTA(mpModel->GetExperimentalVoltageAtTimeT(time), -8.55863245e+01, tol);

            // So turn it off again
            mpModel->TurnOffDataClamp();
            TS_ASSERT_DELTA(mpModel->GetParameter("membrane_data_clamp_current_conductance"), 0.0, 1e-12);
            mpModel->TurnOnDataClamp(200.0);
            TS_ASSERT_DELTA(mpModel->GetParameter("membrane_data_clamp_current_conductance"), 200.0, 1e-12);
            mpModel->TurnOffDataClamp();
            TS_ASSERT_DELTA(mpModel->GetParameter("membrane_data_clamp_current_conductance"), 0.0, 1e-12);
            mpModel->TurnOnDataClamp();
            TS_ASSERT_DELTA(mpModel->GetParameter("membrane_data_clamp_current_conductance"), 100.0, 1e-12); // the default

            // Test a couple of times where no interpolation is needed (on data points).
            time = 116.0;
            double v_at_116 = 1.53670634e+01;
            TS_ASSERT_DELTA(mpModel->GetExperimentalVoltageAtTimeT(time), v_at_116, tol);

            time = 116.2;
            double v_at_116_2 = 1.50089546e+01;
            TS_ASSERT_DELTA(mpModel->GetExperimentalVoltageAtTimeT(time), v_at_116_2, tol);

            // Now test a time where interpolation is required.
            time = 116.1;
            TS_ASSERT_DELTA(mpModel->GetExperimentalVoltageAtTimeT(time), 0.5*(v_at_116 + v_at_116_2), tol);

            // Test ends
            TS_ASSERT_DELTA(mpModel->GetExperimentalVoltageAtTimeT(0.0), expt_data[0], 1e-4);
            TS_ASSERT_DELTA(mpModel->GetExperimentalVoltageAtTimeT(end_time), expt_data.back(), 1e-4);

            // Test exceptions
            TS_ASSERT_THROWS_CONTAINS(mpModel->GetExperimentalVoltageAtTimeT(-1e-12),
                                      "is outside the times stored in the data clamp");
            TS_ASSERT_THROWS_CONTAINS(mpModel->GetExperimentalVoltageAtTimeT(end_time+1e-12),
                                      "is outside the times stored in the data clamp");

            //std::cout << "membrane_data_clamp_current_conductance = " << mpModel->GetParameter("membrane_data_clamp_current_conductance") << std::endl << std::flush;
            //std::cout << "mpModel->GetExperimentalVoltageAtTimeT(time) = " << mpModel->GetExperimentalVoltageAtTimeT(time) << std::endl << std::flush;

            unsigned how_many = 10000u;
            Timer::Reset();
            for (unsigned i=0; i<how_many; i++)
            {
                mpModel->GetExperimentalVoltageAtTimeT(time);
            }
            Timer::PrintAndReset("GetExperimentalVoltageAtTimeT");
        }

        // Generate some noisier data - more like a real cell.
        out_stream experimental_voltage_results_file = handler.OpenOutputFile("Shannon_noisy_data.dat");
        double experimental_noise_sd = 0.25; // directly from Teun's AP data trace of a stationary-looking bit
        for (unsigned i=0; i<expt_data.size(); i++)
        {
            double random_number = RandomNumberGenerator::Instance()->StandardNormalRandomDeviate();
//.........这里部分代码省略.........

void RungeKuttaFehlbergIvpOdeSolver::InternalSolve(OdeSolution& rSolution,
                                                AbstractOdeSystem* pOdeSystem,
                                                std::vector<double>& rYValues,
                                                std::vector<double>& rWorkingMemory,
                                                double startTime,
                                                double endTime,
                                                double maxTimeStep,
                                                double minTimeStep,
                                                double tolerance,
                                                bool outputSolution)
{
    const unsigned number_of_variables = pOdeSystem->GetNumberOfStateVariables();
    mError.resize(number_of_variables);
    mk1.resize(number_of_variables);
    mk2.resize(number_of_variables);
    mk3.resize(number_of_variables);
    mk4.resize(number_of_variables);
    mk5.resize(number_of_variables);
    mk6.resize(number_of_variables);
    myk2.resize(number_of_variables);
    myk3.resize(number_of_variables);
    myk4.resize(number_of_variables);
    myk5.resize(number_of_variables);
    myk6.resize(number_of_variables);

    double current_time = startTime;
    double time_step = maxTimeStep;
    bool got_to_end = false;
    bool accurate_enough = false;
    unsigned number_of_time_steps = 0;

    if (outputSolution)
    {   // Write out ICs
        rSolution.rGetTimes().push_back(current_time);
        rSolution.rGetSolutions().push_back(rYValues);
    }

    // should never get here if this bool has been set to true;
    assert(!mStoppingEventOccurred);
    while (!got_to_end)
    {
        //std::cout << "New timestep\n" << std::flush;
        while (!accurate_enough)
        {
            accurate_enough = true; // assume it is OK until we check and find otherwise

            // Function that calls the appropriate one-step solver
            CalculateNextYValue(pOdeSystem,
                                time_step,
                                current_time,
                                rYValues,
                                rWorkingMemory);

            // Find the maximum error in this vector
            double max_error = -DBL_MAX;
            for (unsigned i=0; i<number_of_variables; i++)
            {
                if (mError[i] > max_error)
                {
                    max_error = mError[i];
                }
            }

            if (max_error > tolerance)
            {// Reject the step-size and do it again.
                accurate_enough = false;
                //std::cout << "Approximation rejected\n" << std::flush;
            }
            else
            {
                // step forward the time since step has now been made
                current_time = current_time + time_step;
                //std::cout << "Approximation accepted with time step = "<< time_step << "\n" << std::flush;
                //std::cout << "max_error = " << max_error << " tolerance = " << tolerance << "\n" << std::flush;
                if (outputSolution)
                {   // Write out ICs
                    //std::cout << "In solver Time = " << current_time << " y = " << rWorkingMemory[0] << "\n" << std::flush;
                    rSolution.rGetTimes().push_back(current_time);
                    rSolution.rGetSolutions().push_back(rWorkingMemory);
                    number_of_time_steps++;
                }
            }

            // Set a new step size based on the accuracy here
            AdjustStepSize(time_step, max_error, tolerance, maxTimeStep, minTimeStep);
        }

        // For the next timestep check the step doesn't go past the end...
        if (current_time + time_step > endTime)
        {   // Allow a smaller timestep for the final step.
            time_step = endTime - current_time;
        }

        if ( pOdeSystem->CalculateStoppingEvent(current_time,
                                                rWorkingMemory) == true )
        {
            mStoppingTime = current_time;
            mStoppingEventOccurred = true;
        }

//.........这里部分代码省略.........

//.........这里部分代码省略.........
         * CVODE errors when trying to run simulations, it can be worth switching off the analytic Jacobian and resorting
         * to a numerical approximation (as happens by default if no analytic Jacobian is available). This can be done with the
         * following command:
         *
         * {{{p_model->ForceUseOfNumericalJacobian();}}}
         *
         */

        /*
         * == Changing Parameters in the Cell Model ==
         *
         * You can also change any parameters that are labelled in the cell model.
         *
         * Instructions for annotating parameters can be found at [wiki:ChasteGuides/CodeGenerationFromCellML]
         *
         * Here we show how to change the parameter dictating the maximal conductance of the IKs current.
         * Note this call actually leaves it unchanged from the default,
         * you can experiment with changing it and examine the impact on APD.
         */
        p_model->SetParameter("membrane_slow_delayed_rectifier_potassium_current_conductance", 0.07);

        /*
         * == Running model to steady state ==
         *
         * Now we run the model to steady state.
         * You can detect for steady state alternans by giving it true as a second parameter
         * {{{SteadyStateRunner steady_runner(p_model, true);}}}
         *
         * You may change the number of maximum paces the runner takes. The default is 1e5.
         */
        SteadyStateRunner steady_runner(p_model);
        steady_runner.SetMaxNumPaces(100u);
        bool result;
        result = steady_runner.RunToSteadyState();

        /*
         * Check that the model has NOT reached steady state
         * (this model needs more than 100 paces to reach steady state).
         *
         */
        TS_ASSERT_EQUALS(result,false);

        /*
         * == Getting detail for paces of interest ==
         *
         * Now we solve for the number of paces we are interested in.
         *
         * The absolute values of start time and end time are typically only relevant for the stimulus, in general
         * nothing else on the right-hand side of the equations uses time directly.
         *
         * i.e. if you have a `RegularStimulus` of period 1000ms then you would get exactly the same results
         * calling Solve(0,1000,...) twice, as you would calling Solve(0,1000,...) and Solve(1000,2000,...).
         *
         * Single cell results can be very sensitive to the sampling time step, because of the steepness of the upstroke.
         *
         * For example, try changing the line below to 1 ms. The upstroke velocity that is detected will change
         * from 339 mV/ms to around 95 mV/ms. APD calculations will only ever be accurate to sampling timestep
         * for the same reason.
         */
        double max_timestep = 0.1;
        p_model->SetMaxTimestep(max_timestep);

        double sampling_timestep = max_timestep;
        double start_time = 0.0;
        double end_time = 1000.0;
        OdeSolution solution = p_model->Compute(start_time, end_time, sampling_timestep);

        /*
         * `p_model` retains the state variables at the end of `Solve`, if you call `Solve` again the state
         * variables will evolve from their new state, not the original initial conditions.
         *
         * Write the data out to a file.
         */
        solution.WriteToFile("TestCvodeCells","Shannon2004Cvode","ms");

        /*
         * == Calculating APD and Upstroke Velocity ==
         *
         * Calculate APD and upstroke velocity using {{{CellProperties}}}
         */
        unsigned voltage_index = p_model->GetSystemInformation()->GetStateVariableIndex("membrane_voltage");
        std::vector<double> voltages = solution.GetVariableAtIndex(voltage_index);
        CellProperties cell_props(voltages, solution.rGetTimes());

        double apd = cell_props.GetLastActionPotentialDuration(90);
        double upstroke_velocity = cell_props.GetLastMaxUpstrokeVelocity();
        /*
         * Here we just check that the values are equal to the ones we expect,
         * with appropriate precision to pass on different versions of CVODE.
         */
        TS_ASSERT_DELTA(apd, 212.41, 1e-2);
        TS_ASSERT_DELTA(upstroke_velocity, 338, 1.25);

        /* CVODE is still an optional dependency for Chaste.
         * If CVODE is not installed this tutorial will
         * not do anything, but we can at least alert the user to this.*/
#else
        std::cout << "Cvode is not enabled.\n";
#endif
    }