Hi Ben,

I have been working on this and I have two questions. Let me explain first. I decided to do it in the following way.

In lightcone_construction.h I changed

Konsole output size_tNlines= 0,Nblocks=0; while(!feof(pFile)){ // loop through blocks

//long Nbatch = unit.scan_block(blocksize,pFile); longNbatch= unit.scan_blockH5(blocksize,snap_filenames[i_file]);

so that Nbatch is changed to the new version I included to read hdf5. I also added to the struct ASCII_XMRT12 a new way of reading:

Konsole output structASCII_XMRRT12:public ASCII_XMRRT{ size_tscan_block(size_tblocksize,FILE*pFile){

// read in a block of points size_ti=0; inttmp;

  points.resize(blocksize);

while(i < blocksize && fscanf(pFile,"%lf %lf %lf %lf %lf %lf %i" ,&points[i].x[0],&points[i].x[1],&points[i].x[2] ,&points[i].mass,&points[i].r_max,&points[i].r_scale,&tmp) != EOF){ if(tmp == 1 || tmp == 2) ++i; } points.resize(i);

for(auto&h: points){ h.r_max /= 1.0e3; h.r_scale /= 1.0e3; } returni; }

//CHANGED ***** MARCOS size_tscan_blockH5(size_tblocksize,H5std_stringfilename){ // read in a block of points size_ti=0, j=0; intk=6; hsize_tnrows, ncols=7, offset=0; doubledata; data = common_scan_blockH5(filename, "/dset", blocksize, ncols, offset, &nrows); points.resize(nrows); for(i = 0; i < nrows; i++) { if(((size_t)data[incols+k] == 1) || ((size_t)data[incols+k] == 2)) { points[j].x[0] = data[incols ]; points[j].x[1] = data[incols+1]; points[j].x[2] = data[incols+2]; points[j].mass = data[incols+3]; points[j].r_max = data[incols+4]; points[j].r_scale = data[i*ncols+5]; j++; } } points.resize(j); deletedata;

for(auto&h: points){ h.r_max /= 1.0e3; h.r_scale /= 1.0e3; } returnj; }

Konsole output doublecommon_scan_blockH5(H5std_stringfilename , H5std_stringdataset_name , hsize_tn_rows , hsize_tn_cols , hsize_toffset_rows , hsize_tread_rows );

}; //CHANGED ***** MARCOS

Which calls a function called common_scan_blockH5() that I included at the end of MultiPlane/lightcone_construction.cpp. I am able to compile glamer but, when I try to compile cones.cpp I get some error that it's not finding the function common_scan_blockH5().

My questions are the following:

1) Where should I put the function common_scan_blockH5() so that lightcone_construction.h sees it. I also need to add some header ( and things like that) to be able to compile it, which place do you think is the best one?

2) I need to include the following flags

Konsole output -lhdf5_cpp -lhdf5

so that the compiler sees the libraries I need to use hdf5. I have been unable to find the place to include them in the Makefile, where do you put them?

I could pull the code to github but as it is not working yet, I thought it was better not to.

Cheers! Marcos.

PD. I attach the files I changed.

On 07/03/2017 09:08 PM, Ben Metcalf wrote:

@mpellejero https://github.com/mpellejero lets continue our discussion of this here.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/glenco/SLsimLib/issues/131, or mute the thread https://github.com/notifications/unsubscribe-auth/AVsIK2TKyXemjYSONUj9A49uJY4TAkLoks5sKUomgaJpZM4OMpVb.

ADVERTENCIA: Sobre la privacidad y cumplimiento de la Ley de Protección de Datos, acceda a http://www.iac.es/disclaimer.php WARNING: For more information on privacy and fulfilment of the Law concerning the Protection of Data, consult http://www.iac.es/disclaimer.php?lang=en

// // lightcone_construction.h // GLAMER // // Created by Ben Metcalf on 26/10/16. // //

ifndef GLAMERlightcone_construction__

define GLAMERlightcone_construction__

include

include "point.h"

include "lens_halos.h"

include "geometry.h"

ifdef OLD_HEADER_FILENAME

include

else

include

endif

using std::cout; using std::endl;

include

include "H5Cpp.h"

using namespace H5;

/* \brief The LightCones namespace is for classes and functions related to making light cones and weak lensing maps from 3D snapshots of particles or halos. / namespace LightCones{

struct DataRockStar{ unsigned int id; Point_3d x; double mass; double Rvir; // in comoving units double Rscale; // in comoving units double Vmax; };

struct DatumXM{ Point_3d x; double mass; }; struct DatumXMRmRs{ Point_3d x; double mass; double r_max; double r_scale; }; struct DatumXMR{ Point_3d x; double mass; double r; };

void ReadLightConeNFW(std::string filename,COSMOLOGY &cosmo ,std::vector<LensHalo* > &lensVec ,PosType &theta_max);

void ReadLightConeParticles(std::string filename,COSMOLOGY &cosmo ,std::vector<LensHaloParticles *> &lensVec ,int Nplanes ,float particle_mass ,float particle_size ,bool angular_sizes = false ,bool verbose = false);

void ReadLightConeParticles(std::string filename,COSMOLOGY &cosmo ,std::vector<LensHaloMassMap *> &lensVec ,int Nplanes ,float particle_mass ,float angular_resolution ,bool verbose = false);

/** \brief Class for constructing light cones form simulation boxes. / class LightCone{ public: LightCone(double angular_radius);

void WriteLightCone(std::string filename,std::vector<DataRockStar> &vec);
void WriteLightCone(std::string filename,std::vector<Point_3d> &vec);

friend class MultiLightCone;

template <typename T>
void select(Point_3d xo,Point_3d v,double Length,double rlow,double rhigh
            ,T* begin
            ,T* end
            ,Utilities::LockableContainer<std::vector<T> > &incone
            ,bool periodic_boundaries = true);

private:

double r_theta; // angular radius of light cone (radians)
double sin_theta_sqrt;

}; /// select the halos from the box that are within the light cone template void LightCone::select(Point_3d xo,Point_3d v,double Length,double rlow,double rhigh ,T begin ,T end ,Utilities::LockableContainer<std::vector > &incone ,bool periodic_boundaries){

if(begin == end) return;

Point_3d dx,x;
double r2,xp;
Point_3d n;
double rhigh2 = rhigh*rhigh;
double rlow2 = rlow*rlow;
T b;

v /= v.length();

// standard coordinates so cone is always centered on the x-axis
Point_3d y_axis,z_axis;

y_axis[0] = -v[1];
y_axis[1] =  v[0];
y_axis[2] =  0;

y_axis /= y_axis.length();
z_axis = v.cross(y_axis);

//std::cout << v*y_axis << " " << v*z_axis << " " << y_axis*z_axis << std::endl;

// we can do better than this !!!!
int n0[2] = {(int)((rhigh + xo[0])/Length),(int)((xo[0] - rhigh)/Length) - 1 };
int n1[2] = {(int)((rhigh + xo[1])/Length),(int)((xo[1] - rhigh)/Length) - 1 };
int n2[2] = {(int)((rhigh + xo[2])/Length),(int)((xo[2] - rhigh)/Length) - 1 };
if(!periodic_boundaries){
  n0[0] = n0[1] = 0;
  n1[0] = n1[1] = 0;
  n2[0] = n2[1] = 0;
}
for(auto a = begin ; a != end ; ++a){
  dx[0] = (*a).x[0] - xo[0];
  dx[1] = (*a).x[1] - xo[1];
  dx[2] = (*a).x[2] - xo[2];

  int count = 0;

  for(n[0] = n0[1]; n[0] <= n0[0] ; ++n[0]){
    for(n[1] = n1[1]; n[1] <= n1[0] ; ++n[1]){
      for(n[2] = n2[1]; n[2] <= n2[0] ; ++n[2]){

        x = dx + n*Length;
        xp = x*v;
        if(xp > 0){
          r2 = x.length_sqr();

          if(r2 > rlow2 && r2 < rhigh2){
            if( r2*sin_theta_sqrt > r2 - xp*xp){
              //std::cout << n << " | " << xo << " | "
              //          << v << " |  " << x << std::endl;

              b = *a;
              //b.x = x;
              // rotate to standard reference frame
              b.x[0] = xp;
              b.x[1] = y_axis*x;
              b.x[2] = z_axis*x;

              //std::cout << "b = " << b.x << std::endl;

              incone.push_back(b);
              ++count;
            }
          }
        }
      }
    }
  }
  if(count > 1) std::cout << " Warning: Same halo appears " << count << " times in light-cone." << std::endl;
}

}

class MultiLightCone{ public: MultiLightCone( double angular_radius ,const std::vector &observers /// postion of observers within the simulation box ,const std::vector &directions /// direction of light cones ): xos(observers),vs(directions) { assert(observers.size() == directions.size()); for(int i=0;i<observers.size(); ++i) cones.push_back(angular_radius); }

/** Read the points in from a snapshot of the halos created by RockStar.
 The output will be the halos in the cone in coordinates where the x-axis is
 the center of the cone.  Changing the observer, xo, and direction of view, V,
 effectively translates and rotates the box, respectively.  Only halos with radii
 between rlow and rhigh are added.
 */
void ReadBoxRockStar(std::string filename
                     ,double rlow,double rhigh
                     ,std::vector<std::vector<LightCones::DataRockStar> > &conehalos
                     ,bool periodic_boundaries = true
                     ,bool allow_subhalos = false);

void ReadBoxXYZ(std::string filename
                ,double rlow,double rhigh
                ,std::vector<std::vector<Point_3d> > &conehalos
                ,double hubble
                ,double BoxLength
                ,bool periodic_boundaries = true
                );

void ReadBoxToMaps(std::string filename
                   ,double rlow,double rhigh
                   ,std::vector<std::vector<PixelMap>  > &maps
                   ,double hubble
                   ,double BoxLength
                   ,bool periodic_boundaries = true
                   );

/* double common_scan_blockH5(H5std_string filename
            , H5std_string dataset_name
            , hsize_t n_rows
            , hsize_t n_cols
            , hsize_t offset_rows
            , hsize_t *read_rows
            ); */

private:

std::vector<Point_3d> xos;
std::vector<Point_3d> vs;
std::vector<LightCone> cones;
/* double common_scan_blockH5(H5std_string filename
                            , H5std_string dataset_name
                            , hsize_t n_rows 
                            , hsize_t n_cols  
                            , hsize_t offset_rows 
                            , hsize_t *read_rows
                            ); */

};

using Utilities::Geometry::Quaternion; using Utilities::Geometry::SphericalPoint;

/ templated functions for projecting mass onto planes /

/ templated functions for reading a block from a file /

template size_t scan_block(size_t blocksize,std::vector &points,FILE *pFile){ double tmpf; // read in a block of points size_t i=0;

points.resize(blocksize);
while(i < blocksize &&
      fscanf(pFile,"%lf %lf %lf %lf %lf %lf"
             ,&points[i][0],&points[i][1],&points[i][2],&tmpf,&tmpf,&tmpf) != EOF)
  ++i;
points.resize(i);

return i;

}

/ structures to implement templated FastLightCones<>() /

struct ASCII_XV{

std::vector<Point_3d> points;

size_t scan_block(size_t blocksize,FILE *pFile){
  double tmpf;
  // read in a block of points
  size_t i=0;

  points.resize(blocksize);
  while(i < blocksize &&
        fscanf(pFile,"%lf %lf %lf %lf %lf %lf"
               ,&points[i][0],&points[i][1],&points[i][2],&tmpf,&tmpf,&tmpf) != EOF)
    ++i;
  points.resize(i);

  return i;
}

static void fastplanes_parallel(
                                Point_3d *begin
                                ,Point_3d *end
                                ,const COSMOLOGY &cosmo
                                ,std::vector<std::vector<Point_3d> > &boxes
                                ,std::vector<Point_3d> &observers
                                ,std::vector<Quaternion> &rotationQs
                                ,std::vector<double> &dsources
                                ,std::vector<std::vector<PixelMap> > &maps
                                ,double dmin
                                ,double dmax
                                ,double BoxLength
                                );

}; struct ASCII_XM{

std::vector<DatumXM> points;
size_t scan_block(size_t blocksize,FILE *pFile){

  // read in a block of points
  size_t i=0;

  points.resize(blocksize);
  while(i < blocksize &&
        fscanf(pFile,"%lf %lf %lf %lf"
               ,&points[i].x[0],&points[i].x[1],&points[i].x[2],&points[i].mass) != EOF)
    ++i;
  points.resize(i);

  return i;
}
static void fastplanes_parallel(
                                DatumXM *begin
                                ,DatumXM *end
                                ,const COSMOLOGY &cosmo
                                ,std::vector<std::vector<Point_3d> > &boxes
                                ,std::vector<Point_3d> &observers
                                ,std::vector<Quaternion> &rotationQs
                                ,std::vector<double> &dsources
                                ,std::vector<std::vector<PixelMap> > &maps
                                ,double dmin
                                ,double dmax
                                ,double BoxLength
                                );

};

struct ASCII_XMR{

std::vector<DatumXMR> points;

size_t scan_block(size_t blocksize,FILE *pFile){

  // read in a block of points
  size_t i=0;

  points.resize(blocksize);
  while(i < blocksize &&
        fscanf(pFile,"%lf %lf %lf %lf %lf"
               ,&points[i].x[0],&points[i].x[1],&points[i].x[2]
               ,&points[i].mass,&points[i].r) != EOF)
    ++i;
  points.resize(i);

  return i;
}
static void fastplanes_parallel(
                                DatumXMR *begin
                                ,DatumXMR *end
                                ,const COSMOLOGY &cosmo
                                ,std::vector<std::vector<Point_3d> > &boxes
                                ,std::vector<Point_3d> &observers
                                ,std::vector<Quaternion> &rotationQs
                                ,std::vector<double> &dsources
                                ,std::vector<std::vector<PixelMap> > &maps
                                ,double dmin
                                ,double dmax
                                ,double BoxLength
                                );

//CHANGED ********************************************************* MARCOS
/*size_t scan_blockHDF5(size_t blocksize,H5std_string filename){
  // read in a block of points
  size_t i=0, j=0;
  int k=5;
  hsize_t nrows, ncols=7, offset=0;
  double *data;
  data = common_scan_blockH5(filename, "/dset", blocksize, ncols, offset, &nrows);
  points.resize(nrows);
  for (i = 0; i < nrows; i++) {
if (((size_t)data[i*ncols+k] == 1) || ((size_t)data[i*ncols+k] == 2)) {
  points[j].x[0]    = data[i*ncols  ];
  points[j].x[1]    = data[i*ncols+1];
  points[j].x[2]    = data[i*ncols+2];
  points[j].mass    = data[i*ncols+3];
  points[j].r_max   = data[i*ncols+4];
  //points[j].r_scale = data[i*ncols+5];
  j++;
}
  }
  points.resize(j);
  delete data;

  for(auto &h: points){
h.r_max /= 1.0e3;
h.r_scale /= 1.0e3;
  }
  return j;
  }*/

};

struct ASCII_XMRRT{

std::vector<DatumXMRmRs> points;
size_t scan_block(size_t blocksize,FILE *pFile){

  // read in a block of points
  size_t i=0;
  int tmp;

  points.resize(blocksize);
  while(i < blocksize &&
        fscanf(pFile,"%lf %lf %lf %lf %lf %lf %i"
               ,&points[i].x[0],&points[i].x[1],&points[i].x[2]
               ,&points[i].mass,&points[i].r_max,&points[i].r_scale,&tmp) != EOF)
    ++i;
  points.resize(i);

  for(auto &h: points){
    h.r_max /= 1.0e3;
    h.r_scale /= 1.0e3;
  }
  return i;
}

static void fastplanes_parallel(
                                DatumXMRmRs *begin
                                ,DatumXMRmRs *end
                                ,const COSMOLOGY &cosmo
                                ,std::vector<std::vector<Point_3d> > &boxes
                                ,std::vector<Point_3d> &observers
                                ,std::vector<Quaternion> &rotationQs
                                ,std::vector<double> &dsources
                                ,std::vector<std::vector<PixelMap> > &maps
                                ,double dmin
                                ,double dmax
                                ,double BoxLength
                                );

};

struct ASCII_XMRRT12:public ASCII_XMRRT{ size_t scan_block(size_t blocksize,FILE *pFile){

  // read in a block of points
  size_t i=0;
  int tmp;

  points.resize(blocksize);
  while(i < blocksize &&
        fscanf(pFile,"%lf %lf %lf %lf %lf %lf %i"
               ,&points[i].x[0],&points[i].x[1],&points[i].x[2]
               ,&points[i].mass,&points[i].r_max,&points[i].r_scale,&tmp) != EOF){
          if(tmp == 1 || tmp == 2) ++i;
        }
  points.resize(i);

  for(auto &h: points){
    h.r_max /= 1.0e3;
    h.r_scale /= 1.0e3;
  }
  return i;
}

//CHANGED ********************************************************* MARCOS
size_t scan_blockH5(size_t blocksize,H5std_string filename){
  // read in a block of points
  size_t i=0, j=0;
  int k=6;
  hsize_t nrows, ncols=7, offset=0;
  double *data;
  data = common_scan_blockH5(filename, "/dset", blocksize, ncols, offset, &nrows);
  points.resize(nrows);
  for (i = 0; i < nrows; i++) {
    if (((size_t)data[i*ncols+k] == 1) || ((size_t)data[i*ncols+k] == 2)) {
      points[j].x[0]    = data[i*ncols  ];
      points[j].x[1]    = data[i*ncols+1];
      points[j].x[2]    = data[i*ncols+2];
      points[j].mass    = data[i*ncols+3];
      points[j].r_max   = data[i*ncols+4];
      points[j].r_scale = data[i*ncols+5]; 
      j++;
    }
  }
  points.resize(j); 
  delete data;

  for(auto &h: points){
h.r_max /= 1.0e3;
h.r_scale /= 1.0e3;
  }
  return j;
}

double *common_scan_blockH5(H5std_string filename
               , H5std_string dataset_name
               , hsize_t n_rows
               , hsize_t n_cols
               , hsize_t offset_rows
               , hsize_t *read_rows
               );

}; //CHANGED ***** MARCOS

/ /

void random_observers(std::vector &observers ,std::vector &directions ,int Ncones ,double BoxLength ,double cone_opening_radius ,Utilities::RandomNumbers_NR &ran );

struct DkappaDz{ DkappaDz(const COSMOLOGY &cos,double zsource):cosmo(cos),zs(zsource){ rho = cosmo.getOmega_matter()*cosmo.rho_crit(0); };

double operator()(double z){
  double x = 1+z;
  return cosmo.drdz(x)*rho*x*x*x/cosmo.SigmaCrit(z,zs)/cosmo.getHubble();
}

const COSMOLOGY &cosmo;
double zs;
double rho;

};

using Utilities::Geometry::Quaternion; using Utilities::Geometry::SphericalPoint; /** \brief This function goes directly from snapshots to lensing maps with Born approximation and linear propogations.

The template parameter allows this function to use different input data formats and different ways of distributing the mass into grid cells. The current options are LightCones::ASCII_XV -- for 6 column ASCII file with position and velocity LightCones::ASCII_XM -- for 5 column ASCII file with position and mass LightCones::ASCII_XMR -- for 6 column ASCII file with position, mass and size LightCones::ASCII_XMRRT -- for 7 column ASCII file with position, mass,Rmax,Rscale and an integer denoting type, the halos are rendered as truncated NFWs LightCone::ASCII_XMRRT12 -- same as LightCones::ASCII_XMRRT but only takes entries with the 7th column equal to 1 or 2 <\p> */ template void FastLightCones( const COSMOLOGY &cosmo ,const std::vector &zsources /// vector of source redshifts ,std::vector > &maps /// output kappa maps, do not need to be allocated initially ,double range /// angular range in radians of maps ,double angular_resolution /// angular resolution in radians of maps ,std::vector &observers /// position of observers within the simulation box coordinates 0 < x < Lengths ,std::vector &directions /// direction of light cones ,const std::vector &snap_filenames /// names of the data files ,const std::vector &snap_redshifts /// the redshift of the data files ,double BoxLength ,double mass_units ,bool verbose = false ,bool addtocone = false /// if false the maps will be cleared and new maps made, if true particles are added to the existing maps ){ assert(cosmo.getOmega_matter() + cosmo.getOmega_lambda() == 1.0); // set coordinate distances for planes int Nmaps = zsources.size(); // number of source planes per cone int Ncones = observers.size(); // number of cones if(directions.size() != observers.size()){ std::cerr << "Size of direction and observers must match." << std::endl; throw std::invalid_argument(""); } double zs_max=0; for(auto z: zsources) zs_max = (z > zs_max) ? z : zs_max; Point_2d center; if(!addtocone){ // allocate mamory for all maps maps.clear(); maps.resize(Ncones); for(auto &map_v : maps){ // loop through cones map_v.reserve(Nmaps); for(int i=0;i > > map_pack(Utilities::GetNThreads()); for(auto &maps : map_pack){ maps.resize(Ncones); for(auto &map_v : maps){ // loop through cones map_v.reserve(Nmaps); for(int i=0;i z_unique(1,snap_redshifts[0]); for(auto z : snap_redshifts) if(z != z_unique.back() ) z_unique.push_back(z); std::sort(z_unique.begin(),z_unique.end()); // find redshift ranges for each snapshot // shift the redshifts to between the snapshots std::vector abins(z_unique.size() + 1); abins[0] = 1.0; for(int i=1;i dbins(abins.size()); for(int i=0 ; i dsources(zsources.size()); for(int i=0 ; i rotationQs(Ncones); for(int i = 0 ; i > boxes(Ncones); for(int icone=0;icone dmin || p2.length() > dmin || (p1 + Point_3d(BoxLength,0,0)).length() > dmin || (p1 + Point_3d(0,BoxLength,0)).length() > dmin || (p1 + Point_3d(0,0,BoxLength)).length() > dmin || (p1 + Point_3d(0,BoxLength,BoxLength)).length() > dmin || (p1 + Point_3d(BoxLength,0,BoxLength)).length() > dmin || (p1 + Point_3d(BoxLength,BoxLength,0)).length() > dmin ) boxes[icone].push_back(n); } } } } } //open file if(verbose) std::cout <<" Opening " << snap_filenames[i_file] << std::endl; FILE *pFile = fopen(snap_filenames[i_file].c_str(),"r"); if(pFile == nullptr){ std::cout << "Can't open file " << snap_filenames[i_file] << std::endl; ERROR_MESSAGE(); throw std::runtime_error(" Cannot open file."); } //points.resize(blocksize); size_t Nlines = 0,Nblocks=0; while(!feof(pFile)){ // loop through blocks //long Nbatch = unit.scan_block(blocksize,pFile); long Nbatch = unit.scan_blockH5(blocksize,snap_filenames[i_file]); if(Nbatch > 0){ // multi-thread the sorting into cones and projection onto planes int nthreads = Utilities::GetNThreads(); int chunk_size; do{ chunk_size = unit.points.size()/nthreads; if(chunk_size == 0) nthreads /= 2; }while(chunk_size == 0); if(nthreads == 0) nthreads = 1; int remainder = unit.points.size()%chunk_size; assert(nthreads*chunk_size + remainder == unit.points.size() ); std::vector thr(nthreads); for(int ii =0; ii< nthreads ; ++ii){ //std::cout << ii*chunk_size << " " << n << std::endl; thr[ii] = std::thread(&unit.fastplanes_parallel ,unit.points.data() + ii*chunk_size ,unit.points.data() + (ii+1)*chunk_size + (ii==nthreads-1)*remainder ,cosmo,std::ref(boxes) ,std::ref(observers),std::ref(rotationQs) ,std::ref(dsources),std::ref(map_pack[ii]) ,dmin,dmax,BoxLength); } for(int ii = 0; ii < nthreads ;++ii){ if(thr[ii].joinable() ) thr[ii].join();} } Nlines += blocksize; ++Nblocks; std::cout << "=" << std::flush ; if(Nblocks%10 == 0) std::cout << " " << std::flush; if(Nblocks%50 == 0) std::cout << std::endl; if(Nblocks%100 == 0) std::cout << std::endl; //std::cout << std::endl; } fclose(pFile); std::cout << std::endl; } for(auto &tmaps : map_pack){ for(int isource=0;isource(dkappadz,0,zsources[isource],1.0e-6); for(int icone=0 ; icone

Hi Ben,

I finally have a working code inside glamer that reads HDF5 in chuncks of a million (that can be changed) where I followed your advice and used another class for reading this:

Konsole output structHDF5_XMRRT12:public ASCII_XMRRT{

using scan_block(). Because of the structure of how H5 is coded I needed to change a bit the arguments I passed to scan_block() so now I am passing

Konsole output scan_block(size_tblocksize,FILEpFile,H5std_stringfilename,hsize_toffset_rows,boolH5eof)

I made sure that it still works in the other ACII_ structures so I guess this is not a big problem (it might give some warnings because we are not using some arguments when reading from ASCII but that won't be a problem, I agree it might not be the most beautiful way).

Something that I don't like is that I call a function (common_scan_blockH5() that reads in H5 format a set of data) from the structure HDF5_XMRRT12 and, as I didn't know where to put it, it is now at lightcones_construction.h, just before overriding scan_block(). Having a function at a .h file is very ugly. Could you tell me where can I put that function so that lightcones_construction sees it and can call it?

Summarizing, I've got a code, and it works, but it is still ugly, it needs some polishing. Can I "git pull" it? Is it as simple as typing "git pull" at the SLsimlib folder? I am attaching the code anyway.

Cheers! Marcos.

** On 07/07/2017 11:59 AM, Nicolas Tessore wrote:

Ok, to test this, try and do the following three steps:

1.
update and checkout the test branch for the *glamer* folder:

|$ cd glamer $ git pull && git checkout nt/hdf5 |
2.
update and checkout the test branch for the *SLsimLib* folder:

|$ cd SLsimLib $ git pull && git checkout nt/hdf5 |
3.
add the (unfortunately necessary) definitions to *your project*'s
CMakeLists.txt:

|# somewhere in file CMakeLists.txt ... find_package(GLAMER
NO_MODULE REQUIRED) include_directories(${GLAMER_INCLUDE_DIRS})
add_definitions(${GLAMER_DEFINITIONS}) # <-- this line is new ... |
Then you should be able to use |cmake .. -DENABLE_HDF5| when building GLAMER to find the HDF5 library. Subsequently building your project should no longer need any intervention.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/glenco/SLsimLib/issues/131#issuecomment-313651746, or mute the thread https://github.com/notifications/unsubscribe-auth/AVsIK4ry4UjmcxU430IeDXJf2cdlFhR6ks5sLg9_gaJpZM4OMpVb.

// // lightcone_construction.h // GLAMER // // Created by Ben Metcalf on 26/10/16. // //

ifndef GLAMERlightcone_construction__

define GLAMERlightcone_construction__

include

include "point.h"

include "lens_halos.h"

include "geometry.h"

ifdef OLD_HEADER_FILENAME

include

else

include

endif

using std::cout; using std::endl;

include

include "H5Cpp.h"

using namespace H5;

/* \brief The LightCones namespace is for classes and functions related to making light cones and weak lensing maps from 3D snapshots of particles or halos. / namespace LightCones{

struct DataRockStar{ unsigned int id; Point_3d x; double mass; double Rvir; // in comoving units double Rscale; // in comoving units double Vmax; };

struct DatumXM{ Point_3d x; double mass; }; struct DatumXMRmRs{ Point_3d x; double mass; double r_max; double r_scale; }; struct DatumXMR{ Point_3d x; double mass; double r; };

void ReadLightConeNFW(std::string filename,COSMOLOGY &cosmo ,std::vector<LensHalo* > &lensVec ,PosType &theta_max);

void ReadLightConeParticles(std::string filename,COSMOLOGY &cosmo ,std::vector<LensHaloMassMap *> &lensVec ,int Nplanes ,float particle_mass ,float angular_resolution ,bool verbose = false);

/** \brief Class for constructing light cones form simulation boxes. / class LightCone{ public: LightCone(double angular_radius);

void WriteLightCone(std::string filename,std::vector<DataRockStar> &vec);
void WriteLightCone(std::string filename,std::vector<Point_3d> &vec);

friend class MultiLightCone;

template <typename T>
void select(Point_3d xo,Point_3d v,double Length,double rlow,double rhigh
            ,T* begin
            ,T* end
            ,Utilities::LockableContainer<std::vector<T> > &incone
            ,bool periodic_boundaries = true);

private:

double r_theta; // angular radius of light cone (radians)
double sin_theta_sqrt;

if(begin == end) return;

Point_3d dx,x;
double r2,xp;
Point_3d n;
double rhigh2 = rhigh*rhigh;
double rlow2 = rlow*rlow;
T b;

v /= v.length();

// standard coordinates so cone is always centered on the x-axis
Point_3d y_axis,z_axis;

y_axis[0] = -v[1];
y_axis[1] =  v[0];
y_axis[2] =  0;

y_axis /= y_axis.length();
z_axis = v.cross(y_axis);

//std::cout << v*y_axis << " " << v*z_axis << " " << y_axis*z_axis << std::endl;

// we can do better than this !!!!
int n0[2] = {(int)((rhigh + xo[0])/Length),(int)((xo[0] - rhigh)/Length) - 1 };
int n1[2] = {(int)((rhigh + xo[1])/Length),(int)((xo[1] - rhigh)/Length) - 1 };
int n2[2] = {(int)((rhigh + xo[2])/Length),(int)((xo[2] - rhigh)/Length) - 1 };
if(!periodic_boundaries){
  n0[0] = n0[1] = 0;
  n1[0] = n1[1] = 0;
  n2[0] = n2[1] = 0;
}
for(auto a = begin ; a != end ; ++a){
  dx[0] = (*a).x[0] - xo[0];
  dx[1] = (*a).x[1] - xo[1];
  dx[2] = (*a).x[2] - xo[2];

  int count = 0;

  for(n[0] = n0[1]; n[0] <= n0[0] ; ++n[0]){
    for(n[1] = n1[1]; n[1] <= n1[0] ; ++n[1]){
      for(n[2] = n2[1]; n[2] <= n2[0] ; ++n[2]){

        x = dx + n*Length;
        xp = x*v;
        if(xp > 0){
          r2 = x.length_sqr();

          if(r2 > rlow2 && r2 < rhigh2){
            if( r2*sin_theta_sqrt > r2 - xp*xp){
              //std::cout << n << " | " << xo << " | "
              //          << v << " |  " << x << std::endl;

              b = *a;
              //b.x = x;
              // rotate to standard reference frame
              b.x[0] = xp;
              b.x[1] = y_axis*x;
              b.x[2] = z_axis*x;

              //std::cout << "b = " << b.x << std::endl;

              incone.push_back(b);
              ++count;
            }
          }
        }
      }
    }
  }
  if(count > 1) std::cout << " Warning: Same halo appears " << count << " times in light-cone." << std::endl;
}

}

/** Read the points in from a snapshot of the halos created by RockStar.
 The output will be the halos in the cone in coordinates where the x-axis is
 the center of the cone.  Changing the observer, xo, and direction of view, V,
 effectively translates and rotates the box, respectively.  Only halos with radii
 between rlow and rhigh are added.
 */
void ReadBoxRockStar(std::string filename
                     ,double rlow,double rhigh
                     ,std::vector<std::vector<LightCones::DataRockStar> > &conehalos
                     ,bool periodic_boundaries = true
                     ,bool allow_subhalos = false);

void ReadBoxXYZ(std::string filename
                ,double rlow,double rhigh
                ,std::vector<std::vector<Point_3d> > &conehalos
                ,double hubble
                ,double BoxLength
                ,bool periodic_boundaries = true
                );

void ReadBoxToMaps(std::string filename
                   ,double rlow,double rhigh
                   ,std::vector<std::vector<PixelMap>  > &maps
                   ,double hubble
                   ,double BoxLength
                   ,bool periodic_boundaries = true
                   );

private:

std::vector<Point_3d> xos;
std::vector<Point_3d> vs;
std::vector<LightCone> cones;

};

using Utilities::Geometry::Quaternion; using Utilities::Geometry::SphericalPoint;

/ templated functions for projecting mass onto planes /

/ templated functions for reading a block from a file /

template size_t scan_block(size_t blocksize,std::vector &points,FILE *pFile){ double tmpf; // read in a block of points size_t i=0;

points.resize(blocksize);
while(i < blocksize &&
      fscanf(pFile,"%lf %lf %lf %lf %lf %lf"
             ,&points[i][0],&points[i][1],&points[i][2],&tmpf,&tmpf,&tmpf) != EOF)
  ++i;
points.resize(i);

return i;

}

/ structures to implement templated FastLightCones<>() /

struct ASCII_XV{

std::vector<Point_3d> points;

size_t scan_block(size_t blocksize,FILE *pFile,H5std_string filename,hsize_t offset_rows,bool *H5eof){
  double tmpf;
  // read in a block of points
  size_t i=0;

  points.resize(blocksize);
  while(i < blocksize &&
        fscanf(pFile,"%lf %lf %lf %lf %lf %lf"
               ,&points[i][0],&points[i][1],&points[i][2],&tmpf,&tmpf,&tmpf) != EOF)
    ++i;
  points.resize(i);

  return i;
}

static void fastplanes_parallel(
                                Point_3d *begin
                                ,Point_3d *end
                                ,const COSMOLOGY &cosmo
                                ,std::vector<std::vector<Point_3d> > &boxes
                                ,std::vector<Point_3d> &observers
                                ,std::vector<Quaternion> &rotationQs
                                ,std::vector<double> &dsources
                                ,std::vector<std::vector<PixelMap> > &maps
                                ,double dmin
                                ,double dmax
                                ,double BoxLength
                                );

}; struct ASCII_XM{

std::vector<DatumXM> points;
size_t scan_block(size_t blocksize,FILE *pFile,H5std_string filename,hsize_t offset_rows,bool *H5eof){

  // read in a block of points
  size_t i=0;

  points.resize(blocksize);
  while(i < blocksize &&
        fscanf(pFile,"%lf %lf %lf %lf"
               ,&points[i].x[0],&points[i].x[1],&points[i].x[2],&points[i].mass) != EOF)
    ++i;
  points.resize(i);

  return i;
}
static void fastplanes_parallel(
                                DatumXM *begin
                                ,DatumXM *end
                                ,const COSMOLOGY &cosmo
                                ,std::vector<std::vector<Point_3d> > &boxes
                                ,std::vector<Point_3d> &observers
                                ,std::vector<Quaternion> &rotationQs
                                ,std::vector<double> &dsources
                                ,std::vector<std::vector<PixelMap> > &maps
                                ,double dmin
                                ,double dmax
                                ,double BoxLength
                                );

};

struct ASCII_XMR{

std::vector<DatumXMR> points;

size_t scan_block(size_t blocksize,FILE *pFile,H5std_string filename,hsize_t offset_rows,bool *H5eof){

  // read in a block of points
  size_t i=0;

  points.resize(blocksize);
  while(i < blocksize &&
        fscanf(pFile,"%lf %lf %lf %lf %lf"
               ,&points[i].x[0],&points[i].x[1],&points[i].x[2]
               ,&points[i].mass,&points[i].r) != EOF)
    ++i;
  points.resize(i);

  return i;
}
static void fastplanes_parallel(
                                DatumXMR *begin
                                ,DatumXMR *end
                                ,const COSMOLOGY &cosmo
                                ,std::vector<std::vector<Point_3d> > &boxes
                                ,std::vector<Point_3d> &observers
                                ,std::vector<Quaternion> &rotationQs
                                ,std::vector<double> &dsources
                                ,std::vector<std::vector<PixelMap> > &maps
                                ,double dmin
                                ,double dmax
                                ,double BoxLength
                                );

};

struct ASCII_XMRRT{

std::vector<DatumXMRmRs> points;
size_t scan_block(size_t blocksize,FILE *pFile,H5std_string filename,hsize_t offset_rows,bool *H5eof){

  // read in a block of points
  size_t i=0;
  int tmp;

  points.resize(blocksize);
  while(i < blocksize &&
        fscanf(pFile,"%lf %lf %lf %lf %lf %lf %i"
               ,&points[i].x[0],&points[i].x[1],&points[i].x[2]
               ,&points[i].mass,&points[i].r_max,&points[i].r_scale,&tmp) != EOF)
    ++i;
  points.resize(i);

  for(auto &h: points){
    h.r_max /= 1.0e3;
    h.r_scale /= 1.0e3;
  }
  return i;
}

static void fastplanes_parallel(
                                DatumXMRmRs *begin
                                ,DatumXMRmRs *end
                                ,const COSMOLOGY &cosmo
                                ,std::vector<std::vector<Point_3d> > &boxes
                                ,std::vector<Point_3d> &observers
                                ,std::vector<Quaternion> &rotationQs
                                ,std::vector<double> &dsources
                                ,std::vector<std::vector<PixelMap> > &maps
                                ,double dmin
                                ,double dmax
                                ,double BoxLength
                                );

};

struct ASCII_XMRRT12:public ASCII_XMRRT{ size_t scan_block(size_t blocksize,FILE pFile,H5std_string filename,hsize_t offset_rows,bool H5eof){

  // read in a block of points
  size_t i=0;
  int tmp;

  points.resize(blocksize);
  while(i < blocksize &&
        fscanf(pFile,"%lf %lf %lf %lf %lf %lf %i"
               ,&points[i].x[0],&points[i].x[1],&points[i].x[2]
               ,&points[i].mass,&points[i].r_max,&points[i].r_scale,&tmp) != EOF){
          if(tmp == 1 || tmp == 2) ++i;
        }
  points.resize(i);

  for(auto &h: points){
    h.r_max /= 1.0e3;
    h.r_scale /= 1.0e3;
  }
  return i;
}

};

struct HDF5_XMRRT12:public ASCII_XMRRT{

// added by Marcos Pellejero & Antonio Dorta*** /Modified by Marcos and Antonio 4th July 2017/

double common_scan_blockH5(H5std_string filename, H5std_string dataset_name, hsize_t n_rows, hsize_t n_cols, hsize_t offset_rows, hsize_t read_rows) { const int RANK = 2; double data_out = NULL; read_rows = 0; // Try block to detect exceptions raised by any of the calls inside it
try { /*

Turn off the auto-printing when failure occurs so that we can
handle the errors appropriately
/ Exception::dontPrint(); /
Open the file and the dataset.
/ H5File file( filename, H5F_ACC_RDONLY ); DataSet dataset = file.openDataSet(dataset_name); /
Get filespace for rank and dimension
/ DataSpace filespace = dataset.getSpace(); /
Get number of dimensions in the file dataspace
/ int rank = filespace.getSimpleExtentNdims(); /
Get and print the dimension sizes of the file dataspace
/ hsize_t dims[RANK]; // dataset dimensions
hsize_t count_rows; hsize_t offset[RANK]; rank = filespace.getSimpleExtentDims( dims ); cout << "dataset rank = " << rank << ", dimensions " << (hsize_t)(dims[0]) << " x " << (hsize_t)(dims[1]) << endl; /
Define the memory space to read dataset.
*/ hsize_t file_rows = dims[0], file_cols = dims[1]; bool readALL = false;

if (n_cols != file_cols) { cout << "ERROR!!! You specified " << n_cols << " cols, but file contains " << dims[1] << " cols.\n"; return data_out; }
```
if (offset_rows > file_rows) {
```
// Reading out of data
cout << "ERROR!!! File has " << file_rows << " rows and you wanted to read till row " << offset_rows + n_rows << "\n"; return data_out; } else if ((offset_rows == 0) && (file_rows >= n_rows)) { // Read ALL data
count_rows = n_rows; readALL = (count_rows == file_rows); cout << "Reading " << count_rows << "\n"; } else if (offset_rows + n_rows > file_rows) { // We are out of bound, limit the max rows
count_rows = file_rows - offset_rows; cout << "Reading LIMIT ROWS: " << count_rows << " FROM ROW " << offset_rows << "\n"; } else { // Inside bounds
count_rows = n_rows; cout << "Reading INSIDE\n"; }

/*
Read dataset back and display.
*/

// Allocate only the required memory data_out = new double [count_rowsfile_cols]; if (data_out == NULL) { cout << "NOT ENOUGH MEMORY!!!\n"; read_rows = 0; return NULL; }

// READ DATA!!! (we are going to read count_rows x file_cols) if (readALL) { DataSpace memspace(RANK, dims); dataset.read(data_out, PredType::NATIVE_DOUBLE, memspace, filespace); } else { hsize_t count[RANK]; count[0] = count_rows; count[1] = file_cols; offset[0] = offset_rows; offset[1] = 0; DataSpace memspace(RANK, count); filespace.selectHyperslab(H5S_SELECT_SET, count, offset); dataset.read(data_out, PredType::NATIVE_DOUBLE, memspace, filespace); } *read_rows = count_rows;

} // end of try block // catch failure caused by the H5File operations catch( FileIException error ) { error.printError(); read_rows = 0; return NULL; } // catch failure caused by the DataSet operations catch( DataSetIException error ) { error.printError(); read_rows = 0; return NULL; } // catch failure caused by the DataSpace operations catch( DataSpaceIException error ) { error.printError(); *read_rows = 0; return NULL; } return data_out; }

//CHANGED ***** MARCOS // read in a block of points size_t scan_block(size_t blocksize,FILE pFile,H5std_string filename,hsize_t offset_rows,bool H5eof){ // HDF5 needs filename // Offset (rows) // Since it does NOT use pointer to file, we need to tell when the EOF is reached. // We could return -1 when reaching EOF, but since site_t could be unsigned, that would NOT work // Instead of that, we return H5eof -> true when EOF size_t i=0, j=0; int k=6; H5eof = false; hsize_t nrows, ncols=7; double data; cout << "Reading " << blocksize << " from offset " << offset_rows << "\n"; H5eof = false; data = common_scan_blockH5(filename, "/dset", blocksize, ncols, offset_rows, &nrows); if ((data == NULL) || (nrows == 0)) { H5eof = true; return 0; } points.resize(nrows); for (i = 0; i < nrows; i++) { if (((size_t)data[incols+k] == 1) || ((size_t)data[incols+k] == 2)) { points[j].x[0] = data[incols ]; points[j].x[1] = data[incols+1]; points[j].x[2] = data[incols+2]; points[j].mass = data[incols+3]; points[j].r_max = data[incols+4]; points[j].r_scale = data[incols+5]; j++; } } points.resize(j); delete data;

for(auto &h: points){ h.r_max /= 1.0e3; h.r_scale /= 1.0e3; } return j; }

}; //CHANGED ***** MARCOS

/ /

void random_observers(std::vector &observers ,std::vector &directions ,int Ncones ,double BoxLength ,double cone_opening_radius ,Utilities::RandomNumbers_NR &ran );

struct DkappaDz{ DkappaDz(const COSMOLOGY &cos,double zsource):cosmo(cos),zs(zsource){ rho = cosmo.getOmega_matter()*cosmo.rho_crit(0); };

double operator()(double z){ double x = 1+z; return cosmo.drdz(x)rhoxxx/cosmo.SigmaCrit(z,zs)/cosmo.getHubble(); }

const COSMOLOGY &cosmo; double zs; double rho; };

using Utilities::Geometry::Quaternion; using Utilities::Geometry::SphericalPoint; /** \brief This function goes directly from snapshots to lensing maps with Born approximation and linear propogations.

The template parameter allows this function to use different input data formats and different ways of distributing the mass into grid cells. The current options are LightCones::ASCII_XV -- for 6 column ASCII file with position and velocity LightCones::ASCII_XM -- for 5 column ASCII file with position and mass LightCones::ASCII_XMR -- for 6 column ASCII file with position, mass and size LightCones::ASCII_XMRRT -- for 7 column ASCII file with position, mass,Rmax,Rscale and an integer denoting type, the halos are rendered as truncated NFWs LightCone::ASCII_XMRRT12 -- same as LightCones::ASCII_XMRRT but only takes entries with the 7th column equal to 1 or 2 <\p> */ template void FastLightCones( const COSMOLOGY &cosmo ,const std::vector &zsources /// vector of source redshifts ,std::vector > &maps /// output kappa maps, do not need to be allocated initially ,double range /// angular range in radians of maps ,double angular_resolution /// angular resolution in radians of maps ,std::vector &observers /// position of observers within the simulation box coordinates 0 < x < Lengths ,std::vector &directions /// direction of light cones ,const std::vector &snap_filenames /// names of the data files ,const std::vector &snap_redshifts /// the redshift of the data files ,double BoxLength ,double mass_units ,bool verbose = false ,bool addtocone = false /// if false the maps will be cleared and new maps made, if true particles are added to the existing maps ){ assert(cosmo.getOmega_matter() + cosmo.getOmega_lambda() == 1.0); // set coordinate distances for planes int Nmaps = zsources.size(); // number of source planes per cone int Ncones = observers.size(); // number of cones if(directions.size() != observers.size()){ std::cerr << "Size of direction and observers must match." << std::endl; throw std::invalid_argument(""); } double zs_max=0; for(auto z: zsources) zs_max = (z > zs_max) ? z : zs_max; Point_2d center; if(!addtocone){ // allocate mamory for all maps maps.clear(); maps.resize(Ncones); for(auto &map_v : maps){ // loop through cones map_v.reserve(Nmaps); for(int i=0;i > > map_pack(Utilities::GetNThreads()); for(auto &maps : map_pack){ maps.resize(Ncones); for(auto &map_v : maps){ // loop through cones map_v.reserve(Nmaps); for(int i=0;i z_unique(1,snap_redshifts[0]); for(auto z : snap_redshifts) if(z != z_unique.back() ) z_unique.push_back(z); std::sort(z_unique.begin(),z_unique.end()); // find redshift ranges for each snapshot // shift the redshifts to between the snapshots std::vector abins(z_unique.size() + 1); abins[0] = 1.0; for(int i=1;i dbins(abins.size()); for(int i=0 ; i dsources(zsources.size()); for(int i=0 ; i rotationQs(Ncones); for(int i = 0 ; i > boxes(Ncones); for(int icone=0;icone dmin || p2.length() > dmin || (p1 + Point_3d(BoxLength,0,0)).length() > dmin || (p1 + Point_3d(0,BoxLength,0)).length() > dmin || (p1 + Point_3d(0,0,BoxLength)).length() > dmin || (p1 + Point_3d(0,BoxLength,BoxLength)).length() > dmin || (p1 + Point_3d(BoxLength,0,BoxLength)).length() > dmin || (p1 + Point_3d(BoxLength,BoxLength,0)).length() > dmin ) boxes[icone].push_back(n); } } } } } //open file if(verbose) std::cout <<" Opening " << snap_filenames[i_file] << std::endl; FILE *pFile = fopen(snap_filenames[i_file].c_str(),"r"); if(pFile == nullptr){ std::cout << "Can't open file " << snap_filenames[i_file] << std::endl; ERROR_MESSAGE(); throw std::runtime_error(" Cannot open file."); } //points.resize(blocksize); size_t Nlines = 0,Nblocks=0; size_t Nbatch = 1; bool H5eof = false; while(!feof(pFile) && !H5eof) { // loop through blocks Nbatch = unit.scan_block(blocksize,pFile,snap_filenames[i_file],Nlines,&H5eof); cout << "READ " << Nbatch << " "<< H5eof << "!!\n"; //long Nbatch = unit.scan_blockH5(blocksize,snap_filenames[i_file]); if(Nbatch > 0){ // multi-thread the sorting into cones and projection onto planes int nthreads = Utilities::GetNThreads(); int chunk_size; do{ chunk_size = unit.points.size()/nthreads; if(chunk_size == 0) nthreads /= 2; }while(chunk_size == 0); if(nthreads == 0) nthreads = 1; int remainder = unit.points.size()%chunk_size; assert(nthreads*chunk_size + remainder == unit.points.size() ); std::vector thr(nthreads); for(int ii =0; ii< nthreads ; ++ii){ //std::cout << ii*chunk_size << " " << n << std::endl; thr[ii] = std::thread(&unit.fastplanes_parallel ,unit.points.data() + ii*chunk_size ,unit.points.data() + (ii+1)*chunk_size + (ii==nthreads-1)*remainder ,cosmo,std::ref(boxes) ,std::ref(observers),std::ref(rotationQs) ,std::ref(dsources),std::ref(map_pack[ii]) ,dmin,dmax,BoxLength); } for(int ii = 0; ii < nthreads ;++ii){ if(thr[ii].joinable() ) thr[ii].join();} } Nlines += blocksize; ++Nblocks; std::cout << "=" << std::flush ; if(Nblocks%10 == 0) std::cout << " " << std::flush; if(Nblocks%50 == 0) std::cout << std::endl; if(Nblocks%100 == 0) std::cout << std::endl; //std::cout << std::endl; } cout << "FILE READ!!\n"; fclose(pFile); std::cout << std::endl; } for(auto &tmaps : map_pack){ for(int isource=0;isource(dkappadz,0,zsources[isource],1.0e-6); for(int icone=0 ; icone

}

endif / defined(GLAMERlightcone_construction__) /

glenco / SLsimLib

hdf5 for FastLightCones #131

ifndef GLAMERlightcone_construction__

define GLAMERlightcone_construction__

include

include "point.h"

include "lens_halos.h"

include "geometry.h"

ifdef OLD_HEADER_FILENAME

include

else

include

endif

include

include "H5Cpp.h"

ifndef GLAMERlightcone_construction__

define GLAMERlightcone_construction__

include

include "point.h"

include "lens_halos.h"

include "geometry.h"

ifdef OLD_HEADER_FILENAME

include

else

include

endif

include

include "H5Cpp.h"

endif / defined(GLAMERlightcone_construction__) /