Warning, /epic-analysis/src/interp/Interpolate.ipp is written in an unsupported language. File is not indexed.
0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023 Duane Byer
0003
0004 #ifndef Interpolate_IPP_
0005 #define Interpolate_IPP_
0006
0007 #include "Interpolate.h"
0008
0009 #include <cmath>
0010 #include <cstddef>
0011 #include <limits>
0012 #include <stdexcept>
0013 #include <string>
0014 #include <type_traits>
0015
0016 // Version of `modf` with correct sign behaviour.
0017 template<typename T>
0018 T modf(T x, T* iptr) {
0019 if (!(x < 0.)) {
0020 return std::modf(x, iptr);
0021 } else {
0022 T frac = std::modf(x, iptr);
0023 frac += 1;
0024 *iptr -= 1;
0025 return frac;
0026 }
0027 }
0028
0029 template<typename T, std::size_t N>
0030 GridView<T, N>::GridView(T const* data, Coords<T, N> coords) :
0031 _data(data),
0032 _coords(coords) {
0033 _count_total = 1;
0034 for (std::size_t idx = 0; idx < N; ++idx) {
0035 _count_total *= _coords[idx].size();
0036 if (_coords[idx].size() <= 1) {
0037 throw SingularDimensionError(idx);
0038 }
0039 }
0040 }
0041
0042 template<typename T, std::size_t N>
0043 GridView<T, N - 1> GridView<T, N>::operator[](std::size_t idx) const {
0044 std::size_t width = _count_total / count(0);
0045 std::array<std::vector<T>, N - 1> coords_red;
0046 std::copy(_coords.begin() + 1, _coords.end(), coords_red.begin());
0047 return GridView<T, N - 1>(_data + idx * width, coords_red);
0048 }
0049
0050 template<typename T, std::size_t N>
0051 T const& GridView<T, N>::operator[](CellIndex<N> idx) const {
0052 CellIndex<N - 1> idx_red;
0053 std::copy(idx.begin() + 1, idx.end(), idx_red.begin());
0054 return (*this)[idx[0]][idx_red];
0055 }
0056
0057 template<typename T, std::size_t N>
0058 Grid<T, N>::Grid(T const* data, Coords<T, N> coords) :
0059 _data(),
0060 _coords(coords) {
0061 _count_total = 1;
0062 for (std::size_t idx = 0; idx < N; ++idx) {
0063 _count_total *= _coords[idx].size();
0064 if (_coords[idx].size() <= 1) {
0065 throw SingularDimensionError(idx);
0066 }
0067 }
0068 _data = std::vector<T>(data, data + _count_total);
0069 }
0070
0071 template<typename T, std::size_t N>
0072 T GridView<T, N>::space_from_grid(std::size_t dim, std::size_t idx) const {
0073 return _coords[dim][idx];
0074 }
0075 template<typename T, std::size_t N>
0076 std::size_t GridView<T, N>::grid_from_space(std::size_t dim, T x) const {
0077 T frac;
0078 return grid_from_space(dim, x, &frac);
0079 }
0080 template<typename T, std::size_t N>
0081 std::size_t GridView<T, N>::grid_from_space(std::size_t dim, T x, T* frac) const {
0082 std::size_t lower = 0;
0083 std::size_t upper = _coords[dim].size() - 1;
0084 T lower_x = _coords[dim][lower];
0085 T upper_x = _coords[dim][upper];
0086 while (lower + 1 < upper) {
0087 T diff = (x - lower_x) / (upper_x - lower_x);
0088 if (!(diff >= 0.) || !(diff <= 1.)) {
0089 return std::numeric_limits<std::size_t>::max();
0090 }
0091 T idx_float;
0092 modf(diff * (upper - lower - 1), &idx_float);
0093 std::size_t idx = static_cast<std::size_t>(idx_float) + lower;
0094 T lower_x_new = _coords[dim][idx];
0095 T upper_x_new = _coords[dim][idx + 1];
0096 if (x >= lower_x_new && x < upper_x_new) {
0097 lower = idx;
0098 upper = idx + 1;
0099 lower_x = lower_x_new;
0100 upper_x = upper_x_new;
0101 } else if (x >= upper_x_new) {
0102 lower = idx + 1;
0103 lower_x = upper_x_new;
0104 } else if (x < lower_x_new) {
0105 upper = idx;
0106 upper_x = lower_x_new;
0107 }
0108 }
0109 *frac = x - lower_x;
0110 return lower;
0111 }
0112
0113 template<typename T, std::size_t N>
0114 Point<T, N> GridView<T, N>::space_from_grid(CellIndex<N> idx) const {
0115 Point<T, N> result;
0116 for (std::size_t dim = 0; dim < N; ++dim) {
0117 result[dim] = space_from_grid(dim, idx[dim]);
0118 }
0119 return result;
0120 }
0121 template<typename T, std::size_t N>
0122 CellIndex<N> GridView<T, N>::grid_from_space(Point<T, N> x) const {
0123 CellIndex<N> result;
0124 for (std::size_t dim = 0; dim < N; ++dim) {
0125 result[dim] = grid_from_space(dim, x[dim]);
0126 }
0127 return result;
0128 }
0129 template<typename T, std::size_t N>
0130 CellIndex<N> GridView<T, N>::grid_from_space(Point<T, N> x, Point<T, N>* frac) const {
0131 CellIndex<N> result;
0132 for (std::size_t dim = 0; dim < N; ++dim) {
0133 result[dim] = grid_from_space(dim, x[dim], &frac[dim]);
0134 }
0135 return result;
0136 }
0137
0138 template<typename T, std::size_t N>
0139 T LinearView<T, N>::operator()(Point<T, N> x) const {
0140 T frac;
0141 std::size_t idx = _grid.grid_from_space(0, x[0], &frac);
0142 Point<T, N - 1> x_next;
0143 std::copy(x.begin() + 1, x.end(), x_next.begin());
0144 if (idx >= 0 && idx < _grid.count(0) - 1) {
0145 T linear_1 = LinearView<T, N - 1>(_grid.reduce(idx))(x_next);
0146 T linear_2 = LinearView<T, N - 1>(_grid.reduce(idx + 1))(x_next);
0147 T delta = _grid.space_from_grid(0, idx + 1) - _grid.space_from_grid(0, idx);
0148 return linear(linear_1, linear_2, delta, frac);
0149 } else {
0150 return std::numeric_limits<T>::quiet_NaN();
0151 }
0152 }
0153
0154 template<typename T, std::size_t N>
0155 T CubicView<T, N>::operator()(Point<T, N> x) const {
0156 T frac;
0157 std::size_t idx = _grid.grid_from_space(0, x[0], &frac);
0158 Point<T, N - 1> x_next;
0159 std::copy(x.begin() + 1, x.end(), x_next.begin());
0160 if (idx >= 0 && idx < _grid.count(0) - 1) {
0161 T cubic_0 = idx == 0 ? 0. :
0162 CubicView<T, N - 1>(_grid.reduce(idx - 1))(x_next);
0163 T cubic_1 = CubicView<T, N - 1>(_grid.reduce(idx))(x_next);
0164 T cubic_2 = CubicView<T, N - 1>(_grid.reduce(idx + 1))(x_next);
0165 T cubic_3 = idx >= _grid.count(0) - 2 ? 0. :
0166 CubicView<T, N - 1>(_grid.reduce(idx + 2))(x_next);
0167 T delta_0 = idx == 0 ? std::numeric_limits<T>::infinity() :
0168 _grid.space_from_grid(0, idx) - _grid.space_from_grid(0, idx - 1);
0169 T delta_1 = _grid.space_from_grid(0, idx + 1) - _grid.space_from_grid(0, idx);
0170 T delta_2 = idx >= _grid.count(0) - 2 ? std::numeric_limits<T>::infinity() :
0171 _grid.space_from_grid(0, idx + 2) - _grid.space_from_grid(0, idx + 1);
0172 return cubic(cubic_0, cubic_1, cubic_2, cubic_3, delta_0, delta_1, delta_2, frac);
0173 } else {
0174 return std::numeric_limits<T>::quiet_NaN();
0175 }
0176 }
0177
0178 template<typename T, std::size_t N, std::size_t K>
0179 inline std::array<Grid<T, N>, K> read_grids(
0180 std::vector<std::array<T, N + K> > const& raw_data) {
0181 std::vector<Point<T, N> > grid_points(raw_data.size());
0182 std::vector<std::array<T, K> > data(raw_data.size());
0183 for (std::size_t row_idx = 0; row_idx < raw_data.size(); ++row_idx) {
0184 for (std::size_t dim_idx = 0; dim_idx < N; ++dim_idx) {
0185 grid_points[row_idx][dim_idx] = raw_data[row_idx][dim_idx];
0186 }
0187 for (std::size_t col_idx = 0; col_idx < K; ++col_idx) {
0188 data[row_idx][col_idx] = raw_data[row_idx][N + col_idx];
0189 }
0190 }
0191
0192 // If there are not enough grid points to form one unit cell, then give up.
0193 if (grid_points.size() < (1 << N)) {
0194 throw NotEnoughPointsError(grid_points.size(), 1 << N);
0195 }
0196
0197 std::size_t count_total = 1;
0198 std::array<bool, N> filled_cols = { false };
0199 std::array<std::size_t, N> dim_permute_map;
0200 CellIndex<N> count;
0201 std::array<std::vector<T>, N> coords;
0202 for (std::size_t dim_idx = 0; dim_idx < N; ++dim_idx) {
0203 std::size_t sub_count = grid_points.size() / count_total;
0204 if (sub_count * count_total != grid_points.size()) {
0205 throw NotEnoughPointsError(
0206 grid_points.size(),
0207 (sub_count + 1) * count_total);
0208 }
0209 // Find the column that steps by `count_total`.
0210 std::size_t step_dim_idx = 0;
0211 while (true) {
0212 if (step_dim_idx == N) {
0213 throw NotEnoughPointsError(
0214 grid_points.size(),
0215 2 * grid_points.size());
0216 }
0217 if (filled_cols[step_dim_idx]) {
0218 step_dim_idx += 1;
0219 continue;
0220 }
0221 T initial_point = grid_points[0][step_dim_idx];
0222 T final_point = grid_points[count_total][step_dim_idx];
0223 if (initial_point == final_point) {
0224 step_dim_idx += 1;
0225 continue;
0226 }
0227 bool same = true;
0228 for (std::size_t row_idx = 0; row_idx < count_total; ++row_idx) {
0229 if (grid_points[row_idx][step_dim_idx] != initial_point) {
0230 same = false;
0231 break;
0232 }
0233 }
0234 if (!same) {
0235 step_dim_idx += 1;
0236 continue;
0237 }
0238 filled_cols[step_dim_idx] = true;
0239 dim_permute_map[dim_idx] = step_dim_idx;
0240 break;
0241 }
0242
0243 // Use the column to fill a vector of the grid spacings.
0244 std::vector<T> next_coords;
0245 for (std::size_t group_idx = 0; group_idx < sub_count; ++group_idx) {
0246 T next = grid_points[group_idx * count_total][step_dim_idx];
0247 if (next_coords.empty() || next != next_coords[0]) {
0248 next_coords.push_back(next);
0249 } else {
0250 break;
0251 }
0252 }
0253 std::size_t next_count = next_coords.size();
0254 count[step_dim_idx] = next_count;
0255 if (next_count <= 1) {
0256 throw SingularDimensionError(step_dim_idx);
0257 }
0258
0259 // Check that the data in the column is valid.
0260 for (std::size_t row_idx = 0; row_idx < grid_points.size(); ++row_idx) {
0261 T next = grid_points[row_idx][step_dim_idx];
0262 std::size_t group_idx = row_idx / count_total;
0263 if (next != next_coords[group_idx % next_count]) {
0264 throw UnexpectedGridPointError(row_idx);
0265 }
0266 }
0267
0268 // Check that the grid is increasing.
0269 for (std::size_t idx = 1; idx < next_coords.size(); ++idx) {
0270 if (!(next_coords[idx] > next_coords[idx - 1])) {
0271 throw InvalidGridPlanesError();
0272 }
0273 }
0274 coords[step_dim_idx] = next_coords;
0275
0276 // Increase count for next round.
0277 count_total *= next_count;
0278 }
0279
0280 // Use the grid information to reorder the data.
0281 std::array<std::vector<T>, K> data_transposed;
0282 std::array<std::size_t, N> count_totals = { count_total / count[0] };
0283 for (std::size_t dim_idx = 1; dim_idx < N; ++dim_idx) {
0284 count_totals[dim_idx] = count_totals[dim_idx - 1] / count[dim_idx];
0285 }
0286 for (std::size_t col_idx = 0; col_idx < K; ++col_idx) {
0287 data_transposed[col_idx].resize(data.size());
0288 }
0289 for (std::size_t row_idx = 0; row_idx < data.size(); ++row_idx) {
0290 std::size_t new_row_idx = 0;
0291 std::size_t new_count = 1;
0292 for (std::size_t dim_idx = 0; dim_idx < N; ++dim_idx) {
0293 std::size_t new_dim_idx = dim_permute_map[dim_idx];
0294 std::size_t rel_idx = (row_idx / new_count) % count[new_dim_idx];
0295 std::size_t old_count = count_totals[new_dim_idx];
0296 new_row_idx += rel_idx * old_count;
0297 new_count *= count[new_dim_idx];
0298 }
0299 for (std::size_t col_idx = 0; col_idx < K; ++col_idx) {
0300 data_transposed[col_idx][new_row_idx] = data[row_idx][col_idx];
0301 }
0302 }
0303
0304 std::array<Grid<T, N>, K> grids;
0305 for (std::size_t col_idx = 0; col_idx < K; ++col_idx) {
0306 grids[col_idx] = Grid<T, N>(data_transposed[col_idx].data(), coords);
0307 }
0308 return grids;
0309 }
0310
0311 inline NotEnoughPointsError::NotEnoughPointsError(
0312 std::size_t points,
0313 std::size_t expected_points) :
0314 std::runtime_error(
0315 "Not enough points to construct the grid (found "
0316 + std::to_string(points) + ", expected "
0317 + std::to_string(expected_points) + ")"),
0318 points(points),
0319 expected_points(expected_points) { }
0320
0321 inline SingularDimensionError::SingularDimensionError(std::size_t dim) :
0322 std::runtime_error("Grid is singular in dimension " + std::to_string(dim)),
0323 dim(dim) { }
0324
0325 inline InvalidGridPlanesError::InvalidGridPlanesError() :
0326 std::runtime_error("Lower grid bound must be smaller than upper bound") { }
0327
0328 inline InvalidSpacingError::InvalidSpacingError() :
0329 std::runtime_error("Grid must be spaced uniformly") { }
0330
0331 inline UnexpectedGridPointError::UnexpectedGridPointError(
0332 std::size_t line_number) :
0333 std::runtime_error(
0334 "Unexpected grid point at line " + std::to_string(line_number)),
0335 line_number(line_number) { }
0336
0337 #endif
0338
