Line data Source code
1 : #ifdef CRPROPA_HAVE_HDF5
2 :
3 : #include "crpropa/module/HDF5Output.h"
4 : #include "crpropa/Version.h"
5 : #include "crpropa/Random.h"
6 : #include "kiss/logger.h"
7 :
8 : #include <hdf5.h>
9 : #include <cstring>
10 :
11 : const hsize_t RANK = 1;
12 : const hsize_t BUFFER_SIZE = 1024 * 16;
13 :
14 : namespace crpropa {
15 :
16 : // map variant types to H5T_NATIVE
17 0 : hid_t variantTypeToH5T_NATIVE(Variant::Type type) {
18 : if (type == Variant::TYPE_INT64)
19 0 : return H5T_NATIVE_INT64;
20 : else if(type == Variant::TYPE_BOOL)
21 0 : return H5T_NATIVE_HBOOL;
22 : else if(type == Variant::TYPE_CHAR)
23 0 : return H5T_NATIVE_CHAR;
24 : else if(type == Variant::TYPE_UCHAR)
25 0 : return H5T_NATIVE_UCHAR;
26 : else if(type == Variant::TYPE_INT16)
27 0 : return H5T_NATIVE_INT16;
28 : else if(type == Variant::TYPE_UINT16)
29 0 : return H5T_NATIVE_UINT16;
30 : else if(type == Variant::TYPE_INT32)
31 0 : return H5T_NATIVE_INT32;
32 : else if(type == Variant::TYPE_UINT32)
33 0 : return H5T_NATIVE_UINT32;
34 : else if(type == Variant::TYPE_INT64)
35 : return H5T_NATIVE_INT64;
36 : else if(type == Variant::TYPE_UINT64)
37 0 : return H5T_NATIVE_UINT64;
38 : else if(type == Variant::TYPE_FLOAT)
39 0 : return H5T_NATIVE_FLOAT;
40 : else if(type == Variant::TYPE_DOUBLE)
41 0 : return H5T_NATIVE_DOUBLE;
42 : else if(type == Variant::TYPE_STRING)
43 0 : return H5T_C_S1;
44 : else
45 : {
46 0 : KISS_LOG_ERROR << "variantTypeToH5T_NATIVE:: Type: " << Variant::getTypeName(type) << " unknown.";
47 0 : throw std::runtime_error("No matching HDF type for Variant type");
48 : }
49 : }
50 :
51 1 : HDF5Output::HDF5Output() : Output(), filename(), file(-1), sid(-1), dset(-1), dataspace(-1), candidatesSinceFlush(0), flushLimit(std::numeric_limits<unsigned int>::max()) {
52 1 : }
53 :
54 0 : HDF5Output::HDF5Output(const std::string& filename) : Output(), filename(filename), file(-1), sid(-1), dset(-1), dataspace(-1), candidatesSinceFlush(0), flushLimit(std::numeric_limits<unsigned int>::max()) {
55 0 : }
56 :
57 0 : HDF5Output::HDF5Output(const std::string& filename, OutputType outputtype) : Output(outputtype), filename(filename), file(-1), sid(-1), dset(-1), dataspace(-1), candidatesSinceFlush(0), flushLimit(std::numeric_limits<unsigned int>::max()) {
58 : outputtype = outputtype;
59 0 : }
60 :
61 1 : HDF5Output::~HDF5Output() {
62 1 : close();
63 1 : }
64 :
65 0 : herr_t HDF5Output::insertStringAttribute(const std::string &key, const std::string &value){
66 : hid_t strtype, attr_space, version_attr;
67 0 : hsize_t dims = 0;
68 : herr_t status;
69 :
70 0 : strtype = H5Tcopy(H5T_C_S1);
71 0 : status = H5Tset_size(strtype, value.size());
72 :
73 0 : attr_space = H5Screate_simple(0, &dims, NULL);
74 0 : version_attr = H5Acreate2(dset, key.c_str(), strtype, attr_space, H5P_DEFAULT, H5P_DEFAULT);
75 0 : status = H5Awrite(version_attr, strtype, value.c_str());
76 0 : status = H5Aclose(version_attr);
77 0 : status = H5Sclose(attr_space);
78 :
79 0 : return status;
80 : }
81 :
82 0 : herr_t HDF5Output::insertDoubleAttribute(const std::string &key, const double &value){
83 : hid_t type, attr_space, version_attr;
84 0 : hsize_t dims = 0;
85 : herr_t status;
86 :
87 0 : type = H5Tcopy(H5T_NATIVE_DOUBLE);
88 :
89 0 : attr_space = H5Screate_simple(0, &dims, NULL);
90 0 : version_attr = H5Acreate2(dset, key.c_str(), type, attr_space, H5P_DEFAULT, H5P_DEFAULT);
91 0 : status = H5Awrite(version_attr, type, &value);
92 0 : status = H5Aclose(version_attr);
93 0 : status = H5Sclose(attr_space);
94 :
95 0 : return status;
96 : }
97 :
98 :
99 :
100 1 : void HDF5Output::open(const std::string& filename) {
101 1 : file = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
102 1 : if (file < 0)
103 2 : throw std::runtime_error(std::string("Cannot create file: ") + filename);
104 :
105 :
106 0 : sid = H5Tcreate(H5T_COMPOUND, sizeof(OutputRow));
107 0 : if (fields.test(TrajectoryLengthColumn))
108 0 : H5Tinsert(sid, "D", HOFFSET(OutputRow, D), H5T_NATIVE_DOUBLE);
109 0 : if (fields.test(TimeColumn))
110 0 : H5Tinsert(sid, "time", HOFFSET(OutputRow, time), H5T_NATIVE_DOUBLE);
111 0 : if (fields.test(RedshiftColumn))
112 0 : H5Tinsert(sid, "z", HOFFSET(OutputRow, z), H5T_NATIVE_DOUBLE);
113 0 : if (fields.test(SerialNumberColumn))
114 0 : H5Tinsert(sid, "SN", HOFFSET(OutputRow, SN), H5T_NATIVE_UINT64);
115 0 : if (fields.test(CurrentIdColumn))
116 0 : H5Tinsert(sid, "ID", HOFFSET(OutputRow, ID), H5T_NATIVE_INT32);
117 0 : if (fields.test(CurrentEnergyColumn))
118 0 : H5Tinsert(sid, "E", HOFFSET(OutputRow, E), H5T_NATIVE_DOUBLE);
119 0 : if (fields.test(CurrentPositionColumn) && oneDimensional)
120 0 : H5Tinsert(sid, "X", HOFFSET(OutputRow, X), H5T_NATIVE_DOUBLE);
121 0 : if (fields.test(CurrentPositionColumn) && not oneDimensional) {
122 0 : H5Tinsert(sid, "X", HOFFSET(OutputRow, X), H5T_NATIVE_DOUBLE);
123 0 : H5Tinsert(sid, "Y", HOFFSET(OutputRow, Y), H5T_NATIVE_DOUBLE);
124 0 : H5Tinsert(sid, "Z", HOFFSET(OutputRow, Z), H5T_NATIVE_DOUBLE);
125 : }
126 0 : if (fields.test(CurrentDirectionColumn) && not oneDimensional) {
127 0 : H5Tinsert(sid, "Px", HOFFSET(OutputRow, Px), H5T_NATIVE_DOUBLE);
128 0 : H5Tinsert(sid, "Py", HOFFSET(OutputRow, Py), H5T_NATIVE_DOUBLE);
129 0 : H5Tinsert(sid, "Pz", HOFFSET(OutputRow, Pz), H5T_NATIVE_DOUBLE);
130 : }
131 0 : if (fields.test(SerialNumberColumn))
132 0 : H5Tinsert(sid, "SN0", HOFFSET(OutputRow, SN0), H5T_NATIVE_UINT64);
133 0 : if (fields.test(SourceIdColumn))
134 0 : H5Tinsert(sid, "ID0", HOFFSET(OutputRow, ID0), H5T_NATIVE_INT32);
135 0 : if (fields.test(SourceEnergyColumn))
136 0 : H5Tinsert(sid, "E0", HOFFSET(OutputRow, E0), H5T_NATIVE_DOUBLE);
137 0 : if (fields.test(SourcePositionColumn) && oneDimensional)
138 0 : H5Tinsert(sid, "X0", HOFFSET(OutputRow, X0), H5T_NATIVE_DOUBLE);
139 0 : if (fields.test(SourcePositionColumn) && not oneDimensional){
140 0 : H5Tinsert(sid, "X0", HOFFSET(OutputRow, X0), H5T_NATIVE_DOUBLE);
141 0 : H5Tinsert(sid, "Y0", HOFFSET(OutputRow, Y0), H5T_NATIVE_DOUBLE);
142 0 : H5Tinsert(sid, "Z0", HOFFSET(OutputRow, Z0), H5T_NATIVE_DOUBLE);
143 : }
144 0 : if (fields.test(SourceDirectionColumn) && not oneDimensional) {
145 0 : H5Tinsert(sid, "P0x", HOFFSET(OutputRow, P0x), H5T_NATIVE_DOUBLE);
146 0 : H5Tinsert(sid, "P0y", HOFFSET(OutputRow, P0y), H5T_NATIVE_DOUBLE);
147 0 : H5Tinsert(sid, "P0z", HOFFSET(OutputRow, P0z), H5T_NATIVE_DOUBLE);
148 : }
149 0 : if (fields.test(SerialNumberColumn))
150 0 : H5Tinsert(sid, "SN1", HOFFSET(OutputRow, SN1), H5T_NATIVE_UINT64);
151 0 : if (fields.test(CreatedIdColumn))
152 0 : H5Tinsert(sid, "ID1", HOFFSET(OutputRow, ID1), H5T_NATIVE_INT32);
153 0 : if (fields.test(CreatedEnergyColumn))
154 0 : H5Tinsert(sid, "E1", HOFFSET(OutputRow, E1), H5T_NATIVE_DOUBLE);
155 0 : if (fields.test(CreatedPositionColumn) && oneDimensional)
156 0 : H5Tinsert(sid, "X1", HOFFSET(OutputRow, X1), H5T_NATIVE_DOUBLE);
157 0 : if (fields.test(CreatedPositionColumn) && not oneDimensional) {
158 0 : H5Tinsert(sid, "X1", HOFFSET(OutputRow, X1), H5T_NATIVE_DOUBLE);
159 0 : H5Tinsert(sid, "Y1", HOFFSET(OutputRow, Y1), H5T_NATIVE_DOUBLE);
160 0 : H5Tinsert(sid, "Z1", HOFFSET(OutputRow, Z1), H5T_NATIVE_DOUBLE);
161 : }
162 0 : if (fields.test(CreatedDirectionColumn) && not oneDimensional) {
163 0 : H5Tinsert(sid, "P1x", HOFFSET(OutputRow, P1x), H5T_NATIVE_DOUBLE);
164 0 : H5Tinsert(sid, "P1y", HOFFSET(OutputRow, P1y), H5T_NATIVE_DOUBLE);
165 0 : H5Tinsert(sid, "P1z", HOFFSET(OutputRow, P1z), H5T_NATIVE_DOUBLE);
166 : }
167 0 : if (fields.test(WeightColumn))
168 0 : H5Tinsert(sid, "W", HOFFSET(OutputRow, weight), H5T_NATIVE_DOUBLE);
169 :
170 0 : if (fields.test(CandidateTagColumn))
171 0 : H5Tinsert(sid, "tag", HOFFSET(OutputRow, tag), H5T_C_S1);
172 :
173 0 : size_t pos = 0;
174 0 : for(std::vector<Output::Property>::const_iterator iter = properties.begin();
175 0 : iter != properties.end(); ++iter)
176 : {
177 0 : hid_t type = variantTypeToH5T_NATIVE((*iter).defaultValue.getType());
178 0 : if (type == H5T_C_S1)
179 : { // set size of string field to size of default value!
180 0 : type = H5Tcopy(H5T_C_S1);
181 0 : H5Tset_size(type, (*iter).defaultValue.toString().size());
182 : }
183 :
184 0 : H5Tinsert(sid, (*iter).name.c_str(), HOFFSET(OutputRow, propertyBuffer) + pos, type);
185 0 : pos += (*iter).defaultValue.getSize();
186 : }
187 0 : if (pos >= propertyBufferSize)
188 : {
189 0 : KISS_LOG_ERROR << "Using " << pos << " bytes for properties output. Maximum is " << propertyBufferSize << " bytes.";
190 0 : throw std::runtime_error("Size of property buffer exceeded");
191 : }
192 :
193 : // chunked prop
194 0 : hid_t plist = H5Pcreate(H5P_DATASET_CREATE);
195 0 : H5Pset_layout(plist, H5D_CHUNKED);
196 0 : hsize_t chunk_dims[RANK] = {BUFFER_SIZE};
197 0 : H5Pset_chunk(plist, RANK, chunk_dims);
198 0 : H5Pset_deflate(plist, 5);
199 :
200 0 : hsize_t dims[RANK] = {0};
201 0 : hsize_t max_dims[RANK] = {H5S_UNLIMITED};
202 0 : dataspace = H5Screate_simple(RANK, dims, max_dims);
203 :
204 0 : dset = H5Dcreate2(file, "CRPROPA3", sid, dataspace, H5P_DEFAULT, plist, H5P_DEFAULT);
205 :
206 0 : insertStringAttribute("OutputType", outputName);
207 0 : insertStringAttribute("Version", g_GIT_DESC);
208 0 : insertDoubleAttribute("LengthScale", this->lengthScale);
209 0 : insertDoubleAttribute("TimeScale", this->timeScale);
210 0 : insertDoubleAttribute("EnergyScale", this->energyScale);
211 :
212 : // add ranom seeds
213 0 : std::vector< std::vector<uint32_t> > seeds = Random::getSeedThreads();
214 0 : for (size_t i = 0; i < seeds.size(); i++)
215 : {
216 : hid_t type, attr_space, version_attr;
217 : herr_t status;
218 0 : hsize_t dims[] = {1, 0};
219 0 : dims[1] = seeds[i].size();
220 :
221 0 : type = H5Tarray_create(H5T_NATIVE_ULONG, 2, dims);
222 :
223 0 : attr_space = H5Screate_simple(0, dims, NULL);
224 : char nameBuffer[256];
225 : snprintf(nameBuffer, 256, "SEED_%03lu", i);
226 0 : KISS_LOG_DEBUG << "Creating HDF5 attribute: " << nameBuffer << " with dimensions " << dims[0] << "x" << dims[1] ;
227 :
228 0 : version_attr = H5Acreate2(dset, nameBuffer, type, attr_space, H5P_DEFAULT, H5P_DEFAULT);
229 0 : status = H5Awrite(version_attr, type, &seeds[i][0]);
230 0 : status = H5Aclose(version_attr);
231 0 : status = H5Sclose(attr_space);
232 :
233 : }
234 :
235 :
236 0 : H5Pclose(plist);
237 :
238 0 : buffer.reserve(BUFFER_SIZE);
239 0 : time(&lastFlush);
240 0 : }
241 :
242 1 : void HDF5Output::close() {
243 1 : if (file >= 0) {
244 0 : flush();
245 0 : H5Dclose(dset);
246 0 : H5Tclose(sid);
247 0 : H5Sclose(dataspace);
248 0 : H5Fclose(file);
249 0 : file = -1;
250 : }
251 1 : }
252 :
253 0 : void HDF5Output::process(Candidate* candidate) const {
254 0 : #pragma omp critical(HDFOutput)
255 : {
256 0 : if (file == -1)
257 : // This is ugly, but necesary as otherwise the user has to manually open the
258 : // file before processing the first candidate
259 0 : const_cast<HDF5Output*>(this)->open(filename);
260 : }
261 :
262 : OutputRow r;
263 0 : r.D = candidate->getTrajectoryLength() / lengthScale;
264 0 : r.time = candidate->getTime() / timeScale;
265 0 : r.z = candidate->getRedshift();
266 :
267 0 : r.SN = candidate->getSerialNumber();
268 0 : r.ID = candidate->current.getId();
269 0 : r.E = candidate->current.getEnergy() / energyScale;
270 0 : Vector3d v = candidate->current.getPosition() / lengthScale;
271 0 : r.X = v.x;
272 0 : r.Y = v.y;
273 0 : r.Z = v.z;
274 0 : v = candidate->current.getDirection();
275 0 : r.Px = v.x;
276 0 : r.Py = v.y;
277 0 : r.Pz = v.z;
278 :
279 0 : r.SN0 = candidate->getSourceSerialNumber();
280 0 : r.ID0 = candidate->source.getId();
281 0 : r.E0 = candidate->source.getEnergy() / energyScale;
282 0 : v = candidate->source.getPosition() / lengthScale;
283 0 : r.X0 = v.x;
284 0 : r.Y0 = v.y;
285 0 : r.Z0 = v.z;
286 0 : v = candidate->source.getDirection();
287 0 : r.P0x = v.x;
288 0 : r.P0y = v.y;
289 0 : r.P0z = v.z;
290 :
291 0 : r.SN1 = candidate->getCreatedSerialNumber();
292 0 : r.ID1 = candidate->created.getId();
293 0 : r.E1 = candidate->created.getEnergy() / energyScale;
294 0 : v = candidate->created.getPosition() / lengthScale;
295 0 : r.X1 = v.x;
296 0 : r.Y1 = v.y;
297 0 : r.Z1 = v.z;
298 0 : v = candidate->created.getDirection();
299 0 : r.P1x = v.x;
300 0 : r.P1y = v.y;
301 0 : r.P1z = v.z;
302 :
303 0 : r.weight= candidate->getWeight();
304 :
305 0 : r.tag = candidate->getTagOrigin();
306 :
307 : size_t pos = 0;
308 0 : for(std::vector<Output::Property>::const_iterator iter = properties.begin();
309 0 : iter != properties.end(); ++iter)
310 : {
311 0 : Variant v;
312 0 : if (candidate->hasProperty((*iter).name))
313 : {
314 0 : v = candidate->getProperty((*iter).name);
315 : }
316 : else
317 : {
318 0 : v = (*iter).defaultValue;
319 : }
320 0 : pos += v.copyToBuffer(&r.propertyBuffer[pos]);
321 0 : }
322 :
323 0 : #pragma omp critical(HDFOutput)
324 : {
325 0 : const_cast<HDF5Output*>(this)->candidatesSinceFlush++;
326 0 : Output::process(candidate);
327 :
328 0 : buffer.push_back(r);
329 :
330 :
331 0 : if (buffer.size() >= buffer.capacity())
332 : {
333 0 : KISS_LOG_DEBUG << "HDF5Output: Flush due to buffer capacity exceeded";
334 0 : flush();
335 : }
336 0 : else if (candidatesSinceFlush >= flushLimit)
337 : {
338 0 : KISS_LOG_DEBUG << "HDF5Output: Flush due to number of candidates";
339 0 : flush();
340 : }
341 0 : else if (difftime(time(NULL), lastFlush) > 60*10)
342 : {
343 0 : KISS_LOG_DEBUG << "HDF5Output: Flush due to time exceeded";
344 0 : flush();
345 : }
346 : }
347 0 : }
348 :
349 0 : void HDF5Output::flush() const {
350 0 : const_cast<HDF5Output*>(this)->lastFlush = time(NULL);
351 0 : const_cast<HDF5Output*>(this)->candidatesSinceFlush = 0;
352 :
353 : hsize_t n = buffer.size();
354 :
355 0 : if (n == 0)
356 0 : return;
357 :
358 0 : hid_t file_space = H5Dget_space(dset);
359 0 : hsize_t count = H5Sget_simple_extent_npoints(file_space);
360 :
361 : // resize dataset
362 0 : hsize_t new_size[RANK] = {count + n};
363 0 : H5Dset_extent(dset, new_size);
364 :
365 : // get updated filespace
366 0 : H5Sclose(file_space);
367 0 : file_space = H5Dget_space(dset);
368 :
369 0 : hsize_t offset[RANK] = {count};
370 0 : hsize_t cnt[RANK] = {n};
371 :
372 0 : H5Sselect_hyperslab(file_space, H5S_SELECT_SET, offset, NULL, cnt, NULL);
373 0 : hid_t mspace_id = H5Screate_simple(RANK, cnt, NULL);
374 :
375 0 : H5Dwrite(dset, sid, mspace_id, file_space, H5P_DEFAULT, buffer.data());
376 :
377 0 : H5Sclose(mspace_id);
378 0 : H5Sclose(file_space);
379 :
380 : buffer.clear();
381 :
382 0 : H5Fflush(file, H5F_SCOPE_GLOBAL);
383 : }
384 :
385 0 : std::string HDF5Output::getDescription() const {
386 0 : return "HDF5Output";
387 : }
388 :
389 0 : void HDF5Output::setFlushLimit(unsigned int N)
390 : {
391 0 : flushLimit = N;
392 0 : }
393 :
394 : } // namespace crpropa
395 :
396 : #endif // CRPROPA_HAVE_HDF5
|
