Skip to content

Commit

Permalink
FEAT: Support set_interpolation_coefs and get_interpolation_values in…
Browse files Browse the repository at this point in the history 
… Issue AetherModel#99
  • Loading branch information
@oscarhua
oscarhua committed May 31, 2023
1 parent 52878c8 commit 6214038
Show file tree
Hide file tree
Showing 4 changed files with 415 additions and 13 deletions.
73 changes: 65 additions & 8 deletions include/grid.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -183,14 +183,41 @@ class Grid {
bool send_one_face(int64_t iFace);
bool receive_one_face(int64_t iFace);

// interpolation functions
// Estimate the value of the point at (lon_in, lat_in, alt_in)
/**
* \brief Set the interpolation coefficients
* \param Lons The longitude of points
* \param Lats The latitude of points
* \param Alts The altitude of points
* \pre This instance is an geo grid
* \pre Lons, Lats and Alts have the same size
* \return true if the function succeeds, false otherwise.
*/
bool set_interpolation_coefs(const std::vector<precision_t> &Lons,
const std::vector<precision_t> &Lats,
const std::vector<precision_t> &Alts);
/**
* \brief Create a map of geographic locations to data and do the interpolation
* \param data The value at the positions of geoLon, geoLat, and geoAlt
* \pre The size of the data should be the same as the geoLat/Lon/Alt_scgc
* \return A vector of estimated value at the points set by the last
* set_interpolation_coefs function call if the function succeeds,
* an empty vector otherwise.
*/
std::vector<precision_t> get_interpolation_values(const arma_cube &data) const;

/**
* \deprecated !!!Use set_interpolation_coefs and get_interpolation_values instead!!!
* \brief Interpolate the value of the point at (lon_in, lat_in, alt_in)
*/
precision_t interp_linear(const arma_cube &data,
const precision_t lon,
const precision_t lat,
const precision_t alt);
// the position of a vector of points is specified by
// three vectors of lon, lat and alt
/**
* \deprecated !!!Use set_interpolation_coefs and get_interpolation_values instead!!!
* \brief Interpolation. The position of a vector of points is specified by
* three vectors of lon, lat and alt
*/
std::vector<precision_t> interp_linear(const arma_cube &data,
const std::vector<precision_t> &Lons,
const std::vector<precision_t> &Lats,
Expand Down Expand Up @@ -247,15 +274,32 @@ class Grid {
bool col_max_exclusive;
};

// The index and coefficient used for interpolation
// Each point is processed by the function set_interpolation_coefs and stored
// in the form of this structure.
// If the point is out of the grid, in_grid = false and all other members are undefined
struct interp_coef_t {
// The point is inside the cube of [iRow, iRow+1], [iCol, iCol+1], [iAlt, iAlt+1]
uint64_t iRow;
uint64_t iCol;
uint64_t iAlt;
// The coefficients along row, column and altitude
precision_t rRow;
precision_t rCol;
precision_t rAlt;
// Whether the point is within this grid or not
bool in_grid;
};

// Return the index of the last element that has altitude smaller than or euqal to the input
uint64_t search_altitude(const precision_t alt_in);
uint64_t search_altitude(const precision_t alt_in) const;

// Calculate the range of a spherical grid
void get_sphere_grid_range(struct sphere_range &sr);
void get_sphere_grid_range(struct sphere_range &sr) const;
// Calculate the range of a cubesphere grid
void get_cubesphere_grid_range(struct cubesphere_range &cr);
void get_cubesphere_grid_range(struct cubesphere_range &cr) const;

// Helper function for interpolation, so that grid range is only
// Helper function for interp_linear, so that grid range is only
// calculated once no matter how many points
precision_t interp_sphere_linear_helper(const arma_cube &data,
const sphere_range &sr,
Expand All @@ -267,6 +311,19 @@ class Grid {
const precision_t lon_in,
const precision_t lat_in,
const precision_t alt_in);

// Helper function for set_interpolation_coefs
void set_interp_coef_sphere(const sphere_range &sr,
const precision_t lon_in,
const precision_t lat_in,
const precision_t alt_in);
void set_interp_coef_cubesphere(const cubesphere_range &cr,
const precision_t lon_in,
const precision_t lat_in,
const precision_t alt_in);

// Processed interpolation coefficients
std::vector<struct interp_coef_t> interp_coefs;
};

#endif // INCLUDE_GRID_H_
7 changes: 6 additions & 1 deletion include/solvers.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,12 @@ arma_cube calc_gradient_lat(arma_cube value, Grid grid);
arma_cube calc_gradient_alt(arma_cube value, Grid grid);
std::vector<arma_cube> calc_gradient_vector(arma_cube value_scgc, Grid grid);

// interpolation in 3D
// interpolation in 1D
precision_t linear_interpolation(const precision_t y0,
const precision_t y1,
const precision_t ratio);

// interpolation in 3D, data should be a cube of size 2-2-2
precision_t interpolate_unit_cube(const arma_cube &data,
const precision_t xRatio,
const precision_t yRatio,
Expand Down
137 changes: 137 additions & 0 deletions src/main/main_test_interpolation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,6 +96,7 @@ int main() {
nGeoGhosts);
DidWork = gGrid.init_geo_grid(quadtree, planet, input, report);

// !!!This part tests Grid::interp_linear, which is desprecated!!!
{
// What the test do is:
// Given the position of an unknown point, estimate its position :)
Expand Down Expand Up @@ -213,6 +214,142 @@ int main() {
}
}

// This part tests Grid::set_interpolation_coefs and get_interpolation_values
{
// What the test do is:
// Given the position of an unknown point, estimate its position :)
// Thus the output should be equal to the input

// cPI = 3.14159, 2*cPI = 6.28318, cPI/2 = 1.57079
// minAlt = 100,000, maxAlt = 100,000 + 5,000 * 50 = 350,000
// Lon within [0, 2*cPI) = [0, 6.28318)
// Lat within [-cPI/2, cPI/2] = [-1.57079, 1.57079]
// Alt within [100000, 350000]

std::cout << "iProc = " << iProc;

// Edge condition for spherical grid
{
std::cout << "\n\n\nEdge 1\n\n\n";
std::vector<precision_t> lon_in({0, 0, 0, 0, 2*cPI - cSmall, 2*cPI - cSmall, 2*cPI - cSmall, 2*cPI - cSmall});
std::vector<precision_t> lat_in({-cPI/2, -cPI/2, cPI/2, cPI/2, -cPI/2, -cPI/2, cPI/2, cPI/2});
std::vector<precision_t> alt_in({100000, 350000, 100000, 350000, 100000, 350000, 100000, 350000});

gGrid.set_interpolation_coefs(lon_in, lat_in, alt_in);

std::vector<precision_t> ansLon = gGrid.get_interpolation_values(gGrid.geoLon_scgc);
std::vector<precision_t> ansLat = gGrid.get_interpolation_values(gGrid.geoLat_scgc);
std::vector<precision_t> ansAlt = gGrid.get_interpolation_values(gGrid.geoAlt_scgc);

for (int64_t i = 0; i < lon_in.size(); ++i) {
std::cout << "Estimation of (" << lon_in[i] << ", " << lat_in[i] << ", " << alt_in[i] << ") is "
<< "(" << ansLon[i] << ", " << ansLat[i] << ", " << ansAlt[i] << ")\n";
}
is_correct(lon_in, ansLon, lat_in, ansLat, alt_in, ansAlt);
}

// Edge condition for cubeshpere grid
{
std::cout << "\n\n\nEdge 2\n\n\n";
std::vector<precision_t> lon_in({0, 0, cPI / 2, cPI / 2, cPI, cPI, 3 * cPI / 2, 3 * cPI / 2});
std::vector<precision_t> lat_in({-0.61547970867, 0.61547970867, -0.61547970867, 0.61547970867, -0.61547970867, 0.61547970867, -0.61547970867, 0.61547970867});
std::vector<precision_t> alt_in({100000, 350000, 100000, 350000, 100000, 350000, 100000, 350000});

gGrid.set_interpolation_coefs(lon_in, lat_in, alt_in);

std::vector<precision_t> ansLon = gGrid.get_interpolation_values(gGrid.geoLon_scgc);
std::vector<precision_t> ansLat = gGrid.get_interpolation_values(gGrid.geoLat_scgc);
std::vector<precision_t> ansAlt = gGrid.get_interpolation_values(gGrid.geoAlt_scgc);

for (int64_t i = 0; i < lon_in.size(); ++i) {
std::cout << "Estimation of (" << lon_in[i] << ", " << lat_in[i] << ", " << alt_in[i] << ") is "
<< "(" << ansLon[i] << ", " << ansLat[i] << ", " << ansAlt[i] << ")\n";
}
is_correct(lon_in, ansLon, lat_in, ansLat, alt_in, ansAlt);
}

// Error handling
{
std::cout << "\n\n\nError handling 1\n\n\n";
std::vector<precision_t> lon_in({-cSmall, -cSmall, -cSmall, -cSmall, 2*cPI, 2*cPI, 2*cPI, 2*cPI});
std::vector<precision_t> lat_in({-cPI/2, -cPI/2, cPI/2, cPI/2, -cPI/2, -cPI/2, cPI/2, cPI/2});
std::vector<precision_t> alt_in({100000, 350000, 100000, 350000, 100000, 350000, 100000, 350000});

gGrid.set_interpolation_coefs(lon_in, lat_in, alt_in);

std::vector<precision_t> ansLon = gGrid.get_interpolation_values(gGrid.geoLon_scgc);
std::vector<precision_t> ansLat = gGrid.get_interpolation_values(gGrid.geoLat_scgc);
std::vector<precision_t> ansAlt = gGrid.get_interpolation_values(gGrid.geoAlt_scgc);

for (int64_t i = 0; i < lon_in.size(); ++i) {
std::cout << "Estimation of (" << lon_in[i] << ", " << lat_in[i] << ", " << alt_in[i] << ") is "
<< "(" << ansLon[i] << ", " << ansLat[i] << ", " << ansAlt[i] << ")\n";
}
is_correct(lon_in, ansLon, lat_in, ansLat, alt_in, ansAlt);
}


{
std::cout << "\n\n\nError handling 2\n\n\n";
std::vector<precision_t> lon_in({0, 0, 0, 0, 2*cPI - cSmall, 2*cPI - cSmall, 2*cPI - cSmall, 2*cPI - cSmall});
std::vector<precision_t> lat_in({-cPI/2 - cSmall, -cPI/2 - cSmall, cPI/2 + cSmall, cPI/2 + cSmall, -cPI/2 - cSmall, -cPI/2 - cSmall, cPI/2 + cSmall, cPI/2 + cSmall});
std::vector<precision_t> alt_in({100000, 350000, 100000, 350000, 100000, 350000, 100000, 350000});

gGrid.set_interpolation_coefs(lon_in, lat_in, alt_in);

std::vector<precision_t> ansLon = gGrid.get_interpolation_values(gGrid.geoLon_scgc);
std::vector<precision_t> ansLat = gGrid.get_interpolation_values(gGrid.geoLat_scgc);
std::vector<precision_t> ansAlt = gGrid.get_interpolation_values(gGrid.geoAlt_scgc);

for (int64_t i = 0; i < lon_in.size(); ++i) {
std::cout << "Estimation of (" << lon_in[i] << ", " << lat_in[i] << ", " << alt_in[i] << ") is "
<< "(" << ansLon[i] << ", " << ansLat[i] << ", " << ansAlt[i] << ")\n";
}
is_correct(lon_in, ansLon, lat_in, ansLat, alt_in, ansAlt);
}


{
std::cout << "\n\n\nError handling 3\n\n\n";
std::vector<precision_t> lon_in({0, 0, 0, 0, 2*cPI - cSmall, 2*cPI - cSmall, 2*cPI - cSmall, 2*cPI - cSmall});
std::vector<precision_t> lat_in({-cPI/2, -cPI/2, cPI/2, cPI/2, -cPI/2, -cPI/2, cPI/2, cPI/2});
std::vector<precision_t> alt_in({100000 - 1, 350000 + 1, 100000 - 1, 350000 + 1, 100000 - 1, 350000 + 1, 100000 - 1, 350000 + 1});

gGrid.set_interpolation_coefs(lon_in, lat_in, alt_in);

std::vector<precision_t> ansLon = gGrid.get_interpolation_values(gGrid.geoLon_scgc);
std::vector<precision_t> ansLat = gGrid.get_interpolation_values(gGrid.geoLat_scgc);
std::vector<precision_t> ansAlt = gGrid.get_interpolation_values(gGrid.geoAlt_scgc);

for (int64_t i = 0; i < lon_in.size(); ++i) {
std::cout << "Estimation of (" << lon_in[i] << ", " << lat_in[i] << ", " << alt_in[i] << ") is "
<< "(" << ansLon[i] << ", " << ansLat[i] << ", " << ansAlt[i] << ")\n";
}
is_correct(lon_in, ansLon, lat_in, ansLat, alt_in, ansAlt);
}

// Normal condition
{
std::cout << "\n\n\nNormal condition\n\n\n";
std::vector<precision_t> lon_in({0.1, 0.2, 0.3, 0.4, 2*cPI - 0.5, 2*cPI - 0.6, 2*cPI - 0.7, cPI, cPI});
std::vector<precision_t> lat_in({-cPI/2 + 0.5, -cPI/2 + 0.6, 0, 0.1, cPI/2 - 0.5, cPI/2 - 0.6, -0.1, -0.2, 0});
std::vector<precision_t> alt_in({100000, 150000, 200000, 250000, 300000, 350000, 120000, 225000, 225000});

gGrid.set_interpolation_coefs(lon_in, lat_in, alt_in);

std::vector<precision_t> ansLon = gGrid.get_interpolation_values(gGrid.geoLon_scgc);
std::vector<precision_t> ansLat = gGrid.get_interpolation_values(gGrid.geoLat_scgc);
std::vector<precision_t> ansAlt = gGrid.get_interpolation_values(gGrid.geoAlt_scgc);

for (int64_t i = 0; i < lon_in.size(); ++i) {
std::cout << "Estimation of (" << lon_in[i] << ", " << lat_in[i] << ", " << alt_in[i] << ") is "
<< "(" << ansLon[i] << ", " << ansLat[i] << ", " << ansAlt[i] << ")\n";
}
std::cout << "\n\n\n";
is_correct(lon_in, ansLon, lat_in, ansLat, alt_in, ansAlt);
}
}

report.exit(function);
report.times();

Expand Down
Loading

0 comments on commit 6214038

Please sign in to comment.