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

    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();
//.........这里部分代码省略.........