diff --git a/include/grid.h b/include/grid.h index 0846620a..3290697e 100644 --- a/include/grid.h +++ b/include/grid.h @@ -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 &Lons, + const std::vector &Lats, + const std::vector &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 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 interp_linear(const arma_cube &data, const std::vector &Lons, const std::vector &Lats, @@ -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, @@ -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 interp_coefs; }; #endif // INCLUDE_GRID_H_ diff --git a/include/solvers.h b/include/solvers.h index 49fc35f1..2911faa7 100644 --- a/include/solvers.h +++ b/include/solvers.h @@ -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 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, diff --git a/src/main/main_test_interpolation.cpp b/src/main/main_test_interpolation.cpp index d5e9bf6b..8b6af903 100644 --- a/src/main/main_test_interpolation.cpp +++ b/src/main/main_test_interpolation.cpp @@ -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 :) @@ -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 lon_in({0, 0, 0, 0, 2*cPI - cSmall, 2*cPI - cSmall, 2*cPI - cSmall, 2*cPI - cSmall}); + std::vector lat_in({-cPI/2, -cPI/2, cPI/2, cPI/2, -cPI/2, -cPI/2, cPI/2, cPI/2}); + std::vector alt_in({100000, 350000, 100000, 350000, 100000, 350000, 100000, 350000}); + + gGrid.set_interpolation_coefs(lon_in, lat_in, alt_in); + + std::vector ansLon = gGrid.get_interpolation_values(gGrid.geoLon_scgc); + std::vector ansLat = gGrid.get_interpolation_values(gGrid.geoLat_scgc); + std::vector 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 lon_in({0, 0, cPI / 2, cPI / 2, cPI, cPI, 3 * cPI / 2, 3 * cPI / 2}); + std::vector lat_in({-0.61547970867, 0.61547970867, -0.61547970867, 0.61547970867, -0.61547970867, 0.61547970867, -0.61547970867, 0.61547970867}); + std::vector alt_in({100000, 350000, 100000, 350000, 100000, 350000, 100000, 350000}); + + gGrid.set_interpolation_coefs(lon_in, lat_in, alt_in); + + std::vector ansLon = gGrid.get_interpolation_values(gGrid.geoLon_scgc); + std::vector ansLat = gGrid.get_interpolation_values(gGrid.geoLat_scgc); + std::vector 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 lon_in({-cSmall, -cSmall, -cSmall, -cSmall, 2*cPI, 2*cPI, 2*cPI, 2*cPI}); + std::vector lat_in({-cPI/2, -cPI/2, cPI/2, cPI/2, -cPI/2, -cPI/2, cPI/2, cPI/2}); + std::vector alt_in({100000, 350000, 100000, 350000, 100000, 350000, 100000, 350000}); + + gGrid.set_interpolation_coefs(lon_in, lat_in, alt_in); + + std::vector ansLon = gGrid.get_interpolation_values(gGrid.geoLon_scgc); + std::vector ansLat = gGrid.get_interpolation_values(gGrid.geoLat_scgc); + std::vector 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 lon_in({0, 0, 0, 0, 2*cPI - cSmall, 2*cPI - cSmall, 2*cPI - cSmall, 2*cPI - cSmall}); + std::vector 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 alt_in({100000, 350000, 100000, 350000, 100000, 350000, 100000, 350000}); + + gGrid.set_interpolation_coefs(lon_in, lat_in, alt_in); + + std::vector ansLon = gGrid.get_interpolation_values(gGrid.geoLon_scgc); + std::vector ansLat = gGrid.get_interpolation_values(gGrid.geoLat_scgc); + std::vector 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 lon_in({0, 0, 0, 0, 2*cPI - cSmall, 2*cPI - cSmall, 2*cPI - cSmall, 2*cPI - cSmall}); + std::vector lat_in({-cPI/2, -cPI/2, cPI/2, cPI/2, -cPI/2, -cPI/2, cPI/2, cPI/2}); + std::vector 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 ansLon = gGrid.get_interpolation_values(gGrid.geoLon_scgc); + std::vector ansLat = gGrid.get_interpolation_values(gGrid.geoLat_scgc); + std::vector 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 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 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 alt_in({100000, 150000, 200000, 250000, 300000, 350000, 120000, 225000, 225000}); + + gGrid.set_interpolation_coefs(lon_in, lat_in, alt_in); + + std::vector ansLon = gGrid.get_interpolation_values(gGrid.geoLon_scgc); + std::vector ansLat = gGrid.get_interpolation_values(gGrid.geoLat_scgc); + std::vector 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(); diff --git a/src/solver_grid_interpolation.cpp b/src/solver_grid_interpolation.cpp index ad1d1a13..96d3fc73 100644 --- a/src/solver_grid_interpolation.cpp +++ b/src/solver_grid_interpolation.cpp @@ -3,7 +3,7 @@ #include "aether.h" -// Hepler varialbes / function begins. These are only used inside this cpp file and not declared in any other file +// Hepler varialbes / function begins. These are only used inside this cpp file and neither declared nor visible in any other file // The size of a 2*2*2 arma cube const arma::SizeCube unit_cube_size = arma::size(2, 2, 2); @@ -97,11 +97,16 @@ int64_t get_cube_surface_number(const arma_vec &point_in) { // Helper variables / function ends. The following are all member functions of Grid class +// Read grid.h completely before trying to understand the following. +// Currently interp_linear(both 2 overloads), interp_sphere_linear_helper +// and interp_cubesphere_linear_helper are deprecated (i.e. no longer recommended +// for use due to low efficiency, but the functionality is still maintained) + // -------------------------------------------------------------------------- // Return the index of the last element that has altitude smaller than or euqal to the input // -------------------------------------------------------------------------- -uint64_t Grid::search_altitude(const precision_t alt_in) { +uint64_t Grid::search_altitude(const precision_t alt_in) const { // Copy from std::upper_bound. Can't directly use it // mainly because geoAlt_scgc(0, 0, *) can't be formed as an iterator uint64_t first, last, len; @@ -125,7 +130,7 @@ uint64_t Grid::search_altitude(const precision_t alt_in) { // Get the range of a spherical grid // -------------------------------------------------------------------------- -void Grid::get_sphere_grid_range(struct sphere_range &sr) { +void Grid::get_sphere_grid_range(struct sphere_range &sr) const { // Retrieve the range and delta of longitude, latitude and altitude sr.lon_min = geoLon_Corner(nGCs, nGCs, nGCs); sr.lon_max = geoLon_Corner(nLons - nGCs, nLats - nGCs, nAlts - nGCs); @@ -144,7 +149,7 @@ void Grid::get_sphere_grid_range(struct sphere_range &sr) { // Get the range of a cubesphere grid // -------------------------------------------------------------------------- -void Grid::get_cubesphere_grid_range(struct cubesphere_range &cr) { +void Grid::get_cubesphere_grid_range(struct cubesphere_range &cr) const { // Get the location of the lower left corner, one step for row and one step for column arma_vec corner = sphere_to_cube(geoLon_Corner(nGCs, nGCs, nGCs), geoLat_Corner(nGCs, nGCs, nGCs)); @@ -396,3 +401,201 @@ std::vector Grid::interp_linear(const arma_cube &data, return ans; } + + +// -------------------------------------------------------------------------- +// Set interpolation coefficients helper function for spherical grid +// Almost the copy of interp_sphere_linear_helper +// -------------------------------------------------------------------------- + +void Grid::set_interp_coef_sphere(const sphere_range &sr, + const precision_t lon_in, + const precision_t lat_in, + const precision_t alt_in) { + // WARNING: IF WE ARE DEALING WITH LESS THAN THE WHOLE EARTH, THEN ALL THE POINTS WITH + // LONGITUDE = geo_grid_input.lon_max = settings["GeoGrid"]["MaxLon"] + // OR LATITUDE = geo_grid_input.lat_max = settings["GeoGrid"]["MaxLat"] + // ARE EXCLUDED. + // TO FIX IT, EACH GRID SHOULD BE ABLE TO ACCESS THE MaxLon and MaxLat + + // The structure which will be put into the interp_coefs. Initialize in_grid to be false + struct interp_coef_t coef; + coef.in_grid = false; + + // Determine whether the point is inside this grid + // Treat north pole specially because latitude is inclusive for both -cPI/2 and cPI/2 + if (lon_in < sr.lon_min || lon_in >= sr.lon_max || lat_in < sr.lat_min + || lat_in > sr.lat_max || (lat_in == sr.lat_max && sr.lat_max != cPI/2) + || alt_in < sr.alt_min || alt_in > sr.alt_max) { + interp_coefs.push_back(coef); + return; + } + + // ASSUMPTION: LONGITUDE AND LATITUDE ARE LINEARLY SPACED, nGCs >= 1 + // For the cell containing it, directly calculate its x and y index + // Find its z index using binary search + + // The number of dLon between the innermost ghost cell and the given point + coef.rRow = (lon_in - sr.lon_min) / sr.dLon + 0.5; + // Take the integer part + coef.iRow = static_cast(coef.rRow); + // Calculate the fractional part, which is the ratio for Longitude + coef.rRow -= coef.iRow; + // The actual x-axis index of the bottom-left of the cube used for interpolation + coef.iRow += nGCs - 1; + // Do the same for the Latitude + coef.rCol = (lat_in - sr.lat_min) / sr.dLat + 0.5; + coef.iCol = static_cast(coef.rCol); + coef.rCol -= coef.iCol; + coef.iCol += nGCs - 1; + + // The altitude may not be linearly spaced, so use binary search to find + // the first element smaller than or equal to the altitude of the give point + // Implemented in search_altitude + coef.iAlt = search_altitude(alt_in); + coef.rAlt = (alt_in - geoAlt_scgc(0, 0, coef.iAlt)) + / (geoAlt_scgc(0, 0, coef.iAlt + 1) - geoAlt_scgc(0, 0, coef.iAlt)); + + // Put the coefficient into the vector + coef.in_grid = true; + interp_coefs.push_back(coef); +} + +// -------------------------------------------------------------------------- +// Set interpolation coefficients helper function for cubesphere grid +// Almost the copy of interp_cubesphere_linear_helper +// -------------------------------------------------------------------------- + +void Grid::set_interp_coef_cubesphere(const cubesphere_range &cr, + const precision_t lon_in, + const precision_t lat_in, + const precision_t alt_in) { + // ASSUMPTION: THE SURFACES OF THE CUBE IS LINEARLY SPACED + // I.E. init_geo_grid.cpp:106-137 WILL NEVER BE CHANGED + + // The structure which will be put into the interp_coefs. Initialize in_grid to be false + struct interp_coef_t coef; + coef.in_grid = false; + + // Find the projection point onto the cube and its surface number + arma_vec point_in = sphere_to_cube(lon_in, lat_in); + int64_t surface_in = get_cube_surface_number(point_in); + + // Determine whether the projection point is on the surface of the grid + if (surface_in != cr.surface_number) { + interp_coefs.push_back(coef); + return; + } + + // Calculate the theoretical fractional row index and column index + precision_t row_frac_index, col_frac_index, row_in, col_in; + row_in = point_in(cr.row_direction); + col_in = point_in(cr.col_direction); + row_frac_index = (row_in - cr.row_min) / cr.drow; + col_frac_index = (col_in - cr.col_min) / cr.dcol; + + // Determine whether the projection point is out of range + int64_t row_index_max, col_index_max; + row_index_max = nLons - 2 * nGCs; + col_index_max = nLats - 2 * nGCs; + if (row_frac_index < 0 || (row_frac_index == 0 && cr.row_min_exclusive) + || col_frac_index < 0 || (col_frac_index == 0 && cr.col_min_exclusive) + || row_frac_index > row_index_max || (row_frac_index == row_index_max && cr.row_max_exclusive) + || col_frac_index > col_index_max || (col_frac_index == col_index_max && cr.col_max_exclusive) + || alt_in < cr.alt_min || alt_in > cr.alt_max) { + interp_coefs.push_back(coef); + return; + } + + // Get the real integer index and the interpolation coefficient + uint64_t row_index, col_index, alt_index; + precision_t rRow, rCol, rAlt; + // Add 0.5 because the data we have is at the center of the cell rather than corner of the cell + row_frac_index += 0.5; + // Take the integer part + coef.iRow = static_cast(row_frac_index); + // Calculate the fractional part, which is the coefficient + coef.rRow = row_frac_index - coef.iRow; + // The actual index considering the ghost cells + coef.iRow += nGCs - 1; + // Do the same for the column + col_frac_index += 0.5; + coef.iCol = static_cast(col_frac_index); + coef.rCol = col_frac_index - coef.iCol; + coef.iCol += nGCs - 1; + // Use binary search to find the index for altitude + coef.iAlt = search_altitude(alt_in); + coef.rAlt = (alt_in - geoAlt_scgc(0, 0, coef.iAlt)) + / (geoAlt_scgc(0, 0, coef.iAlt + 1) - geoAlt_scgc(0, 0, coef.iAlt)); + + // Put the coefficient into the vector + coef.in_grid = true; + interp_coefs.push_back(coef); +} + +// -------------------------------------------------------------------------- +// Set the interpolation coefficients +// -------------------------------------------------------------------------- + +bool Grid::set_interpolation_coefs(const std::vector &Lons, + const std::vector &Lats, + const std::vector &Alts) { + // If this is not a geo grid, return false + if (!IsGeoGrid) { + return false; + } + + // If the size of Lons, Lats and Alts are not the same, return false + if (Lons.size() != Lats.size() || Lats.size() != Alts.size()) { + return false; + } + + // Clear the previous interpolation coefficients + interp_coefs.clear(); + + // Handle according to whether it is cubesphere or not + if (IsCubeSphereGrid) { + // Calculate the range of the grid + struct cubesphere_range cr; + get_cubesphere_grid_range(cr); + // Calculate the index and coefficients for each point + for (size_t i = 0; i < Lons.size(); ++i) { + set_interp_coef_cubesphere(cr, Lons[i], Lats[i], Alts[i]); + } + } else { + // Calculate the range of the grid + struct sphere_range sr; + get_sphere_grid_range(sr); + // Calculate the index and coefficients for each point + for (size_t i = 0; i < Lons.size(); ++i) { + set_interp_coef_sphere(sr, Lons[i], Lats[i], Alts[i]); + } + } + return true; +} + +std::vector Grid::get_interpolation_values(const arma_cube &data) const { + std::vector ans; + + // If the size of data is not the same as the size of grid, return an empty vector + if (data.n_rows != nLons || data.n_cols != nLats || data.n_slices != nAlts) { + return ans; + } + + for (auto &it : interp_coefs) { + // Do interpolation if in_grid = true. Push cNinf otherwise + if (it.in_grid) { + ans.push_back(interpolate_unit_cube( + data.subcube(it.iRow, it.iCol, it.iAlt, unit_cube_size), + it.rRow, + it.rCol, + it.rAlt + )); + // Add std::cout if needed here + // std::cout << "iProc = " << iProc << " interpolates the point successfully\n"; + } else { + ans.push_back(cNinf); + } + } + return ans; +}