Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
PostProcessingWriter.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, 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
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
37#include "HeartConfig.hpp"
38#include "PostProcessingWriter.hpp"
39#include "PetscTools.hpp"
40#include "OutputFileHandler.hpp"
41#include "DistanceMapCalculator.hpp"
42#include "PseudoEcgCalculator.hpp"
43#include "Version.hpp"
44#include "HeartEventHandler.hpp"
45#include "Hdf5DataWriter.hpp"
46#include "Hdf5ToMeshalyzerConverter.hpp"
47#include "Hdf5ToVtkConverter.hpp"
48
49#include <iostream>
50
51template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
53 const FileFinder& rDirectory,
54 const std::string& rHdf5FileName,
55 const std::string& rVoltageName,
56 hsize_t hdf5DataWriterChunkSize)
57 : mDirectory(rDirectory),
58 mHdf5File(rHdf5FileName),
59 mVoltageName(rVoltageName),
60 mrMesh(rMesh),
61 mHdf5DataWriterChunkSize(hdf5DataWriterChunkSize)
62{
63 mLo = mrMesh.GetDistributedVectorFactory()->GetLow();
64 mHi = mrMesh.GetDistributedVectorFactory()->GetHigh();
67 // Check that the hdf file was generated by simulations from (probably) the same mesh.
68 assert(mpDataReader->GetNumberOfRows() == mrMesh.GetNumNodes());
69}
70
71template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
73{
74 //Check that post-processing is really needed
75 assert(HeartConfig::Instance()->IsPostProcessingRequested());
76
77 // Please note that only the master processor should write to file.
78 // Each of the private methods called here takes care of checking.
79 if (HeartConfig::Instance()->IsApdMapsRequested())
80 {
81 std::vector<std::pair<double,double> > apd_maps;
83 for (unsigned i=0; i<apd_maps.size(); i++)
84 {
85 WriteApdMapFile(apd_maps[i].first, apd_maps[i].second);
86 }
87 }
88
89 if (HeartConfig::Instance()->IsUpstrokeTimeMapsRequested())
90 {
91 std::vector<double> upstroke_time_maps;
92 HeartConfig::Instance()->GetUpstrokeTimeMaps(upstroke_time_maps);
93 for (unsigned i=0; i<upstroke_time_maps.size(); i++)
94 {
95 WriteUpstrokeTimeMap(upstroke_time_maps[i]);
96 }
97 }
98
99 if (HeartConfig::Instance()->IsMaxUpstrokeVelocityMapRequested())
100 {
101 std::vector<double> upstroke_velocity_maps;
102 HeartConfig::Instance()->GetMaxUpstrokeVelocityMaps(upstroke_velocity_maps);
103 for (unsigned i=0; i<upstroke_velocity_maps.size(); i++)
104 {
105 WriteMaxUpstrokeVelocityMap(upstroke_velocity_maps[i]);
106 }
107 }
108
109 if (HeartConfig::Instance()->IsConductionVelocityMapsRequested())
110 {
111 std::vector<unsigned> conduction_velocity_maps;
112 HeartConfig::Instance()->GetConductionVelocityMaps(conduction_velocity_maps);
113
114 //get the mesh here
115 DistanceMapCalculator<ELEMENT_DIM, SPACE_DIM> dist_map_calculator(mrMesh);
116
117 for (unsigned i=0; i<conduction_velocity_maps.size(); i++)
118 {
119 std::vector<double> distance_map;
120 std::vector<unsigned> origin_surface;
121 origin_surface.push_back(conduction_velocity_maps[i]);
122 dist_map_calculator.ComputeDistanceMap(origin_surface, distance_map);
123 WriteConductionVelocityMap(conduction_velocity_maps[i], distance_map);
124 }
125 }
126
127 if (HeartConfig::Instance()->IsAnyNodalTimeTraceRequested())
128 {
129 std::vector<unsigned> requested_nodes;
131 WriteVariablesOverTimeAtNodes(requested_nodes);
132 }
133
134 if (HeartConfig::Instance()->IsPseudoEcgCalculationRequested())
135 {
136 std::vector<ChastePoint<SPACE_DIM> > electrodes;
138
141 delete mpDataReader;
142
143 for (unsigned i=0; i<electrodes.size(); i++)
144 {
146 electrodes[i],
147 mDirectory,
148 mHdf5File,
149 mVoltageName);
150 calculator.WritePseudoEcg();
151 }
152
154 mpDataReader = new Hdf5DataReader(mDirectory, mHdf5File);
155 mpCalculator->SetHdf5DataReader(mpDataReader);
156 }
157}
158
159template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
161{
162 delete mpDataReader;
163 delete mpCalculator;
164}
165
166template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
167void PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::WriteOutputDataToHdf5(const std::vector<std::vector<double> >& rDataPayload,
168 const std::string& rDatasetName,
169 const std::string& rDatasetUnit,
170 const std::string& rUnlimitedVariableName,
171 const std::string& rUnlimitedVariableUnit)
172{
173 DistributedVectorFactory* p_factory = mrMesh.GetDistributedVectorFactory();
175 Hdf5DataWriter writer(*p_factory,
176 mDirectory.GetRelativePath(test_output), // Path relative to CHASTE_TEST_OUTPUT
177 mHdf5File,
178 false, // to wiping
179 true, // to extending
180 rDatasetName); // dataset name
181
182 /* Probe define mode. We asked to extend, so if the writer is currently in define mode it means the
183 * dataset doesn't exist yet and needs creating. If it is NOT in define mode, it means the dataset
184 * exists, so we'll empty it and calculate new postprocessing data. */
185 int apd_id;
186 if (writer.IsInDefineMode())
187 {
188 apd_id = writer.DefineVariable(rDatasetName, rDatasetUnit);
189 writer.DefineFixedDimension(mrMesh.GetNumNodes());
190 writer.DefineUnlimitedDimension(rUnlimitedVariableName, rUnlimitedVariableUnit);
191 if (mHdf5DataWriterChunkSize>0u)
192 {
193 /* Pass target chunk size through to writer. (We don't do the
194 * alignment one as well because that can only be done for a new
195 * file and PostProcessingWriter can only add to an existing file.)
196 */
197 writer.SetTargetChunkSize(mHdf5DataWriterChunkSize);
198 }
199 writer.EndDefineMode();
200 }
201 else
202 {
203 apd_id = writer.GetVariableByName(rDatasetName);
204 writer.EmptyDataset();
205 }
206
207 //Determine the maximum number of paces
208 unsigned local_max_paces = 0u;
209 for (unsigned node_index = 0; node_index < rDataPayload.size(); ++node_index)
210 {
211 if (rDataPayload[node_index].size() > local_max_paces)
212 {
213 local_max_paces = rDataPayload[node_index].size();
214 }
215 }
216
217 unsigned max_paces = 0u;
218 MPI_Allreduce(&local_max_paces, &max_paces, 1, MPI_UNSIGNED, MPI_MAX, PETSC_COMM_WORLD);
219
220 for (unsigned pace_idx = 0; pace_idx < max_paces; pace_idx++)
221 {
222 Vec apd_vec = p_factory->CreateVec();
223 DistributedVector distributed_vector = p_factory->CreateDistributedVector(apd_vec);
224 for (DistributedVector::Iterator index = distributed_vector.Begin();
225 index!= distributed_vector.End();
226 ++index)
227 {
228 unsigned node_idx = index.Local;
229 // pad with -999 if no pace defined at this node
230 if (pace_idx < rDataPayload[node_idx].size() )
231 {
232 distributed_vector[index] = rDataPayload[node_idx][pace_idx];
233 }
234 else
235 {
236 distributed_vector[index] = -999.0;
237 }
238 }
239 writer.PutVector(apd_id, apd_vec);
240 PetscTools::Destroy(apd_vec);
241 writer.PutUnlimitedVariable(pace_idx);
243 }
244 writer.Close();
245}
246
247template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
248void PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::WriteApdMapFile(double repolarisationPercentage, double threshold)
249{
250 std::vector<std::vector<double> > local_output_data = mpCalculator->CalculateAllActionPotentialDurationsForNodeRange(repolarisationPercentage, mLo, mHi, threshold);
251
252 // HDF5 shouldn't have minus signs in the data names..
253 std::stringstream hdf5_dataset_name;
254 hdf5_dataset_name << "Apd_" << repolarisationPercentage;
255
256 WriteOutputDataToHdf5(local_output_data,
257 hdf5_dataset_name.str() + ConvertToHdf5FriendlyString(threshold) + "_Map",
258 "msec");
259}
260
261template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
263{
264 std::vector<std::vector<double> > output_data;
265 //Fill in data
266 for (unsigned node_index = mLo; node_index < mHi; node_index++)
267 {
268 std::vector<double> upstroke_times;
269 try
270 {
271 upstroke_times = mpCalculator->CalculateUpstrokeTimes(node_index, threshold);
272 assert(upstroke_times.size() != 0);
273 }
274 catch(Exception&)
275 {
276 upstroke_times.push_back(0);
277 assert(upstroke_times.size() == 1);
278 }
279 output_data.push_back(upstroke_times);
280 }
281
282 WriteOutputDataToHdf5(output_data,
283 "UpstrokeTimeMap" + ConvertToHdf5FriendlyString(threshold),
284 "msec");
285}
286
287template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
289{
290 std::vector<std::vector<double> > output_data;
291 //Fill in data
292 for (unsigned node_index = mLo; node_index < mHi; node_index++)
293 {
294 std::vector<double> upstroke_velocities;
295 try
296 {
297 upstroke_velocities = mpCalculator->CalculateAllMaximumUpstrokeVelocities(node_index, threshold);
298 assert(upstroke_velocities.size() != 0);
299 }
300 catch(Exception&)
301 {
302 upstroke_velocities.push_back(0);
303 assert(upstroke_velocities.size() ==1);
304 }
305 output_data.push_back(upstroke_velocities);
306 }
307
308 WriteOutputDataToHdf5(output_data,
309 "MaxUpstrokeVelocityMap" + ConvertToHdf5FriendlyString(threshold),
310 "mV_per_msec");
311}
312
313template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
314void PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::WriteConductionVelocityMap(unsigned originNode, std::vector<double> distancesFromOriginNode)
315{
316 //Fill in data
317 std::vector<std::vector<double> > output_data;
318 for (unsigned dest_node = mLo; dest_node < mHi; dest_node++)
319 {
320 std::vector<double> conduction_velocities;
321 try
322 {
323 conduction_velocities = mpCalculator->CalculateAllConductionVelocities(originNode, dest_node, distancesFromOriginNode[dest_node]);
324 assert(conduction_velocities.size() != 0);
325 }
326 catch(Exception&)
327 {
328 conduction_velocities.push_back(0);
329 assert(conduction_velocities.size() == 1);
330 }
331 output_data.push_back(conduction_velocities);
332 }
333 std::stringstream filename_stream;
334 filename_stream << "ConductionVelocityFromNode" << originNode;
335
336 WriteOutputDataToHdf5(output_data, filename_stream.str(), "cm_per_msec");
337}
338
339template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
341{
342 std::vector<std::vector<double> > output_data;
343
344 //Fill in data
345 for (unsigned node_index = mLo; node_index < mHi; node_index++)
346 {
347 std::vector<double> upstroke_velocities;
348 std::vector<unsigned> above_threshold_depolarisations;
349 std::vector<double> output_item;
350 bool no_upstroke_occurred = false;
351
352 try
353 {
354 upstroke_velocities = mpCalculator->CalculateAllMaximumUpstrokeVelocities(node_index, threshold);
355 assert(upstroke_velocities.size() != 0);
356 }
357 catch(Exception&)
358 {
359 upstroke_velocities.push_back(0);
360 assert(upstroke_velocities.size() == 1);
361 no_upstroke_occurred = true;
362 }
363 // This method won't throw any exception, so there is no need to put it into the try/catch:
364 above_threshold_depolarisations = mpCalculator->CalculateAllAboveThresholdDepolarisations(node_index, threshold);
365
366 // Count the total above threshold depolarisations
367 unsigned total_number_of_above_threshold_depolarisations = 0;
368 for (unsigned ead_index = 0; ead_index< above_threshold_depolarisations.size();ead_index++)
369 {
370 total_number_of_above_threshold_depolarisations = total_number_of_above_threshold_depolarisations + above_threshold_depolarisations[ead_index];
371 }
372
373 // For this item, push back the number of upstrokes...
374 if (no_upstroke_occurred)
375 {
376 output_item.push_back(0);
377 }
378 else
379 {
380 output_item.push_back((double)upstroke_velocities.size());
381 }
382 //... and the number of above threshold depolarisations
383 output_item.push_back((double) total_number_of_above_threshold_depolarisations);
384
385 output_data.push_back(output_item);
386 }
387
388 // we just use meshalyzer format so that something is generated in simple column format for use with gnuplot etc.
389 WriteGenericFileToMeshalyzer(output_data, "", "AboveThresholdDepolarisations" + ConvertToHdf5FriendlyString(threshold) +
390 ".dat");
391}
392
393template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
395{
396 std::stringstream dataset_stream;
397 std::string extra_message = "_";
398 if (threshold < 0.0)
399 {
400 extra_message += "minus_";
401 threshold = -threshold;
402 }
403
404 // if threshold is a decimal eg: because using phenomenological model
405 if (threshold - floor(threshold) > 1e-8)
406 {
407 // give the answer to 2dp
408 dataset_stream << extra_message << floor(threshold) << "pt" << floor(threshold*100)-(floor(threshold)*100);
409 }
410 else
411 {
412 dataset_stream << extra_message << threshold;
413 }
414
415 return dataset_stream.str();
416}
417
418template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
420{
421 std::vector<std::string> variable_names = mpDataReader->GetVariableNames();
422
423 //we will write one file per variable in the hdf5 file
424 for (unsigned name_index=0; name_index < variable_names.size(); name_index++)
425 {
426 std::vector<std::vector<double> > output_data;
427 if (PetscTools::AmMaster())//only master process fills the data structure
428 {
429 //allocate memory: NXM matrix where N = number of time steps and M number of requested nodes
430 output_data.resize( mpDataReader->GetUnlimitedDimensionValues().size() );
431 for (unsigned j = 0; j < mpDataReader->GetUnlimitedDimensionValues().size(); j++)
432 {
433 output_data[j].resize(rNodeIndices.size());
434 }
435
436 for (unsigned requested_index = 0; requested_index < rNodeIndices.size(); requested_index++)
437 {
438 unsigned node_index = rNodeIndices[requested_index];
439
440 // Handle permutation, if any
441 if ((mrMesh.rGetNodePermutation().size() != 0) &&
443 {
444 node_index = mrMesh.rGetNodePermutation()[ rNodeIndices[requested_index] ];
445 }
446
447 // Grab the data from the hdf5 file.
448 std::vector<double> time_series = mpDataReader->GetVariableOverTime(variable_names[name_index], node_index);
449 assert ( time_series.size() == mpDataReader->GetUnlimitedDimensionValues().size());
450
451 // Fill the output_data data structure
452 for (unsigned time_step = 0; time_step < time_series.size(); time_step++)
453 {
454 output_data[time_step][requested_index] = time_series[time_step];
455 }
456 }
457 }
458 std::stringstream filename_stream;
459 filename_stream << "NodalTraces_" << variable_names[name_index] << ".dat";
460 // we just use meshalyzer format so that something is generated in simple column format for use with gnuplot etc.
461 WriteGenericFileToMeshalyzer(output_data, "", filename_stream.str());
462 }
463}
464
465template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
466void PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::WriteGenericFileToMeshalyzer(std::vector<std::vector<double> >& rDataPayload, const std::string& rFolder, const std::string& rFileName)
467{
468 OutputFileHandler output_file_handler(HeartConfig::Instance()->GetOutputDirectory() + "/" + rFolder, false);
470 {
471 out_stream p_file=out_stream(NULL);
472 //Open file
474 {
475 // Open the file for the first time
476 p_file = output_file_handler.OpenOutputFile(rFileName);
477 }
478 else
479 {
480 // Append data to the existing file opened by master
481 p_file = output_file_handler.OpenOutputFile(rFileName, std::ios::app);
482 }
483 // Write data
484 for (unsigned line_number=0; line_number<rDataPayload.size(); line_number++)
485 {
486 for (unsigned i = 0; i < rDataPayload[line_number].size(); i++)
487 {
488 *p_file << rDataPayload[line_number][i] << "\t";
489 }
490 *p_file << std::endl;
491 }
492
493 // Last processor appends comment line
495 {
496 std::string comment = "# " + ChasteBuildInfo::GetProvenanceString();
497 *p_file << comment;
498 }
499 p_file->close();
500 }
501 //There's a barrier included here: Process i+1 waits for process i to close the file
503}
504
505// Explicit instantiation
506template class PostProcessingWriter<1,1>;
507template class PostProcessingWriter<1,2>;
508template class PostProcessingWriter<2,2>;
509template class PostProcessingWriter<1,3>;
510template class PostProcessingWriter<3,3>;
static std::string GetProvenanceString()
void ComputeDistanceMap(const std::vector< unsigned > &rSourceNodeIndices, std::vector< double > &rNodeDistances)
DistributedVector CreateDistributedVector(Vec vec, bool readOnly=false)
unsigned GetNumberOfRows()
void SetTargetChunkSize(hsize_t targetSize)
void DefineUnlimitedDimension(const std::string &rVariableName, const std::string &rVariableUnits, unsigned estimatedLength=1)
void PutUnlimitedVariable(double value)
void AdvanceAlongUnlimitedDimension()
int GetVariableByName(const std::string &rVariableName)
void DefineFixedDimension(long dimensionSize)
int DefineVariable(const std::string &rVariableName, const std::string &rVariableUnits)
virtual void EndDefineMode()
void PutVector(int variableID, Vec petscVector)
bool GetOutputUsingOriginalNodeOrdering()
void GetNodalTimeTraceRequested(std::vector< unsigned > &rRequestedNodes) const
void GetApdMaps(std::vector< std::pair< double, double > > &rApdMaps) const
void GetConductionVelocityMaps(std::vector< unsigned > &rConductionVelocityMaps) const
void GetPseudoEcgElectrodePositions(std::vector< ChastePoint< SPACE_DIM > > &rPseudoEcgElectrodePositions) const
void GetUpstrokeTimeMaps(std::vector< double > &rUpstrokeTimeMaps) const
void GetMaxUpstrokeVelocityMaps(std::vector< double > &rUpstrokeVelocityMaps) const
static HeartConfig * Instance()
out_stream OpenOutputFile(const std::string &rFileName, std::ios_base::openmode mode=std::ios::out|std::ios::trunc) const
static void Destroy(Vec &rVec)
static bool AmMaster()
static bool AmTopMost()
static void EndRoundRobin()
static void BeginRoundRobin()
void WriteMaxUpstrokeVelocityMap(double threshold)
void WriteConductionVelocityMap(unsigned originNode, std::vector< double > distancesFromOriginNode)
void WriteVariablesOverTimeAtNodes(std::vector< unsigned > &rNodeIndices)
void WriteUpstrokeTimeMap(double threshold)
void WriteGenericFileToMeshalyzer(std::vector< std::vector< double > > &rDataPayload, const std::string &rFolder, const std::string &rFileName)
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > & mrMesh
std::string ConvertToHdf5FriendlyString(double threshold)
void WriteApdMapFile(double repolarisationPercentage, double threshold)
void WriteOutputDataToHdf5(const std::vector< std::vector< double > > &rDataPayload, const std::string &rDatasetName, const std::string &rDatasetUnit, const std::string &rUnlimitedVariableName="PaceNumber", const std::string &rUnlimitedVariableUnit="dimensionless")
PostProcessingWriter(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rMesh, const FileFinder &rDirectory, const std::string &rHdf5FileName, const std::string &rVoltageName="V", hsize_t hdf5DataWriterChunkSize=0)
void WriteAboveThresholdDepolarisationFile(double threshold)
PropagationPropertiesCalculator * mpCalculator