RunAndCheckIonicModels.hpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 
00030 #ifndef _RUNANDCHECKIONICMODELS_HPP_
00031 #define _RUNANDCHECKIONICMODELS_HPP_
00032 
00033 #include <vector>
00034 #include <string>
00035 
00036 #include "OdeSolution.hpp"
00037 
00038 #include "ColumnDataWriter.hpp"
00039 #include "ColumnDataReader.hpp"
00040 
00041 #include "AbstractCardiacCell.hpp"
00042 #include "HeartConfig.hpp"
00043 
00044 void RunOdeSolverWithIonicModel(AbstractCardiacCell* pOdeSystem,
00045                                 double endTime,
00046                                 std::string filename,
00047                                 int stepPerRow=100,
00048                                 bool doComputeExceptVoltage=true,
00049                                 bool useSamplingInterval=false);
00050 
00051 void CheckCellModelResults(const std::string& rBaseResultsFilename,
00052                            std::string validResultsBasename = "",
00053                            double tolerance = 1e-3);
00054 
00055 std::vector<double> GetVoltages(ColumnDataReader& rReader);
00056 
00057 void CompareCellModelResults(std::string baseResultsFilename1, std::string baseResultsFilename2,
00058                              double tolerance, bool vOnly=false, std::string folderName="TestIonicModels");
00059 
00060 
00061 /*
00062  * Note: we have to have the function definitions here, rather than in a .cpp file, if we're
00063  * to include them in heart/src, as otherwise linking of the heart library fails because
00064  * CxxTest routines are not defined.
00065  */
00066 
00067 #include <cmath>
00068 
00069 void RunOdeSolverWithIonicModel(AbstractCardiacCell* pOdeSystem,
00070                                 double endTime,
00071                                 std::string filename,
00072                                 int stepPerRow,
00073                                 bool doComputeExceptVoltage,
00074                                 bool useSamplingInterval)
00075 {
00076     double start_time = 0.0;
00077 
00078     if (doComputeExceptVoltage)
00079     {
00080         // Store the current system state
00081         std::vector<double> state_variables_copy = pOdeSystem->rGetStateVariables();
00082 
00083         // Test ComputeExceptVoltage
00084         double v_init = pOdeSystem->GetVoltage();
00085         pOdeSystem->ComputeExceptVoltage(start_time, start_time + 100.0 /*ms*/);
00086         double v_end = pOdeSystem->GetVoltage();
00087         TS_ASSERT_DELTA(v_init, v_end, 1e-6);
00088 
00089         // Test SetVoltage
00090         pOdeSystem->SetVoltage(1e6);
00091         TS_ASSERT_DELTA(pOdeSystem->GetVoltage(), 1e6, 1e-6);
00092 
00093         // Reset the system
00094         pOdeSystem->SetStateVariables(state_variables_copy);
00095     }
00096 
00097     // Solve and write to file
00098     if (useSamplingInterval)
00099     {
00100         OdeSolution solution = pOdeSystem->Compute(start_time, endTime, HeartConfig::Instance()->GetOdeTimeStep() * stepPerRow);
00101         solution.WriteToFile("TestIonicModels", filename, "ms", 1, false, 4);
00102     }
00103     else
00104     {
00105         OdeSolution solution = pOdeSystem->Compute(start_time, endTime);
00106         solution.WriteToFile("TestIonicModels", filename, "ms", stepPerRow, false, 4);
00107     }
00108 }
00109 
00110 std::vector<double> GetVoltages(ColumnDataReader& rReader)
00111 {
00112     //Rather Ugly, we can't currently guarantee what the name of the voltage column is,
00113     //hence we try to cover the most common possibilities
00114     std::vector<double> voltages;
00115     if (rReader.HasValues("V"))
00116     {
00117         voltages = rReader.GetValues("V");
00118     }
00119     else if (rReader.HasValues("membrane_voltage"))
00120     {
00121         voltages = rReader.GetValues("membrane_voltage");
00122     }
00123     else if (rReader.HasValues("membrane__V"))
00124     {
00125         voltages = rReader.GetValues("membrane__V");
00126     }
00127     else
00128     {
00129         EXCEPTION("Model membrane voltage not recognised.");
00130     }
00131     return voltages;
00132 }
00133 
00134 std::vector<double> GetIntracellularCalcium(ColumnDataReader& rReader)
00135 {
00136     //Rather Ugly, we can't currently guarantee what the name of the calcium column is,
00137     //hence we try to cover the most common possibilities
00138     std::vector<double> cai;
00139     if (rReader.HasValues("CaI"))
00140     {
00141         cai = rReader.GetValues("CaI");
00142     }
00143     else if (rReader.HasValues("cytosolic_calcium_concentration"))
00144     {
00145         cai = rReader.GetValues("cytosolic_calcium_concentration");
00146     }
00147     else if (rReader.HasValues("Cai"))
00148     {
00149         cai = rReader.GetValues("Cai");
00150     }
00151     else
00152     {
00153         EXCEPTION("Model intracellular calcium is not recognised.");
00154     }
00155     return cai;
00156 }
00157 
00158 std::vector<double> GetHGate(ColumnDataReader& rReader)
00159 {
00160     //Rather Ugly, we can't currently guarantee what the name of the h gate column is,
00161     //hence we try to cover the most common possibilities
00162     std::vector<double> h_values;
00163     if (rReader.HasValues("h"))
00164     {
00165         h_values = rReader.GetValues("h");
00166     }
00167     else if (rReader.HasValues("fast_sodium_current_h_gate__h"))
00168     {
00169         h_values = rReader.GetValues("fast_sodium_current_h_gate__h");
00170     }
00171     else if (rReader.HasValues("membrane_fast_sodium_current_h_gate"))
00172     {
00173         h_values = rReader.GetValues("membrane_fast_sodium_current_h_gate");
00174     }
00175     else
00176     {
00177         EXCEPTION("Model h gate is not recognised.");
00178     }
00179     return h_values;
00180 }
00181 
00182 /*
00183  * Check the cell model against a previous version
00184  * or another source e.g. Alan's COR
00185  */
00186 void CheckCellModelResults(const std::string& rBaseResultsFilename,
00187                            std::string validResultsBasename,
00188                            double tolerance)
00189 {
00190     // read data entries for the new file
00191     ColumnDataReader data_reader("TestIonicModels", rBaseResultsFilename);
00192     std::vector<double> times = data_reader.GetValues("Time");
00193     std::vector<double> voltages = GetVoltages(data_reader);
00194 
00195     if (validResultsBasename == "")
00196     {
00197         validResultsBasename = rBaseResultsFilename;
00198     }
00199 
00200     ColumnDataReader valid_reader("heart/test/data/ionicmodels", validResultsBasename + "ValidData",
00201                                   false);
00202     std::vector<double> valid_times = valid_reader.GetValues("Time");
00203     std::vector<double> valid_voltages = GetVoltages(valid_reader);
00204 
00205     TS_ASSERT_EQUALS(times.size(), valid_times.size());
00206     for (unsigned i=0; i<valid_times.size(); i++)
00207     {
00208         TS_ASSERT_DELTA(times[i], valid_times[i], 1e-12);
00209         TS_ASSERT_DELTA(voltages[i], valid_voltages[i], tolerance);
00210     }
00211 }
00212 
00213 void CompareCellModelResults(std::string baseResultsFilename1, std::string baseResultsFilename2,
00214                              double tolerance, bool vOnly, std::string folderName)
00215 {
00216     // Compare 2 sets of results, e.g. from 2 different solvers for the same model.
00217     // If the time series differ, the finer resolution must be given first.
00218     ColumnDataReader data_reader1(folderName, baseResultsFilename1);
00219     std::vector<double> times1 = data_reader1.GetValues("Time");
00220     std::vector<double> voltages1 = GetVoltages(data_reader1);
00221     std::vector<double> calcium1;
00222     std::vector<double> h1;
00223 
00224     ColumnDataReader data_reader2(folderName, baseResultsFilename2);
00225     std::vector<double> times2 = data_reader2.GetValues("Time");
00226     std::vector<double> voltages2 = GetVoltages(data_reader2);
00227     std::vector<double> calcium2;
00228     std::vector<double> h2;
00229 
00230     if (!vOnly)
00231     {
00232         calcium1 = GetIntracellularCalcium(data_reader1);
00233         h1 = GetHGate(data_reader1);
00234         calcium2 = GetIntracellularCalcium(data_reader2);
00235         h2 = GetHGate(data_reader2);
00236     }
00237 
00238     TS_ASSERT(times1.size() >= times2.size());
00239     double last_v = voltages2[0];
00240     double tol = tolerance;
00241     for (unsigned i=0, j=0; i<times2.size(); i++)
00242     {
00243         // Find corresponding time index
00244         while (j<times1.size() && times1[j] < times2[i] - 1e-12)
00245         {
00246             j++;
00247         }
00248 
00249         // Set tolerance higher in upstroke
00250         if (fabs(voltages2[i] - last_v) > 0.05)
00251         {
00252             tol = tolerance * 25;
00253         }
00254         else
00255         {
00256             tol = tolerance;
00257         }
00258         last_v = voltages2[i];
00259 
00260         TS_ASSERT_DELTA(times1[j], times2[i], 1e-12);
00261         // adjust tol to data
00262         TS_ASSERT_DELTA(voltages1[j], voltages2[i], tol);
00263         if (!vOnly)
00264         {
00265             TS_ASSERT_DELTA(calcium1[j],  calcium2[i],  tol/100);
00266             TS_ASSERT_DELTA(h1[j],        h2[i],        tol/10);
00267         }
00268     }
00269 }
00270 
00271 
00272 #endif //_RUNANDCHECKIONICMODELS_HPP_

Generated on Tue May 31 14:31:41 2011 for Chaste by  doxygen 1.5.5