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 2 : }
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(RedshiftColumn))
110 0 : H5Tinsert(sid, "z", HOFFSET(OutputRow, z), H5T_NATIVE_DOUBLE);
111 0 : if (fields.test(SerialNumberColumn))
112 0 : H5Tinsert(sid, "SN", HOFFSET(OutputRow, SN), H5T_NATIVE_UINT64);
113 0 : if (fields.test(CurrentIdColumn))
114 0 : H5Tinsert(sid, "ID", HOFFSET(OutputRow, ID), H5T_NATIVE_INT32);
115 0 : if (fields.test(CurrentEnergyColumn))
116 0 : H5Tinsert(sid, "E", HOFFSET(OutputRow, E), H5T_NATIVE_DOUBLE);
117 0 : if (fields.test(CurrentPositionColumn) && oneDimensional)
118 0 : H5Tinsert(sid, "X", HOFFSET(OutputRow, X), H5T_NATIVE_DOUBLE);
119 0 : if (fields.test(CurrentPositionColumn) && not oneDimensional) {
120 0 : H5Tinsert(sid, "X", HOFFSET(OutputRow, X), H5T_NATIVE_DOUBLE);
121 0 : H5Tinsert(sid, "Y", HOFFSET(OutputRow, Y), H5T_NATIVE_DOUBLE);
122 0 : H5Tinsert(sid, "Z", HOFFSET(OutputRow, Z), H5T_NATIVE_DOUBLE);
123 : }
124 0 : if (fields.test(CurrentDirectionColumn) && not oneDimensional) {
125 0 : H5Tinsert(sid, "Px", HOFFSET(OutputRow, Px), H5T_NATIVE_DOUBLE);
126 0 : H5Tinsert(sid, "Py", HOFFSET(OutputRow, Py), H5T_NATIVE_DOUBLE);
127 0 : H5Tinsert(sid, "Pz", HOFFSET(OutputRow, Pz), H5T_NATIVE_DOUBLE);
128 : }
129 0 : if (fields.test(SerialNumberColumn))
130 0 : H5Tinsert(sid, "SN0", HOFFSET(OutputRow, SN0), H5T_NATIVE_UINT64);
131 0 : if (fields.test(SourceIdColumn))
132 0 : H5Tinsert(sid, "ID0", HOFFSET(OutputRow, ID0), H5T_NATIVE_INT32);
133 0 : if (fields.test(SourceEnergyColumn))
134 0 : H5Tinsert(sid, "E0", HOFFSET(OutputRow, E0), H5T_NATIVE_DOUBLE);
135 0 : if (fields.test(SourcePositionColumn) && oneDimensional)
136 0 : H5Tinsert(sid, "X0", HOFFSET(OutputRow, X0), H5T_NATIVE_DOUBLE);
137 0 : if (fields.test(SourcePositionColumn) && not oneDimensional){
138 0 : H5Tinsert(sid, "X0", HOFFSET(OutputRow, X0), H5T_NATIVE_DOUBLE);
139 0 : H5Tinsert(sid, "Y0", HOFFSET(OutputRow, Y0), H5T_NATIVE_DOUBLE);
140 0 : H5Tinsert(sid, "Z0", HOFFSET(OutputRow, Z0), H5T_NATIVE_DOUBLE);
141 : }
142 0 : if (fields.test(SourceDirectionColumn) && not oneDimensional) {
143 0 : H5Tinsert(sid, "P0x", HOFFSET(OutputRow, P0x), H5T_NATIVE_DOUBLE);
144 0 : H5Tinsert(sid, "P0y", HOFFSET(OutputRow, P0y), H5T_NATIVE_DOUBLE);
145 0 : H5Tinsert(sid, "P0z", HOFFSET(OutputRow, P0z), H5T_NATIVE_DOUBLE);
146 : }
147 0 : if (fields.test(SerialNumberColumn))
148 0 : H5Tinsert(sid, "SN1", HOFFSET(OutputRow, SN1), H5T_NATIVE_UINT64);
149 0 : if (fields.test(CreatedIdColumn))
150 0 : H5Tinsert(sid, "ID1", HOFFSET(OutputRow, ID1), H5T_NATIVE_INT32);
151 0 : if (fields.test(CreatedEnergyColumn))
152 0 : H5Tinsert(sid, "E1", HOFFSET(OutputRow, E1), H5T_NATIVE_DOUBLE);
153 0 : if (fields.test(CreatedPositionColumn) && oneDimensional)
154 0 : H5Tinsert(sid, "X1", HOFFSET(OutputRow, X1), H5T_NATIVE_DOUBLE);
155 0 : if (fields.test(CreatedPositionColumn) && not oneDimensional) {
156 0 : H5Tinsert(sid, "X1", HOFFSET(OutputRow, X1), H5T_NATIVE_DOUBLE);
157 0 : H5Tinsert(sid, "Y1", HOFFSET(OutputRow, Y1), H5T_NATIVE_DOUBLE);
158 0 : H5Tinsert(sid, "Z1", HOFFSET(OutputRow, Z1), H5T_NATIVE_DOUBLE);
159 : }
160 0 : if (fields.test(CreatedDirectionColumn) && not oneDimensional) {
161 0 : H5Tinsert(sid, "P1x", HOFFSET(OutputRow, P1x), H5T_NATIVE_DOUBLE);
162 0 : H5Tinsert(sid, "P1y", HOFFSET(OutputRow, P1y), H5T_NATIVE_DOUBLE);
163 0 : H5Tinsert(sid, "P1z", HOFFSET(OutputRow, P1z), H5T_NATIVE_DOUBLE);
164 : }
165 0 : if (fields.test(WeightColumn))
166 0 : H5Tinsert(sid, "W", HOFFSET(OutputRow, weight), H5T_NATIVE_DOUBLE);
167 :
168 0 : if (fields.test(CandidateTagColumn))
169 0 : H5Tinsert(sid, "tag", HOFFSET(OutputRow, tag), H5T_C_S1);
170 :
171 0 : size_t pos = 0;
172 0 : for(std::vector<Output::Property>::const_iterator iter = properties.begin();
173 0 : iter != properties.end(); ++iter)
174 : {
175 0 : hid_t type = variantTypeToH5T_NATIVE((*iter).defaultValue.getType());
176 0 : if (type == H5T_C_S1)
177 : { // set size of string field to size of default value!
178 0 : type = H5Tcopy(H5T_C_S1);
179 0 : H5Tset_size(type, (*iter).defaultValue.toString().size());
180 : }
181 :
182 0 : H5Tinsert(sid, (*iter).name.c_str(), HOFFSET(OutputRow, propertyBuffer) + pos, type);
183 0 : pos += (*iter).defaultValue.getSize();
184 : }
185 0 : if (pos >= propertyBufferSize)
186 : {
187 0 : KISS_LOG_ERROR << "Using " << pos << " bytes for properties output. Maximum is " << propertyBufferSize << " bytes.";
188 0 : throw std::runtime_error("Size of property buffer exceeded");
189 : }
190 :
191 : // chunked prop
192 0 : hid_t plist = H5Pcreate(H5P_DATASET_CREATE);
193 0 : H5Pset_layout(plist, H5D_CHUNKED);
194 0 : hsize_t chunk_dims[RANK] = {BUFFER_SIZE};
195 0 : H5Pset_chunk(plist, RANK, chunk_dims);
196 0 : H5Pset_deflate(plist, 5);
197 :
198 0 : hsize_t dims[RANK] = {0};
199 0 : hsize_t max_dims[RANK] = {H5S_UNLIMITED};
200 0 : dataspace = H5Screate_simple(RANK, dims, max_dims);
201 :
202 0 : dset = H5Dcreate2(file, "CRPROPA3", sid, dataspace, H5P_DEFAULT, plist, H5P_DEFAULT);
203 :
204 0 : insertStringAttribute("OutputType", outputName);
205 0 : insertStringAttribute("Version", g_GIT_DESC);
206 0 : insertDoubleAttribute("LengthScale", this->lengthScale);
207 0 : insertDoubleAttribute("EnergyScale", this->energyScale);
208 :
209 : // add ranom seeds
210 0 : std::vector< std::vector<uint32_t> > seeds = Random::getSeedThreads();
211 0 : for (size_t i = 0; i < seeds.size(); i++)
212 : {
213 : hid_t type, attr_space, version_attr;
214 : herr_t status;
215 0 : hsize_t dims[] = {1, 0};
216 0 : dims[1] = seeds[i].size();
217 :
218 0 : type = H5Tarray_create(H5T_NATIVE_ULONG, 2, dims);
219 :
220 0 : attr_space = H5Screate_simple(0, dims, NULL);
221 : char nameBuffer[256];
222 : sprintf(nameBuffer, "SEED_%03lu", i);
223 0 : KISS_LOG_DEBUG << "Creating HDF5 attribute: " << nameBuffer << " with dimensions " << dims[0] << "x" << dims[1] ;
224 :
225 0 : version_attr = H5Acreate2(dset, nameBuffer, type, attr_space, H5P_DEFAULT, H5P_DEFAULT);
226 0 : status = H5Awrite(version_attr, type, &seeds[i][0]);
227 0 : status = H5Aclose(version_attr);
228 0 : status = H5Sclose(attr_space);
229 :
230 : }
231 :
232 :
233 0 : H5Pclose(plist);
234 :
235 0 : buffer.reserve(BUFFER_SIZE);
236 0 : time(&lastFlush);
237 0 : }
238 :
239 1 : void HDF5Output::close() {
240 1 : if (file >= 0) {
241 0 : flush();
242 0 : H5Dclose(dset);
243 0 : H5Tclose(sid);
244 0 : H5Sclose(dataspace);
245 0 : H5Fclose(file);
246 0 : file = -1;
247 : }
248 1 : }
249 :
250 0 : void HDF5Output::process(Candidate* candidate) const {
251 0 : #pragma omp critical
252 : {
253 0 : if (file == -1)
254 : // This is ugly, but necesary as otherwise the user has to manually open the
255 : // file before processing the first candidate
256 0 : const_cast<HDF5Output*>(this)->open(filename);
257 : }
258 :
259 : OutputRow r;
260 0 : r.D = candidate->getTrajectoryLength() / lengthScale;
261 0 : r.z = candidate->getRedshift();
262 :
263 0 : r.SN = candidate->getSerialNumber();
264 0 : r.ID = candidate->current.getId();
265 0 : r.E = candidate->current.getEnergy() / energyScale;
266 0 : Vector3d v = candidate->current.getPosition() / lengthScale;
267 0 : r.X = v.x;
268 0 : r.Y = v.y;
269 0 : r.Z = v.z;
270 0 : v = candidate->current.getDirection();
271 0 : r.Px = v.x;
272 0 : r.Py = v.y;
273 0 : r.Pz = v.z;
274 :
275 0 : r.SN0 = candidate->getSourceSerialNumber();
276 0 : r.ID0 = candidate->source.getId();
277 0 : r.E0 = candidate->source.getEnergy() / energyScale;
278 0 : v = candidate->source.getPosition() / lengthScale;
279 0 : r.X0 = v.x;
280 0 : r.Y0 = v.y;
281 0 : r.Z0 = v.z;
282 0 : v = candidate->source.getDirection();
283 0 : r.P0x = v.x;
284 0 : r.P0y = v.y;
285 0 : r.P0z = v.z;
286 :
287 0 : r.SN1 = candidate->getCreatedSerialNumber();
288 0 : r.ID1 = candidate->created.getId();
289 0 : r.E1 = candidate->created.getEnergy() / energyScale;
290 0 : v = candidate->created.getPosition() / lengthScale;
291 0 : r.X1 = v.x;
292 0 : r.Y1 = v.y;
293 0 : r.Z1 = v.z;
294 0 : v = candidate->created.getDirection();
295 0 : r.P1x = v.x;
296 0 : r.P1y = v.y;
297 0 : r.P1z = v.z;
298 :
299 0 : r.weight= candidate->getWeight();
300 :
301 0 : r.tag = candidate->getTagOrigin();
302 :
303 : size_t pos = 0;
304 0 : for(std::vector<Output::Property>::const_iterator iter = properties.begin();
305 0 : iter != properties.end(); ++iter)
306 : {
307 0 : Variant v;
308 0 : if (candidate->hasProperty((*iter).name))
309 : {
310 0 : v = candidate->getProperty((*iter).name);
311 : }
312 : else
313 : {
314 0 : v = (*iter).defaultValue;
315 : }
316 0 : pos += v.copyToBuffer(&r.propertyBuffer[pos]);
317 0 : }
318 :
319 0 : #pragma omp critical
320 : {
321 0 : const_cast<HDF5Output*>(this)->candidatesSinceFlush++;
322 0 : Output::process(candidate);
323 :
324 0 : buffer.push_back(r);
325 :
326 :
327 0 : if (buffer.size() >= buffer.capacity())
328 : {
329 0 : KISS_LOG_DEBUG << "HDF5Output: Flush due to buffer capacity exceeded";
330 0 : flush();
331 : }
332 0 : else if (candidatesSinceFlush >= flushLimit)
333 : {
334 0 : KISS_LOG_DEBUG << "HDF5Output: Flush due to number of candidates";
335 0 : flush();
336 : }
337 0 : else if (difftime(time(NULL), lastFlush) > 60*10)
338 : {
339 0 : KISS_LOG_DEBUG << "HDF5Output: Flush due to time exceeded";
340 0 : flush();
341 : }
342 : }
343 0 : }
344 :
345 0 : void HDF5Output::flush() const {
346 0 : const_cast<HDF5Output*>(this)->lastFlush = time(NULL);
347 0 : const_cast<HDF5Output*>(this)->candidatesSinceFlush = 0;
348 :
349 : hsize_t n = buffer.size();
350 :
351 0 : if (n == 0)
352 0 : return;
353 :
354 0 : hid_t file_space = H5Dget_space(dset);
355 0 : hsize_t count = H5Sget_simple_extent_npoints(file_space);
356 :
357 : // resize dataset
358 0 : hsize_t new_size[RANK] = {count + n};
359 0 : H5Dset_extent(dset, new_size);
360 :
361 : // get updated filespace
362 0 : H5Sclose(file_space);
363 0 : file_space = H5Dget_space(dset);
364 :
365 0 : hsize_t offset[RANK] = {count};
366 0 : hsize_t cnt[RANK] = {n};
367 :
368 0 : H5Sselect_hyperslab(file_space, H5S_SELECT_SET, offset, NULL, cnt, NULL);
369 0 : hid_t mspace_id = H5Screate_simple(RANK, cnt, NULL);
370 :
371 0 : H5Dwrite(dset, sid, mspace_id, file_space, H5P_DEFAULT, buffer.data());
372 :
373 0 : H5Sclose(mspace_id);
374 0 : H5Sclose(file_space);
375 :
376 0 : buffer.clear();
377 :
378 0 : H5Fflush(file, H5F_SCOPE_GLOBAL);
379 : }
380 :
381 0 : std::string HDF5Output::getDescription() const {
382 0 : return "HDF5Output";
383 : }
384 :
385 0 : void HDF5Output::setFlushLimit(unsigned int N)
386 : {
387 0 : flushLimit = N;
388 0 : }
389 :
390 : } // namespace crpropa
391 :
392 : #endif // CRPROPA_HAVE_HDF5
|
