Chaste  Release::3.4
PropagationPropertiesCalculator.cpp
1 /*
2 
3 Copyright (c) 2005-2016, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #include "UblasIncludes.hpp"
37 #include "PropagationPropertiesCalculator.hpp"
38 #include "CellProperties.hpp"
39 #include "Exception.hpp"
40 #include <sstream>
41 #include "HeartEventHandler.hpp"
42 
44  const std::string voltageName)
45  : mpDataReader(pDataReader),
46  mVoltageName(voltageName),
47  mTimes(mpDataReader->GetUnlimitedDimensionValues()),
48  mCachedNodeGlobalIndex(UNSIGNED_UNSET)
49 {}
50 
52 {
53  // We don't own the data reader, so we don't destroy it.
54 }
55 
57 {
58  std::vector<double>& r_voltages = rGetCachedVoltages(globalNodeIndex);
59  CellProperties cell_props(r_voltages, mTimes);
60  return cell_props.GetLastMaxUpstrokeVelocity();
61 }
62 
63 std::vector<double> PropagationPropertiesCalculator::CalculateAllMaximumUpstrokeVelocities(unsigned globalNodeIndex, double threshold)
64 {
65  std::vector<double>& r_voltages = rGetCachedVoltages(globalNodeIndex);
66  CellProperties cell_props(r_voltages, mTimes, threshold);
67  return cell_props.GetMaxUpstrokeVelocities();
68 }
69 
70 std::vector<double> PropagationPropertiesCalculator::CalculateUpstrokeTimes(unsigned globalNodeIndex, double threshold)
71 {
72  std::vector<double>& r_voltages = rGetCachedVoltages(globalNodeIndex);
73  CellProperties cell_props(r_voltages, mTimes, threshold);
74  return cell_props.GetTimesAtMaxUpstrokeVelocity();
75 }
76 
78  unsigned globalNodeIndex)
79 {
80  if (percentage < 1.0 || percentage >= 100.0)
81  {
82  EXCEPTION("First argument of CalculateActionPotentialDuration() is expected to be a percentage");
83  }
84  std::vector<double>& r_voltages = rGetCachedVoltages(globalNodeIndex);
85  CellProperties cell_props(r_voltages, mTimes);
86  return cell_props.GetLastActionPotentialDuration(percentage);
87 }
88 
90  unsigned globalNodeIndex,
91  double threshold)
92 {
93  std::vector<double>& r_voltages = rGetCachedVoltages(globalNodeIndex);
94  CellProperties cell_props(r_voltages, mTimes, threshold);
95  return cell_props.GetAllActionPotentialDurations(percentage);
96 }
97 
99  const double percentage,
100  unsigned lowerNodeIndex,
101  unsigned upperNodeIndex,
102  double threshold)
103 {
104  std::vector<std::vector<double> > output_data;
105  output_data.reserve(upperNodeIndex-lowerNodeIndex+1);
106  unsigned num_nodes_per_data_block = 100; // number of nodes
107  unsigned num_complete_blocks = (upperNodeIndex-lowerNodeIndex) / num_nodes_per_data_block;
108  unsigned size_last_block = (upperNodeIndex-lowerNodeIndex) % num_nodes_per_data_block;
109 
110  for (unsigned block_num=0;
111  block_num<num_complete_blocks+1;
112  block_num++)
113  {
114  unsigned num_nodes_to_read;
115  if (block_num != num_complete_blocks)
116  {
117  num_nodes_to_read = num_nodes_per_data_block;
118  }
119  else
120  {
121  num_nodes_to_read = size_last_block;
122  }
123 
124  if (num_nodes_to_read > 0)
125  {
126  // Read a big block of data
127  unsigned low_node = lowerNodeIndex + block_num*num_nodes_per_data_block;
128  unsigned high_node = low_node + num_nodes_to_read;
129  std::vector<std::vector<double> > voltages = mpDataReader->GetVariableOverTimeOverMultipleNodes(mVoltageName, low_node, high_node);
130 
131  for (unsigned node_within_block=0;
132  node_within_block < num_nodes_to_read;
133  node_within_block++)
134  {
135  std::vector<double>& r_voltages = voltages[node_within_block];
136  CellProperties cell_props(r_voltages, mTimes, threshold);
137  std::vector<double> apds;
138  try
139  {
140  apds = cell_props.GetAllActionPotentialDurations(percentage);
141  assert(apds.size() != 0);
142  }
143  catch (Exception& e)
144  {
145  assert(e.GetShortMessage()=="No full action potential was recorded" ||
146  e.GetShortMessage()=="AP did not occur, never exceeded threshold voltage.");
147  apds.push_back(0);
148  assert(apds.size() == 1);
149  }
150  output_data.push_back(apds);
151  }
152  }
153  }
154  return output_data;
155 }
156 
158 {
159  std::vector<double>& r_voltages = rGetCachedVoltages(globalNodeIndex);
160  double max = -DBL_MAX;
161  for (unsigned i=0; i<r_voltages.size(); i++)
162  {
163  if (r_voltages[i]>max)
164  {
165  max = r_voltages[i];
166  }
167  }
168  return max;
169 }
170 
172  unsigned globalFarNodeIndex,
173  const double euclideanDistance)
174 {
175  double t_near = 0;
176  double t_far = 0;
177  std::vector<double>& r_near_voltages = rGetCachedVoltages(globalNearNodeIndex);
178  std::vector<double> far_voltages = mpDataReader->GetVariableOverTime(mVoltageName, globalFarNodeIndex);
179 
180  CellProperties near_cell_props(r_near_voltages, mTimes);
181  CellProperties far_cell_props(far_voltages, mTimes);
182 
183  //The size of each vector is the number of APs that reached that node
184  unsigned aps_near_node = near_cell_props.GetMaxUpstrokeVelocities().size();
185  unsigned aps_far_node = far_cell_props.GetMaxUpstrokeVelocities().size();
186 
187  //These should never be empty. If so, an exception should have been thrown in the GetMaxUpstrokeVelocities() method.
188  assert(aps_near_node > 0);
189  assert(aps_far_node > 0);
190 
191  //if the same number of APs reached both nodes, get the last one...
192  if (aps_near_node == aps_far_node)
193  {
194  t_near = near_cell_props.GetTimeAtLastMaxUpstrokeVelocity();
195  t_far = far_cell_props.GetTimeAtLastMaxUpstrokeVelocity();
196  }
197  //...otherwise get the one with the smallest value, which is the last AP to reach both nodes
198  //This prevents possible calculation of negative conduction velocities
199  //for repeated stimuli
200  else if (aps_near_node > aps_far_node)
201  {
202  t_near = near_cell_props.GetTimesAtMaxUpstrokeVelocity()[aps_far_node-1];
203  t_far = far_cell_props.GetTimesAtMaxUpstrokeVelocity()[aps_far_node-1];
204  }
205  else
206  {
207  t_near = near_cell_props.GetTimesAtMaxUpstrokeVelocity()[aps_near_node-1];
208  t_far = far_cell_props.GetTimesAtMaxUpstrokeVelocity()[aps_near_node-1];
209  }
210 
212  if ((globalNearNodeIndex == globalFarNodeIndex) || ( fabs(t_far - t_near) < 1e-8))
213  {
214  // globalNearNodeIndex and globalFarNodeIndex are the same node, preventing a 0/0
215  // or
216  // AP number i is happening at the same time at nodes globalNearNodeIndex and globalFarNodeIndex
217  return 0.0;
218  }
219  else
220  {
221  return euclideanDistance / (t_far - t_near);
222  }
223 
224 
225 }
226 
227 std::vector<double> PropagationPropertiesCalculator::CalculateAllConductionVelocities(unsigned globalNearNodeIndex,
228  unsigned globalFarNodeIndex,
229  const double euclideanDistance)
230 {
231  std::vector<double> conduction_velocities;
232 
233  std::vector<double> t_near;
234  std::vector<double> t_far;
235  unsigned number_of_aps = 0;
236 
237  std::vector<double>& r_near_voltages = rGetCachedVoltages(globalNearNodeIndex);
238  std::vector<double> far_voltages = mpDataReader->GetVariableOverTime(mVoltageName, globalFarNodeIndex);
239 
240  CellProperties near_cell_props(r_near_voltages, mTimes);
241  CellProperties far_cell_props(far_voltages, mTimes);
242 
243  t_near = near_cell_props.GetTimesAtMaxUpstrokeVelocity();
244  t_far = far_cell_props.GetTimesAtMaxUpstrokeVelocity();
245 
246  //exception should have been thrown within the GetTimesAtMaxUpstrokeVelocity method if the threshold is never reached
247  //and these vectors are empty
248  assert(t_near.size() !=0);
249  assert(t_far.size() !=0);
250 
251  //Check the node where the least number of aps is reached.
252  //We will calculate only where AP reached both nodes
253  if (t_near.size() > t_far.size())
254  {
255  number_of_aps = t_far.size();
256  }
257  else
258  {
259  number_of_aps = t_near.size();
260  }
261  //now fill the vector
262 
263  if (globalNearNodeIndex == globalFarNodeIndex)
264  {
265  // globalNearNodeIndex and globalFarNodeIndex are the same node, preventing a 0/0
266  for (unsigned i = 0 ; i < number_of_aps;i++)
267  {
268  conduction_velocities.push_back(0.0);
269  }
270  }
271  else
272  {
273  for (unsigned i = 0 ; i < number_of_aps;i++)
274  {
276  if ( fabs(t_far[i] - t_near[i]) < 1e-8)
277  {
278  // AP number i is happening at the same time at nodes globalNearNodeIndex and globalFarNodeIndex
279  conduction_velocities.push_back(0.0);
280  }
281  else
282  {
283  conduction_velocities.push_back(euclideanDistance / (t_far[i] - t_near[i]));
284  }
285  }
286  }
287 
288  return conduction_velocities;
289 }
290 
291 
292 std::vector<unsigned> PropagationPropertiesCalculator::CalculateAllAboveThresholdDepolarisations(unsigned globalNodeIndex, double threshold)
293 {
294  std::vector<double>& r_voltages = rGetCachedVoltages(globalNodeIndex);
295  CellProperties cell_props(r_voltages, mTimes, threshold);
297 }
298 
299 
301 {
302  std::vector<double>& r_voltages = rGetCachedVoltages(globalNodeIndex);
303  CellProperties cell_props(r_voltages, mTimes, threshold);
305 }
306 
307 
308 std::vector<double>& PropagationPropertiesCalculator::rGetCachedVoltages(unsigned globalNodeIndex)
309 {
310  if (globalNodeIndex != mCachedNodeGlobalIndex)
311  {
313  mCachedNodeGlobalIndex = globalNodeIndex;
314  }
315  return mCachedVoltages;
316 }
317 
319 {
320  mpDataReader = pDataReader;
321 }
322 
323 
std::vector< double > GetAllActionPotentialDurations(const double percentage)
double GetLastMaxUpstrokeVelocity()
void SetHdf5DataReader(Hdf5DataReader *pDataReader)
PropagationPropertiesCalculator(Hdf5DataReader *pDataReader, const std::string voltageName="V")
double CalculatePeakMembranePotential(unsigned globalNodeIndex)
#define EXCEPTION(message)
Definition: Exception.hpp:143
std::vector< double > CalculateAllMaximumUpstrokeVelocities(unsigned globalNodeIndex, double threshold)
std::vector< double > GetMaxUpstrokeVelocities()
unsigned GetNumberOfAboveThresholdDepolarisationsForLastAp()
std::vector< double > CalculateAllConductionVelocities(unsigned globalNearNodeIndex, unsigned globalFarNodeIndex, const double euclideanDistance)
double CalculateMaximumUpstrokeVelocity(unsigned globalNodeIndex)
std::vector< unsigned > CalculateAllAboveThresholdDepolarisations(unsigned globalNodeIndex, double threshold)
double CalculateActionPotentialDuration(const double percentage, unsigned globalNodeIndex)
std::vector< std::vector< double > > GetVariableOverTimeOverMultipleNodes(const std::string &rVariableName, unsigned lowerIndex, unsigned upperIndex)
std::vector< double > & rGetCachedVoltages(unsigned globalNodeIndex)
std::vector< double > GetTimesAtMaxUpstrokeVelocity()
double GetLastActionPotentialDuration(const double percentage)
std::vector< double > CalculateAllActionPotentialDurations(const double percentage, unsigned globalNodeIndex, double threshold)
const unsigned UNSIGNED_UNSET
Definition: Exception.hpp:52
double CalculateConductionVelocity(unsigned globalNearNodeIndex, unsigned globalFarNodeIndex, const double euclideanDistance)
std::vector< unsigned > GetNumberOfAboveThresholdDepolarisationsForAllAps()
std::vector< double > GetVariableOverTime(const std::string &rVariableName, unsigned nodeIndex)
std::vector< double > CalculateUpstrokeTimes(unsigned globalNodeIndex, double threshold)
unsigned CalculateAboveThresholdDepolarisationsForLastAp(unsigned globalNodeIndex, double threshold)
std::string GetShortMessage() const
Definition: Exception.cpp:92
std::vector< std::vector< double > > CalculateAllActionPotentialDurationsForNodeRange(const double percentage, unsigned lowerNodeIndex, unsigned upperNodeIndex, double threshold)
double GetTimeAtLastMaxUpstrokeVelocity()