From df08d3fc3880f8442ed7efa8c5a3b3e52a2872bb Mon Sep 17 00:00:00 2001 From: Ben Haller Date: Mon, 6 Nov 2023 19:31:22 -0500 Subject: [PATCH] implement 2D bicubic interpolation with periodic boundaries --- VERSIONS | 1 + core/spatial_map.cpp | 90 +++++++++++++++++++++++++++++++++----------- 2 files changed, 69 insertions(+), 22 deletions(-) diff --git a/VERSIONS b/VERSIONS index 775e295c..06c3dcb5 100644 --- a/VERSIONS +++ b/VERSIONS @@ -163,6 +163,7 @@ development head (in the master branch): add InteractionType method setConstraints(), generalizing sex-specificity with a whole raft of other possible constraints on receivers and/or exerters; interally, now keep two k-d trees, one constrained, one not fixed a bug in passing: drawByStrength() and nearestInteractingNeighbors() were not applying sex-specificity to receivers, for the returnDict=T case (never released publicly) fixed a bug in passing: partial sex-specificity (receivers but not exerters, or exerters but not receivers) was not applied flexibly in multispecies models where one species is sexual and one is not + implement 2D bicubic interpolation with periodic boundaries for the interpolate() method version 4.0.1 (Eidos version 3.0.1): diff --git a/core/spatial_map.cpp b/core/spatial_map.cpp index 276520f7..23f5a97a 100644 --- a/core/spatial_map.cpp +++ b/core/spatial_map.cpp @@ -1712,49 +1712,94 @@ EidosValue_SP SpatialMap::ExecuteMethod_interpolate(EidosGlobalStringID p_method case 2: { // bicubic interpolation - // FIXME the GSL does not support 2D periodic boundaries, but we could do it, I think, - // by replicating the map data into a 3x3 repeated pattern, doing the interpolation, and - // then taking the middle. Not unreasonable, but I will wait until someone requests it. - // When this is done, note that one of the periodic edges in each dimension must be omitted, - // since the left/right and top/bottom are (or should be) duplicates. - if (periodic) - EIDOS_TERMINATION << "ERROR (SpatialMap::ExecuteMethod_interpolate): bicubic interpolation does not currently support periodic boundaries; please open a feature request if you need this." << EidosTerminate(nullptr); - if ((grid_size_[0] < 4) || (grid_size_[1] < 4)) EIDOS_TERMINATION << "ERROR (SpatialMap::ExecuteMethod_interpolate): bicubic interpolation requires a starting map with a grid size at least 4x4." << EidosTerminate(nullptr); + // the periodic boundaries case is similar to the non-periodic case, except that we perform + // the (non-periodic) interpolation on an expanded grid with margins of 20 grid points on + // all sides (20 comes from here: https://stackoverflow.com/a/25106574/2752221 as a suggested + // margin to produce acceptably small error). Then we use the central part of that grid, + // minus the margins. + + int64_t margin = periodic ? 20 : 0; + int64_t gs0_with_margins = grid_size_[0] + margin * 2; + int64_t gs1_with_margins = grid_size_[1] + margin * 2; + + // dim_a and dim_b are the dimensions of the final grid we want, without margins; new_values is the final values int64_t dim_a = (factor * (grid_size_[0] - 1)) + 1, dim_b = (factor * (grid_size_[1] - 1)) + 1; double *new_values = (double *)malloc(dim_a * dim_b * sizeof(double)); - double *x = (double *)malloc(grid_size_[0] * sizeof(double)); - double *y = (double *)malloc(grid_size_[1] * sizeof(double)); - double *z = (double *)malloc(grid_size_[0] * grid_size_[1] * sizeof(double)); + + // x and y are the coordinates of the grid with margins; z is the original values to interpolate, with margins + double *x = (double *)malloc(gs0_with_margins * sizeof(double)); + double *y = (double *)malloc(gs1_with_margins * sizeof(double)); + double *z = (double *)malloc(gs0_with_margins * gs1_with_margins * sizeof(double)); if (!new_values || !x || !y || !z) EIDOS_TERMINATION << "ERROR (SpatialMap::ExecuteMethod_interpolate): allocation failed; you may need to raise the memory limit for SLiM." << EidosTerminate(nullptr); - // set up coordinates on our grid, not in user coordinates, for simplicity - for (int i = 0; i < grid_size_[0]; ++i) + // set up coordinates on our grid with margins, not in user coordinates, for simplicity + for (int i = 0; i < gs0_with_margins; ++i) x[i] = i; - for (int i = 0; i < grid_size_[1]; ++i) + for (int i = 0; i < gs1_with_margins; ++i) y[i] = i; const gsl_interp2d_type *T = gsl_interp2d_bicubic; - gsl_spline2d *spline = gsl_spline2d_alloc(T, grid_size_[0], grid_size_[1]); + gsl_spline2d *spline = gsl_spline2d_alloc(T, gs0_with_margins, gs1_with_margins); gsl_interp_accel *xacc = gsl_interp_accel_alloc(); gsl_interp_accel *yacc = gsl_interp_accel_alloc(); double *new_values_ptr = new_values; double scale = 1.0 / factor; - for (int64_t b = 0; b < grid_size_[1]; ++b) - for (int64_t a = 0; a < grid_size_[0]; ++a) - gsl_spline2d_set(spline, z, a, b, values_[a + b * grid_size_[0]]); + if (!periodic) + { + // in the non-periodic case, there are no margins so we can use our grid values directly + // the periodic case with margin==0 reduces to this, so this is just optimization + for (int64_t b = 0; b < gs1_with_margins; ++b) + for (int64_t a = 0; a < gs0_with_margins; ++a) + gsl_spline2d_set(spline, z, a, b, values_[a + b * grid_size_[0]]); + } + else + { + // in the periodic case, we have to add the margins, so there is some futzing around + // note that we repeat (grid_size_[0] - 1) and (grid_size_[1] - 1) elements, because + // the last column/row ought to be duplicates of the first column/row (we don't check) + // the "+ (grid_size_[1] - 1) * 10)" term is to make things positive so % is not wonky + for (int64_t b = 0; b < gs1_with_margins; ++b) + { + int64_t original_grid_b = (b - margin + (grid_size_[1] - 1) * 10) % (grid_size_[1] - 1); + + for (int64_t a = 0; a < gs0_with_margins; ++a) + { + int64_t original_grid_a = (a - margin + (grid_size_[0] - 1) * 10) % (grid_size_[0] - 1); + + gsl_spline2d_set(spline, z, a, b, values_[original_grid_a + original_grid_b * grid_size_[0]]); + } + } + } - gsl_spline2d_init(spline, x, y, z, grid_size_[0], grid_size_[1]); + gsl_spline2d_init(spline, x, y, z, gs0_with_margins, gs1_with_margins); // FIXME: TO BE PARALLELIZED - for (int64_t b = 0; b < dim_b; ++b) - for (int64_t a = 0; a < dim_a; ++a) - *(new_values_ptr++) = gsl_spline2d_eval(spline, a * scale, b * scale, xacc, yacc); + if (!periodic) + { + // in the non-periodic case, there are no margins so we can extract grid values directly + // the periodic case with margin==0 reduces to this, so this is just optimization + for (int64_t b = 0; b < dim_b; ++b) + for (int64_t a = 0; a < dim_a; ++a) + *(new_values_ptr++) = gsl_spline2d_eval(spline, a * scale, b * scale, xacc, yacc); + } + else + { + // in the periodic case, we want to extract grid values from the central area, within the margins + // recall that (factor-1) values are inserted between each grid value, so for interpolate() with + // b==0 with margin==2 and factor==3, we want to start at 6: M**M**X, X is at position 6; when + // non-periodic, margin==0 and so offset==0. + int64_t offset = margin * factor; + + for (int64_t b = 0; b < dim_b; ++b) + for (int64_t a = 0; a < dim_a; ++a) + *(new_values_ptr++) = gsl_spline2d_eval(spline, (a + offset) * scale, (b + offset) * scale, xacc, yacc); + } gsl_spline2d_free(spline); gsl_interp_accel_free(xacc); @@ -1765,6 +1810,7 @@ EidosValue_SP SpatialMap::ExecuteMethod_interpolate(EidosGlobalStringID p_method int64_t dims[2] = {dim_a, dim_b}; TakeOverMallocedValues(new_values, 2, dims); // takes new_values from us + break; } case 3: