Gateau User Manual
Atmospheric simulation of astronomical signals
fio.h
Go to the documentation of this file.
1 /*! \file fio.h
2  * \brief File input/output operations.
3  **/
4 
5 #include <iostream>
6 #include <iomanip>
7 #include <cmath>
8 #include <fstream>
9 #include <string>
10 #include <sstream>
11 #include <filesystem>
12 #include <vector>
13 
14 #include "hdf5.h"
15 
16 #include "structs.h"
17 
18 #ifndef __FILEIO_H
19 #define __FILEIO_H
20 
21 namespace fs = std::filesystem;
22 
23 #define NPWVATM 55
24 #define NFREQ 8301
25 #define NATMGRID 3
26 #define COLBUFF 16
27 #define FMTBUFF 4
28 
29 #define OBSATTRS_NAME "OBSATTRS"
30 #define SPAX_NAME "SPAXEL"
31 #define NSPAX_NAME "num_spax"
32 #define FREQS_NAME "frequencies"
33 #define ETA_AP_NAME "eta_ap"
34 #define TIME_NAME "times"
35 #define AZ_NAME "az"
36 #define EL_NAME "el"
37 #define PWV_NAME "pwv"
38 
39 #define AZ_SPAX_NAME "az_spax"
40 #define EL_SPAX_NAME "el_spax"
41 #define OUT_NAME "data"
42 
43 #define CHBUFF 100
44 #define RANK1D 1
45 #define RANK2D 2
46 
47 void readAtmMeta(int **meta, std::string path);
48 
49 template <typename T, typename U>
50 void readEtaATM(
51  T **eta_array,
52  U *pwv_atm,
53  U *freq_atm
54  );
55 
56 //template <typename T, typename U>
57 //void readAtmScreen(
58 // T **PWV_screen,
59 // U *x_spec,
60 // U *y_spec,
61 // std::string path,
62 // std::vector<std::string> datp
63 // );
64 template <typename T, typename U>
65 void readAtmScreen(
66  T **PWV_screen,
67  U *x_spec,
68  U *y_spec,
69  std::string path,
70  std::string datp
71  );
72 
73 class OutputFile
74 {
75  private:
76  const char *filename;
77  bool dry_run;
78 
79  // Handles for hdf5 file and groups
80  hid_t file_id;
81  hid_t obsattrs_id;
82  hid_t spax_id;
83 
84  // Handles for hdf5 dataspace/set creation
85  // Special one for PWV because its easier
86  hid_t dspace_id, dspace_pwv_id, dspace_slab_id, dspace_pwv_slab_id, dset_id, dset_pwv_id;
87 
88  // Output dimensions
89  hsize_t dims_1D[RANK1D],
90  dims_2D[RANK2D],
91  dims_2D_chunk[RANK2D],
92  dims_2D_stride[RANK2D];
93 
94  hsize_t dims_2D_null[2] = {0,0};
95 
96  // Hyperslab dimensions
97  hsize_t start[2], count[2], start_pink[2], count_pink[2], start_pwv[1], count_pwv[1];
98 
99  int ntimes, nfreqs;
100  int offset_times = 0;
101  int offset_freqs = 0;
102  int offset_times_pwv = 0;
103 
104  void check_API_call_status(
105  herr_t status,
106  int loc
107  )
108  {
109  if(status < 0) {printf("HDF5 API error occured on line %d\n", loc);}
110  }
111 
112  public:
113  OutputFile(
114  const char *filename,
115  bool dry_run,
116  int nspax,
117  int ntimes,
118  int nfreqs,
119  float *freqs,
120  float *times,
121  float *eta_ap,
122  float *az,
123  float *el
124  )
125  {
126  this->filename = filename;
127  this->dry_run = dry_run;
128  this->ntimes = ntimes;
129  this->nfreqs = nfreqs;
130  printf("%d\n", this->dry_run);
131 
132  if (this->dry_run) {return;}
133 
134  dims_2D[0] = nfreqs;
135  dims_2D[1] = ntimes;
136 
137  start[0] = 0;
138  count[0] = nfreqs;
139 
140  start_pink[1] = 0;
141  count_pink[1] = ntimes;
142 
143  // Make file and obsattrs group
144  file_id = H5Fcreate(
145  filename,
146  H5F_ACC_TRUNC,
147  H5P_DEFAULT,
148  H5P_DEFAULT
149  );
150 
151  obsattrs_id = H5Gcreate(
152  file_id,
153  OBSATTRS_NAME,
154  H5P_DEFAULT,
155  H5P_DEFAULT,
156  H5P_DEFAULT
157  );
158 
159  // Write dataspace with number of spaxels
160  dspace_id = H5Screate(H5S_SCALAR);
161 
162  dset_id = H5Dcreate(
163  obsattrs_id,
164  NSPAX_NAME,
165  H5T_NATIVE_INT,
166  dspace_id,
167  H5P_DEFAULT,
168  H5P_DEFAULT,
169  H5P_DEFAULT
170  );
171 
172  check_API_call_status(
173  H5Dwrite(
174  dset_id,
175  H5T_NATIVE_INT,
176  H5S_ALL,
177  H5S_ALL,
178  H5P_DEFAULT,
179  &nspax
180  ),
181  __LINE__
182  );
183 
184  check_API_call_status(H5Dclose(dset_id), __LINE__);
185 
186  // Write dataspace for frequencies
187  dims_1D[0] = nfreqs;
188 
189  dspace_id = H5Screate_simple(
190  RANK1D,
191  dims_1D,
192  NULL
193  );
194 
195  dset_id = H5Dcreate(
196  obsattrs_id,
197  FREQS_NAME,
198  H5T_IEEE_F32LE,
199  dspace_id,
200  H5P_DEFAULT,
201  H5P_DEFAULT,
202  H5P_DEFAULT
203  );
204 
205  check_API_call_status(
206  H5Dwrite(
207  dset_id,
208  H5T_IEEE_F32LE,
209  H5S_ALL,
210  H5S_ALL,
211  H5P_DEFAULT,
212  freqs
213  ),
214  __LINE__
215  );
216  check_API_call_status(H5Dclose(dset_id), __LINE__);
217 
218  // Write eta_ap
219  dset_id = H5Dcreate(
220  obsattrs_id,
221  ETA_AP_NAME,
222  H5T_IEEE_F32LE,
223  dspace_id,
224  H5P_DEFAULT,
225  H5P_DEFAULT,
226  H5P_DEFAULT
227  );
228 
229  check_API_call_status(
230  H5Dwrite(
231  dset_id,
232  H5T_IEEE_F32LE,
233  H5S_ALL,
234  H5S_ALL,
235  H5P_DEFAULT,
236  eta_ap
237  ),
238  __LINE__
239  );
240  check_API_call_status(H5Dclose(dset_id), __LINE__);
241  check_API_call_status(H5Sclose(dspace_id), __LINE__);
242 
243  // Write time and az-el arrays
244  dims_1D[0] = ntimes;
245 
246  dspace_id = H5Screate_simple(
247  RANK1D,
248  dims_1D,
249  NULL
250  );
251 
252  dset_id = H5Dcreate(
253  obsattrs_id,
254  TIME_NAME,
255  H5T_IEEE_F32LE,
256  dspace_id,
257  H5P_DEFAULT,
258  H5P_DEFAULT,
259  H5P_DEFAULT
260  );
261 
262  check_API_call_status(
263  H5Dwrite(
264  dset_id,
265  H5T_IEEE_F32LE,
266  H5S_ALL,
267  H5S_ALL,
268  H5P_DEFAULT,
269  times
270  ),
271  __LINE__
272  );
273  check_API_call_status(H5Dclose(dset_id), __LINE__);
274 
275  dset_id = H5Dcreate(
276  obsattrs_id,
277  AZ_NAME,
278  H5T_IEEE_F32LE,
279  dspace_id,
280  H5P_DEFAULT,
281  H5P_DEFAULT,
282  H5P_DEFAULT
283  );
284 
285  check_API_call_status(
286  H5Dwrite(
287  dset_id,
288  H5T_IEEE_F32LE,
289  H5S_ALL,
290  H5S_ALL,
291  H5P_DEFAULT,
292  az
293  ),
294  __LINE__
295  );
296  check_API_call_status(H5Dclose(dset_id), __LINE__);
297 
298  dset_id = H5Dcreate(
299  obsattrs_id,
300  EL_NAME,
301  H5T_IEEE_F32LE,
302  dspace_id,
303  H5P_DEFAULT,
304  H5P_DEFAULT,
305  H5P_DEFAULT
306  );
307 
308  check_API_call_status(
309  H5Dwrite(
310  dset_id,
311  H5T_IEEE_F32LE,
312  H5S_ALL,
313  H5S_ALL,
314  H5P_DEFAULT,
315  el
316  ),
317  __LINE__
318  );
319  check_API_call_status(H5Dclose(dset_id), __LINE__);
320  check_API_call_status(H5Sclose(dspace_id), __LINE__);
321 
322  // Allocate output array for zenith PWV
323  dspace_pwv_id = H5Screate_simple(
324  RANK1D,
325  dims_1D,
326  NULL
327  );
328 
329  dset_pwv_id = H5Dcreate(
330  obsattrs_id,
331  PWV_NAME,
332  H5T_IEEE_F32LE,
333  dspace_pwv_id,
334  H5P_DEFAULT,
335  H5P_DEFAULT,
336  H5P_DEFAULT
337  );
338  }
339  void write_chunk_to_pwv(
340  int ntimes_chunk,
341  float *data
342  )
343  {
344  if (this->dry_run) {return;}
345  start_pwv[0] = offset_times_pwv;
346  count_pwv[0] = ntimes_chunk;
347 
348  dims_1D[0] = ntimes_chunk;
349 
350  check_API_call_status(
351  H5Sselect_hyperslab(
352  dspace_pwv_id,
353  H5S_SELECT_SET,
354  start_pwv,
355  NULL,
356  count_pwv,
357  NULL
358  ),
359  __LINE__
360  );
361 
362  dspace_pwv_slab_id = H5Screate_simple(
363  RANK1D,
364  dims_1D,
365  NULL
366  );
367 
368  check_API_call_status(
369  H5Dwrite(
370  dset_pwv_id,
371  H5T_IEEE_F32LE,
372  dspace_pwv_slab_id,
373  dspace_pwv_id,
374  H5P_DEFAULT,
375  data
376  ),
377  __LINE__
378  );
379 
380  offset_times_pwv += ntimes_chunk;
381  check_API_call_status(H5Sclose(dspace_pwv_slab_id), __LINE__);
382  }
383  void close_obsattrs() {
384  if (this->dry_run) {return;}
385  check_API_call_status(H5Dclose(dset_pwv_id), __LINE__);
386  check_API_call_status(H5Sclose(dspace_pwv_id), __LINE__);
387  check_API_call_status(H5Gclose(obsattrs_id), __LINE__);
388  }
389 
390  void open_spaxel(
391  int spax_index,
392  float az_spax,
393  float el_spax
394  )
395  {
396  if (this->dry_run) {return;}
397  offset_times = 0;
398  offset_freqs = 0;
399 
400  char spax_name[CHBUFF] = SPAX_NAME;
401  char buffer[CHBUFF];
402 
403  sprintf(buffer, "%d", spax_index);
404  strcat(spax_name, buffer);
405 
406  spax_id = H5Gcreate(
407  file_id,
408  spax_name,
409  H5P_DEFAULT,
410  H5P_DEFAULT,
411  H5P_DEFAULT
412  );
413 
414  dspace_id = H5Screate(H5S_SCALAR);
415 
416  dset_id = H5Dcreate(
417  spax_id,
418  AZ_SPAX_NAME,
419  H5T_IEEE_F32LE,
420  dspace_id,
421  H5P_DEFAULT,
422  H5P_DEFAULT,
423  H5P_DEFAULT
424  );
425 
426  check_API_call_status(
427  H5Dwrite(
428  dset_id,
429  H5T_IEEE_F32LE,
430  H5S_ALL,
431  H5S_ALL,
432  H5P_DEFAULT,
433  &az_spax
434  ),
435  __LINE__
436  );
437  check_API_call_status(H5Dclose(dset_id), __LINE__);
438 
439  dset_id = H5Dcreate(
440  spax_id,
441  EL_SPAX_NAME,
442  H5T_IEEE_F32LE,
443  dspace_id,
444  H5P_DEFAULT,
445  H5P_DEFAULT,
446  H5P_DEFAULT
447  );
448 
449  check_API_call_status(
450  H5Dwrite(
451  dset_id,
452  H5T_IEEE_F32LE,
453  H5S_ALL,
454  H5S_ALL,
455  H5P_DEFAULT,
456  &el_spax
457  ),
458  __LINE__
459  );
460  check_API_call_status(H5Dclose(dset_id), __LINE__);
461  check_API_call_status(H5Sclose(dspace_id), __LINE__);
462 
463  // Allocate output array for this spaxel
464  dspace_id = H5Screate_simple(
465  RANK2D,
466  dims_2D,
467  NULL
468  );
469 
470  dset_id = H5Dcreate(
471  spax_id,
472  OUT_NAME,
473  H5T_IEEE_F32LE,
474  dspace_id,
475  H5P_DEFAULT,
476  H5P_DEFAULT,
477  H5P_DEFAULT
478  );
479  }
480 
481  void write_pink_chunk_to_spaxel(
482  int k_ch,
483  float *data
484  )
485  {
486  if (this->dry_run) {return;}
487  start_pink[0] = k_ch;
488  count_pink[0] = 1;
489 
490  dims_1D[0] = ntimes;
491 
492  check_API_call_status(
493  H5Sselect_hyperslab(
494  dspace_id,
495  H5S_SELECT_SET,
496  start_pink,
497  NULL,
498  count_pink,
499  NULL
500  ),
501  __LINE__
502  );
503 
504  dspace_slab_id = H5Screate_simple(
505  RANK1D,
506  dims_1D,
507  NULL
508  );
509 
510  check_API_call_status(
511  H5Dwrite(
512  dset_id,
513  H5T_IEEE_F32LE,
514  dspace_slab_id,
515  dspace_id,
516  H5P_DEFAULT,
517  data
518  ),
519  __LINE__
520  );
521 
522  check_API_call_status(H5Sclose(dspace_slab_id), __LINE__);
523  }
524 
525  void write_chunk_to_spaxel(
526  int ntimes_chunk,
527  float *data
528  )
529  {
530  if (this->dry_run) {return;}
531  start[1] = offset_times;
532  count[1] = ntimes_chunk;
533 
534  dims_1D[0] = ntimes_chunk * nfreqs;
535 
536  check_API_call_status(
537  H5Sselect_hyperslab(
538  dspace_id,
539  H5S_SELECT_SET,
540  start,
541  NULL,
542  count,
543  NULL
544  ),
545  __LINE__
546  );
547 
548  dspace_slab_id = H5Screate_simple(
549  RANK1D,
550  dims_1D,
551  NULL
552  );
553 
554  float *buffer = new float[dims_1D[0]];
555  check_API_call_status(
556  H5Dread(
557  dset_id,
558  H5T_IEEE_F32LE,
559  dspace_slab_id,
560  dspace_id,
561  H5P_DEFAULT,
562  buffer
563  ),
564  __LINE__
565  );
566 
567  for(int ii=0; ii < dims_1D[0]; ii++)
568  {
569  buffer[ii] += data[ii];
570  }
571 
572  check_API_call_status(
573  H5Dwrite(
574  dset_id,
575  H5T_IEEE_F32LE,
576  dspace_slab_id,
577  dspace_id,
578  H5P_DEFAULT,
579  buffer
580  ),
581  __LINE__
582  );
583 
584  delete[] buffer;
585 
586  offset_times += ntimes_chunk;
587  check_API_call_status(H5Sclose(dspace_slab_id), __LINE__);
588  }
589 
590  void close_spaxel()
591  {
592  if (this->dry_run) {return;}
593  check_API_call_status(H5Dclose(dset_id), __LINE__);
594  check_API_call_status(H5Sclose(dspace_id), __LINE__);
595  check_API_call_status(H5Gclose(spax_id), __LINE__);
596  }
597  ~OutputFile()
598  {
599  if (this->dry_run) {return;}
600  check_API_call_status(H5Fclose(file_id), __LINE__);
601  }
602 };
603 #endif
604 
605 void readAtmMeta(
606  int **meta,
607  std::string path
608  )
609 {
610  fs::path dir(path);
611  fs::path file("atm_meta.datp");
612  fs::path abs_loc = dir / file;
613 
614  *meta = new int[NATMGRID];
615 
616  std::string store;
617 
618  std::ifstream myfile(abs_loc);
619  std::string line;
620 
621  int idx = 0;
622 
623  if (!myfile)
624  {
625  std::cerr
626  << "Could not open the file at "
627  << abs_loc
628  << std::endl;
629  exit(5);
630  }
631  else
632  {
633  while(std::getline(myfile, line))
634  {
635  std::istringstream iss(line);
636  while(std::getline(iss, store, ' '))
637  {
638  if (store=="") {continue;}
639  (*meta)[idx] = std::stoi(store);
640  idx++;
641  }
642  }
643  myfile.close();
644  }
645 }
646 
647 template <typename T, typename U>
648 void readEtaATM(
649  T **eta_array,
650  U *pwv_atm,
651  U *freq_atm,
652  const char* filepath
653  )
654 {
655  // TODO read these in from file? Ask Arend
656  pwv_atm->start = 0.1;
657  pwv_atm->step = 0.1;
658  pwv_atm->num = NPWVATM;
659 
660  freq_atm->start = 70.e9;
661  freq_atm->step = 0.1e9;
662  freq_atm->num = NFREQ;
663 
664  *eta_array = new T[NPWVATM * NFREQ];
665 
666  std::string store;
667  //std::cout << abi::__cxa_demangle(typeid(store).name(), NULL, NULL, &status) << std::endl;
668 
669  std::ifstream myfile(filepath);
670  std::string line;
671 
672  int line_nr = 0;
673  int idx = 0;
674 
675  if (!myfile)
676  {
677  std::cerr
678  << "Could not open the resource file at "
679  << filepath
680  << std::endl;
681  exit(5);
682  /* TODO Standardize exit codes */
683  }
684 
685  while(std::getline(myfile, line))
686  {
687  std::istringstream iss(line);
688  if(!line_nr)
689  {
690  line_nr++;
691  continue;
692  }
693 
694  while(std::getline(iss, store, ' '))
695  {
696  if(!idx)
697  {
698  idx++;
699  continue;
700  }
701  else if (store=="")
702  {
703  continue;
704  }
705 
706  (*eta_array)[NFREQ * (idx-1) + (line_nr - 1)] = std::stof(store);
707  idx++;
708  }
709  line_nr++;
710  idx = 0;
711  }
712  myfile.close();
713 }
714 
715  template <typename T, typename U>
716 void readAtmScreen(
717  T **PWV_screen,
718  U *x_spec,
719  U *y_spec,
720  std::string path,
721  std::string datp
722  )
723 {
724  *PWV_screen = new T[x_spec->num * y_spec->num];
725 
726  std::string store;
727  std::string line;
728 
729  fs::path dir(path);
730 
731  fs::path file(datp);
732  fs::path abs_loc = dir / file;
733 
734  std::ifstream myfile(abs_loc);
735  int line_nr = 0;
736  int idx = 0;
737 
738  if (!myfile)
739  {
740  std::cerr
741  << "Could not open the file!"
742  << std::endl;
743  }
744  else
745  {
746  while(std::getline(myfile, line))
747  {
748  std::istringstream iss(line);
749  while(std::getline(iss, store, ' '))
750  {
751  if (store=="") {continue;}
752  (*PWV_screen)[y_spec->num * line_nr + idx] = std::stof(store);
753  idx++;
754  }
755  line_nr++;
756  idx = 0;
757  }
758  myfile.close();
759  }
760 }
761 
762 //template <typename T, typename U>
763 //void readAtmScreen(
764 // T **PWV_screen,
765 // U *x_spec,
766 // U *y_spec,
767 // std::string path,
768 // std::vector<std::string> datp
769 // )
770 //{
771 // int line_nr = 0;
772 // int idx = 0;
773 // int num_screen = datp.size();
774 //
775 // *PWV_screen = new T[num_screen * x_spec->num * y_spec->num];
776 //
777 // std::string store;
778 // std::string line;
779 //
780 // fs::path dir(path);
781 //
782 // for(int i=0; i<num_screen; i++)
783 // {
784 // fs::path file(datp[i]);
785 // fs::path abs_loc = dir / file;
786 //
787 // std::ifstream myfile(abs_loc);
788 //
789 // if (!myfile)
790 // {
791 // std::cerr
792 // << "Could not open the file!"
793 // << std::endl;
794 // }
795 // else
796 // {
797 // while(std::getline(myfile, line))
798 // {
799 // std::istringstream iss(line);
800 // while(std::getline(iss, store, ' '))
801 // {
802 // if (store=="") {continue;}
803 // (*PWV_screen)[y_spec->num * line_nr + idx] = std::stof(store);
804 // idx++;
805 // }
806 // line_nr++;
807 // idx = 0;
808 // }
809 // myfile.close();
810 // }
811 // }
812 //}
structs.h
Data structures for receiving data from Python interface.
OutputFile
Definition: fio.h:73