Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Loading Grids with properties #485

Merged
merged 11 commits into from
Jun 27, 2024
Merged

Loading Grids with properties #485

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
7930249
add getDescription function for GridProperties
JulienDoerner May 4, 2023
3e307e2
update interpolation type for grids
JulienDoerner May 4, 2023
331dd69
add possibility to store gridprop for TXT saving
JulienDoerner May 4, 2023
91f373c
loading function for Grid1f with Gridproperties
JulienDoerner May 4, 2023
5632e8f
loading function for Grid3f with Gridproperties
JulienDoerner May 4, 2023
6c9ff3e
add tests for grid loading with properties
JulienDoerner May 8, 2023
fd86a0a
Merge remote-tracking branch 'upstream/master' into GridProps
JulienDoerner Apr 23, 2024
5649102
add clip volume to reading function; update test for reading grid pro…
JulienDoerner Apr 23, 2024
e27102b
update example with grid property loading
JulienDoerner Apr 23, 2024
81c34f3
add missing throw statements
JulienDoerner Apr 23, 2024
a1adcde
update error messages.
JulienDoerner Jun 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
113 changes: 47 additions & 66 deletions doc/pages/example_notebooks/density/density_grid_sampling.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

31 changes: 31 additions & 0 deletions include/crpropa/Grid.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,8 @@

#include <vector>
#include <type_traits>
#include <string>
#include <sstream>
#if HAVE_SIMD
#include <immintrin.h>
#include <smmintrin.h>
Expand Down Expand Up @@ -128,6 +130,20 @@ class GridProperties: public Referenced {
void setClipVolume(bool b) {
clipVolume = b;
}

/** show all GridProperty parameters
* @param unit unit for the lengthscale (origin, spacing). Default is 1 = SI units
*/
std::string getDescription(double unit = 1) const {
std::stringstream ss;
ss << "GridProperties:\torigin: " << origin / unit
<< "\t" << "Nx: " << Nx << " Ny: " << Ny << " Nz: " << Nz
<< "\t" << "spacing: " << spacing / unit
<< "\t" << "refletive: " << reflective
<< "\t" << "interpolation: " << ipol
<< "\n";
return ss.str();
}
};

/**
Expand Down Expand Up @@ -247,6 +263,21 @@ class Grid: public Referenced {
}
}

interpolationType getInterpolationType() {
return ipolType;
}

std::string getInterpolationTypeName() {
if (ipolType == TRILINEAR)
return "TRILINEAR";
if (ipolType == TRICUBIC)
return "TRICUBIC";
if (ipolType == NEAREST_NEIGHBOUR)
return "NEAREST_NEIGHBOUR";

return "NOT_UNDERSTOOD";
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you add a more meaningful error message. Something like "ipolType not well defined: should be either TRILINEAR, TRICUBIC, or NEAREST_NEIGHBOUR".

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is not meant to be an error message. It is the conversion from the data-type to a string for saving.
To be consistent with the reading it has to be without spaces.

The alternative would be raising a runtime error at this point. In this case the saving would not run at all.

}

/** returns the position of the lower left front corner of the volume */
Vector3d getOrigin() const {
return origin;
Expand Down
18 changes: 16 additions & 2 deletions include/crpropa/GridTools.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -117,6 +117,12 @@ void dumpGrid(ref_ptr<Grid3f> grid, std::string filename,
void dumpGrid(ref_ptr<Grid1f> grid, std::string filename,
double conversion = 1);

/** Load a Grid3f from a plain text file based on the gridproperties stored in the header
@param filename name of the input file
@param conversion multiply every point in a grid by a conversion factor
*/
ref_ptr<Grid3f> loadGrid3fFromTxt(std::string filename, double conversion = 1);

/** Load a Grid3f grid from a plain text file.
@param grid a vector grid (Grid3f) to which the points will be loaded
@param filename name of input file
Expand All @@ -125,6 +131,12 @@ void dumpGrid(ref_ptr<Grid1f> grid, std::string filename,
void loadGridFromTxt(ref_ptr<Grid3f> grid, std::string filename,
double conversion = 1);

/** Load a Grid1f from a plain text file based on the gridproperties stored in the header
@param filename name of the input file
@param conversion multiply every point in a grid by a conversion factor
*/
ref_ptr<Grid1f> loadGrid1fFromTxt(std::string filename, double conversion = 1);

/** Load a Grid1f from a plain text file
@param grid a scalar grid (Grid1f) to which the points will be loaded
@param filename name of input file
Expand All @@ -137,17 +149,19 @@ void loadGridFromTxt(ref_ptr<Grid1f> grid, std::string filename,
@param grid a vector grid (Grid3f)
@param filename name of output file
@param conversion multiply every point in grid by a conversion factor
@param storeProperties if true the grid properties are stored as a comment
*/
void dumpGridToTxt(ref_ptr<Grid3f> grid, std::string filename,
double conversion = 1);
double conversion = 1, bool storeProperties = false);

/** Dump a Grid1f to a plain text file.
@param grid a scalar grid (Grid1f)
@param filename name of output file
@param conversion multiply every point in grid by a conversion factor
@param storeProperties if true the grid properties are stored as a comment
*/
void dumpGridToTxt(ref_ptr<Grid1f> grid, std::string filename,
double conversion = 1);
double conversion = 1, bool storeProperties = false);

#ifdef CRPROPA_HAVE_FFTW3F
/**
Expand Down
166 changes: 163 additions & 3 deletions src/GridTools.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -260,15 +260,84 @@ void loadGridFromTxt(ref_ptr<Grid3f> grid, std::string filename, double c) {
fin.close();
}

ref_ptr<Grid3f> loadGrid3fFromTxt(std::string filename, double c) {
std::ifstream fin(filename.c_str());
if (!fin) {
std::stringstream ss;
ss << "load Grid3f: " << filename << " not found";
throw std::runtime_error(ss.str());
}

// search in header lines for GridProperties
while (fin.peek() == '#') {
std::string line;
std::getline(fin, line);

// find gridproperties in the header
if (line.rfind("GridProperties:") == 2) {
GridProperties gp(Vector3d(0.), 1, 1, 1, 1.); // simple grid properties for default
std::stringstream ss(line);

// skip first names and check type
std::string name, type;
ss >> name >> name >> name >> type;
if (type != "Grid3f")
throw std::runtime_error("Tried to load Grid3f, but Gridproperties assume grid type " + type);

// grid origin
double x, y, z;
ss >> name >> x >> y >> z ;
gp.origin = Vector3d(x, y, z);

// grid size
ss >> name >> gp.Nx >> gp.Ny >> gp.Nz;

// spacing
double dX, dY, dZ;
ss >> name >> dX >> dY >> dZ;
gp.spacing = Vector3d(dX, dY, dZ);

// reflective
ss >> name >> gp.reflective;

// clip volume
bool clip;
ss >> name >> clip;
gp.setClipVolume(clip);

// interpolation type
ss >> name >> type;
if (type == "TRICUBIC")
gp.setInterpolationType(TRICUBIC);
else if (type == "NEAREST_NEIGHBOUR")
gp.setInterpolationType(NEAREST_NEIGHBOUR);
else
gp.setInterpolationType(TRILINEAR);

// create new grid
ref_ptr<Grid3f> grid = new Grid3f(gp);
fin.close();

// load data for grid
loadGridFromTxt(grid, filename, c);

return grid;
}
}
throw std::runtime_error("could not find GridProperties in file " + filename);
}


void loadGridFromTxt(ref_ptr<Grid1f> grid, std::string filename, double c) {
std::ifstream fin(filename.c_str());
if (!fin) {
std::stringstream ss;
ss << "load Grid1f: " << filename << " not found";
throw std::runtime_error(ss.str());
}

// skip header lines
while (fin.peek() == '#')
while (fin.peek() == '#')
fin.ignore(std::numeric_limits<std::streamsize>::max(), '\n');

for (int ix = 0; ix < grid->getNx(); ix++) {
Expand All @@ -285,13 +354,92 @@ void loadGridFromTxt(ref_ptr<Grid1f> grid, std::string filename, double c) {
fin.close();
}

void dumpGridToTxt(ref_ptr<Grid3f> grid, std::string filename, double c) {
ref_ptr<Grid1f> loadGrid1fFromTxt(std::string filename, double c) {
std::ifstream fin(filename.c_str());
if (!fin) {
std::stringstream ss;
ss << "load Grid1f: " << filename << " not found";
throw std::runtime_error(ss.str());
}

// search in header lines for GridProperties
while (fin.peek() == '#') {
std::string line;
std::getline(fin, line);

// find gridproperties in the header
if (line.rfind("GridProperties:") == 2) {
GridProperties gp(Vector3d(0.), 1, 1, 1, 1.); // simple grid properties for default
std::stringstream ss(line);

// skip first names and check type
std::string name, type;
ss >> name >> name >> name >> type;
if (type != "Grid1f")
throw std::runtime_error("Tried to load Grid1f, but Gridproperties assume grid type " + type);

// grid origin
double x, y, z;
ss >> name >> x >> y >> z ;
gp.origin = Vector3d(x, y, z);

// grid size
ss >> name >> gp.Nx >> gp.Ny >> gp.Nz;

// spacing
double dX, dY, dZ;
ss >> name >> dX >> dY >> dZ;
gp.spacing = Vector3d(dX, dY, dZ);

// reflective
ss >> name >> gp.reflective;

// clip volume
bool clip;
ss >> name >> clip;
gp.setClipVolume(clip);

// interpolation type
ss >> name >> type;
if (type == "TRICUBIC")
gp.setInterpolationType(TRICUBIC);
else if (type == "NEAREST_NEIGHBOUR")
gp.setInterpolationType(NEAREST_NEIGHBOUR);
else
gp.setInterpolationType(TRILINEAR);

// create new grid
ref_ptr<Grid1f> grid = new Grid1f(gp);
fin.close();

// load data for grid
loadGridFromTxt(grid, filename, c);

return grid;
}
}
throw std::runtime_error("could not find GridProperties in file " + filename);
}

void dumpGridToTxt(ref_ptr<Grid3f> grid, std::string filename, double c, bool saveProp) {
std::ofstream fout(filename.c_str());
if (!fout) {
std::stringstream ss;
ss << "dump Grid3f: " << filename << " not found";
throw std::runtime_error(ss.str());
}

// store the properties in the file as header information
if (saveProp) {
fout << "# GridProperties: Type Grid3f"
<< "\t" << "origin: " << grid -> getOrigin()
<< "\t" << "gridsize: " << grid -> getNx() << " " << grid -> getNy() << " " << grid -> getNz()
<< "\t" << "spacing: " << grid -> getSpacing ()
<< "\t" << "reflective: " << grid -> isReflective()
<< "\t" << "clipVolume: " << grid -> getClipVolume()
<< "\t" << "interpolation: " << grid -> getInterpolationTypeName() << "\n";
}

for (int ix = 0; ix < grid->getNx(); ix++) {
for (int iy = 0; iy < grid->getNy(); iy++) {
for (int iz = 0; iz < grid->getNz(); iz++) {
Expand All @@ -303,13 +451,25 @@ void dumpGridToTxt(ref_ptr<Grid3f> grid, std::string filename, double c) {
fout.close();
}

void dumpGridToTxt(ref_ptr<Grid1f> grid, std::string filename, double c) {
void dumpGridToTxt(ref_ptr<Grid1f> grid, std::string filename, double c, bool saveProp) {
std::ofstream fout(filename.c_str());
if (!fout) {
std::stringstream ss;
ss << "dump Grid1f: " << filename << " not found";
throw std::runtime_error(ss.str());
}

// save properties as header information
if (saveProp) {
fout << "# GridProperties: Type Grid1f"
<< "\t" << "origin: " << grid -> getOrigin()
<< "\t" << "gridsize: " << grid -> getNx() << " " << grid -> getNy() << " " << grid -> getNz()
<< "\t" << "spacing: " << grid -> getSpacing ()
<< "\t" << "reflective: " << grid -> isReflective()
<< "\t" << "clipVolume: " << grid -> getClipVolume()
<< "\t" << "interpolation: " << grid -> getInterpolationTypeName() << "\n";
}

for (int ix = 0; ix < grid->getNx(); ix++) {
for (int iy = 0; iy < grid->getNy(); iy++) {
for (int iz = 0; iz < grid->getNz(); iz++) {
Expand Down
Loading
Loading