– (object<SpatialMap>$)interpolate(integer$ factor, [string$ method = "linear"])
Increases the resolution of the spatial map by factor, changing the dimensions of the spatial map’s grid of values (while leaving its spatial bounds unchanged), by interpolating new values between the existing values. The parameter factor must be an integer in [2, 10001], somewhat arbitrarily. The target spatial map is returned, to allow easy chaining of operations.
For a 1D spatial map, factor-1 new values will be inserted between every pair of values in the original value grid. A factor of 2 would therefore insert one new value between each pair of existing values, thereby increasing the map’s resolution by a factor of two. Note that if the spatial map’s original grid dimension was N, the new grid dimension with a factor of k would be k(N−1)+1, not kN, because new values are inserted only between existing values. For 2D and 3D spatial maps, essentially the same process is conducted along each axis of the map’s spatiality, increasing the resolution of the map by factor in every dimension.
-If method is "linear" (the default), linear (or bilinear or trilinear, for 2D/3D maps) interpolation will be used to interpolate the values for the new grid points. Alternatively, if method is "nearest", the nearest value in the old grid will be used for new grid points; with this method, it is recommended that factor be odd, not even, to avoid artifacts due to rounding of coordinates midway between the original grid positions. If method is "cubic", cubic (or bicubic, for 2D maps) will be used; this generally produces smoother interpolation with fewer artifacts than "linear", but it is not supported for 3D maps. The choice of interpolation method used here is independent of the map’s interpolate property. Note that while the "nearest" and "linear" interpolation methods will leave the range of values in the map unchanged, "cubic" interpolation may produce interpolated values that are outside the original range of values (by design).
+If method is "linear" (the default), linear (or bilinear or trilinear, for 2D/3D maps) interpolation will be used to interpolate the values for the new grid points. Alternatively, if method is "nearest", the nearest value in the old grid will be used for new grid points; with this method, it is recommended that factor be odd, not even, to avoid artifacts due to rounding of coordinates midway between the original grid positions. If method is "cubic", cubic (or bicubic, for 2D maps) will be used; this generally produces smoother interpolation with fewer artifacts than "linear", but it is not supported for 3D maps. The choice of interpolation method used here is independent of the map’s interpolate property. Note that while the "nearest" and "linear" interpolation methods will leave the range of values in the map unchanged, "cubic" interpolation may produce interpolated values that are outside the original range of values (by design). Periodic boundaries are currently supported only for "nearest", "linear", and 1D "cubic" interpolation.
– (string)mapColor(numeric value)
Uses the spatial map’s color-translation machinery (as defined by the valueRange and colors parameters to defineSpatialMap()) to translate each element of value into a corresponding color string. If the spatial map does not have color-translation capabilities, an error will result. See the documentation for defineSpatialMap() for information regarding the details of color translation. See the Eidos manual for further information on color strings.
– (object<Image>$)mapImage([Ni$ width = NULL], [Ni$ height = NULL], [logical$ centers = F], [logical$ color = T])
@@ -728,7 +728,7 @@The sampling of the spatial map can be done in one of two ways, as controlled by the centers parameter. If centers is T, a (width+1) × (height+1) grid of lines that delineates width × height rectangular pixels will be overlaid on top of the spatial map, and values will be sampled from the spatial map at the center of each of these pixels. If centers is F (the default), a width × height grid of lines will be overlaid on top of the spatial map, and values will be sampled from the spatial map at the vertices of the grid. If interpolation is not enabled for the spatial map, these two options will both recover the original matrix of values used to define the spatial map (assuming, here and below, that width and height are NULL). If interpolation is enabled for the spatial map, however, centers == F will recover the original values, but will not capture the “typical” value of each pixel in the image; centers == T, on the other hand, will not recover the original values, but will capture the “typical” value of each pixel in the image (i.e., the value at the center of each pixel, as produced by interpolation). The figures in section 15.11 may be helpful for visualizing the difference between these options; the overlaid grids span the full extent of the spatial map, just as shown in that section.
If color is T (the default), the valueRange and colors parameters supplied to defineSpatialMap() will be used to translate map values to RGB color values as described in the documentation of that method, providing the same appearance as in SLiMgui; of course those parameters must have been supplied, otherwise an error will result. If color is F, on the other hand, a grayscale image will be produced that directly reflects the map values without color translation. In this case, this method needs to translate map values, which can have any float value, into grayscale pixel values that are integers in [0, 255]. To do so, the map values are multiplied by 255.0, clamped to [0.0, 255.0], and then rounded to the nearest integer. This translation scheme essentially assumes that map values are in [0, 1]; for spatial maps that were defined using the floatK channel of a grayscale PNG image, this should recover the original image’s pixel values. (If a different translation scheme is desired, color=T with the desired valueRange and colors should be used.)
– (float)mapValue(float point)
-Uses the spatial map’s mapping machinery (as defined by the gridSize, values, and interpolate parameters to defineSpatialMap()) to translate the coordinates of point into a corresponding map value. The parameter map may specify the map either as a SpatialMap object, or by its string name; in either case, the map must have been added to the subpopulation. The length of point must be equal to the spatiality of the spatial map; in other words, for a spatial map with spatiality "xz", point must be of length 2, specifying the x and z coordinates of the point to be evaluated. Interpolation will automatically be used if it was enabled for the spatial map. Point coordinates are clamped into the range defined by the spatial boundaries, even if the spatial boundaries are periodic; use pointPeriodic() to wrap the point coordinates first if desired. See the documentation for defineSpatialMap() for information regarding the details of value mapping.
+Uses the spatial map’s mapping machinery (as defined by the gridSize, values, and interpolate parameters to defineSpatialMap()) to translate the coordinates of point into a corresponding map value. The length of point must be equal to the spatiality of the spatial map; in other words, for a spatial map with spatiality "xz", point must be of length 2, specifying the x and z coordinates of the point to be evaluated. Interpolation will automatically be used if it was enabled for the spatial map. Point coordinates are clamped into the range defined by the spatial boundaries, even if the spatial boundaries are periodic; use pointPeriodic() to wrap the point coordinates first if desired. See the documentation for defineSpatialMap() for information regarding the details of value mapping.
Beginning in SLiM 3.3, point may contain more than one point to be looked up. In this case, the length of point must be an exact multiple of the spatiality of the spatial map; for a spatial map with spatiality "xz", for example, the length of point must be an exact multiple of 2, and successive pairs of elements from point (elements 0 and 1, then elements 2 and 3, etc.) will be taken as the x and z coordinates of the points to be evaluated. This allows spatialMapValue() to be used in a vectorized fashion.
The spatialMapValue() method of Subpopulation provides the same functionality as this method; it may be more convenient to use, for some usage cases, and it checks that the spatial map is actually added to the subpopulation in question, providing an additional consistency check. However, either method may be used.
– (object<SpatialMap>$)multiply(ifo<SpatialMap> x)
@@ -743,13 +743,13 @@This variant of sampleNearbyPoint() samples a Metropolis–Hastings move on the spatial map. See sampleNearbyPoint() for discussion of the basic idea. This method proposes a nearby point drawn from the given kernel. If the drawn point has a larger map value than the original point, the new point is returned. If the drawn point has a smaller map value than the original point, it is returned with a probability equal to the ratio between its map value and the original map value, otherwise the original point is returned. The distribution of individuals that move (or not) to new locations governed by this method will converge upon the map itself, in a similar manner to how MCMC converges upon the posterior distribution (assuming no other forces, such as birth or death, influence the distribution of individuals). Movement governed by this method is “improved” in the sense that individuals will tend to remain where they are unless the new sampled point is an improvement for them – a higher map value.
Note that unlike sampleNearbyPoint(), this method requires that all map values are non-negative.
– (float)sampleNearbyPoint(float point, float$ maxDistance, string$ functionType, ...)
-For a spatial point supplied in point, returns a nearby point sampled from a kernel weighted by the spatial map’s values. Only points within the maximum distance of the kernel, maxDistance, will be chosen, and the probability that a given point is chosen will be proportional to the density of the kernel at that point multiplied by the value of the map at that point (interpolated, if interpolation is enabled for the map). Negative values of the map will be treated as zero.
+For a spatial point supplied in point, returns a nearby point sampled from a kernel weighted by the spatial map’s values. Only points within the maximum distance of the kernel, maxDistance, will be chosen, and the probability that a given point is chosen will be proportional to the density of the kernel at that point multiplied by the value of the map at that point (interpolated, if interpolation is enabled for the map). Negative values of the map will be treated as zero. The point returned will be within spatial bounds, respecting periodic boundaries if in effect (so there is no need to call pointPeriodic() on the result).
The kernel is specified with a kernel shape, functionType, followed by zero or more ellipsis arguments; see smooth() for further information. For this method, at present only kernel types "f", "l", "e", "n", and "t" are supported, and type "t" is not presently supported for 3D kernels.
This method can be used to find points in the vicinity of individuals that are favorable – possessing more resources, or better environmental conditions, etc. It can also be used to guide the dispersal or foraging behavior of individuals. See sampleImprovedNearbyPoint() for a variant that may be useful for directed movement across a landscape.
– (object<SpatialMap>$)smooth(float$ maxDistance, string$ functionType, ...)
Smooths (or blurs, one could say) the values of the spatial map by convolution with a kernel. The kernel is specified with a maximum distance maxDistance (beyond which the kernel cuts off to a value of zero), a kernel type functionType that should be "f", "l", "e", "n", "c", or "t", and additional parameters in the ellipsis ... that depend upon the kernel type and further specify its shape. The target spatial map is returned, to allow easy chaining of operations.
-The kernel specification is similar to that for the setInteractionType() method of InteractionType, but omits the maximum value of the kernel. Specifically, functionType may be "f", in which case no ellipsis arguments should be supplied; "l", similarly with no ellipsis arguments; "e", in which case the ellipsis should supply a numeric$ lambda (shape) parameter for a negative exponential function; "n", in which case the ellipsis should supply a numeric$ sigma (standard deviation) parameter for a Gaussian function; "c", in which case the ellipsis should supply a numeric$ scale parameter for a Cauchy distribution function; or "t", in which case the ellipsis should supply a numeric$ degrees of freedom and a numeric$ scale parameter for a t-distribution function. See the InteractionType class documentation for discussions of these kernel types.
-Distance metrics specified to this method, such as maxDistance and the additional kernel shape parameters, are measured in the distance scale of the spatial map – the same distance scale in which the spatial bounds of the map are specified. The operation is performed upon the grid values of the spatial map; distances are internally translated into the scale of the value grid. Clipping at the edge of the spatial map is done; in particular, the weights of edge and corner grid values are adjusted for their partial (one-half and one-quarter) coverage.
+The kernel specification is similar to that for the setInteractionType() method of InteractionType, but omits the maximum value of the kernel. Specifically, functionType may be "f", in which case no ellipsis arguments should be supplied; "l", similarly with no ellipsis arguments; "e", in which case the ellipsis should supply a numeric$ lambda (rate) parameter for a negative exponential function; "n", in which case the ellipsis should supply a numeric$ sigma (standard deviation) parameter for a Gaussian function; "c", in which case the ellipsis should supply a numeric$ scale parameter for a Cauchy distribution function; or "t", in which case the ellipsis should supply a numeric$ degrees of freedom and a numeric$ scale parameter for a t-distribution function. See the InteractionType class documentation for discussions of these kernel types.
+Distance metrics specified to this method, such as maxDistance and the additional kernel shape parameters, are measured in the distance scale of the spatial map – the same distance scale in which the spatial bounds of the map are specified. The operation is performed upon the grid values of the spatial map; distances are internally translated into the scale of the value grid. For non-periodic boundaries, clipping at the edge of the spatial map is done; in a 2D map with no periodic boundaries, for example, the weights of edge and corner grid values are adjusted for their partial (one-half and one-quarter) coverage. For periodic boundaries, the smoothing operation will automatically wrap around based upon the assumption that the grid values at the two connected edges of the periodic boundary have identical values (which they should, since by definition they represent the same position in space).
The density scale of the kernel has no effect and will be normalized; this is the reason that smooth(), unlike InteractionType, does not require specification of the maximum value of the kernel. This normalization prevents the kernel from increasing or decreasing the average spatial map value (apart from possible edge effects).
– (object<SpatialMap>$)subtract(ifo<SpatialMap> x)
Subtracts x from the spatial map. One possibility is that x is a singleton integer or float value; in this case, x is subtracted from each grid value of the target spatial map. Another possibility is that x is an integer or float vector/matrix/array of the same dimensions as the target spatial map’s grid; in this case, each value of x is subtracted from the corresponding grid value of the target spatial map. The third possibility is that x is itself a (singleton) spatial map; in this case, each grid value of x is subtracted from the corresponding grid value of the target spatial map (and thus the two spatial maps must match in their spatiality, their spatial bounds, and their grid dimensions). The target spatial map is returned, to allow easy chaining of operations.
diff --git a/SLiMgui/SLiMHelpClasses.rtf b/SLiMgui/SLiMHelpClasses.rtf index b1411677..20c5576f 100644 --- a/SLiMgui/SLiMHelpClasses.rtf +++ b/SLiMgui/SLiMHelpClasses.rtf @@ -5947,7 +5947,13 @@ If \f3\fs18 "linear" \f4\fs20 interpolation methods will leave the range of values in the map unchanged, \f3\fs18 "cubic" -\f4\fs20 interpolation may produce interpolated values that are outside the original range of values (by design).\ +\f4\fs20 interpolation may produce interpolated values that are outside the original range of values (by design). Periodic boundaries are currently supported only for +\f3\fs18 "nearest" +\f4\fs20 , +\f3\fs18 "linear" +\f4\fs20 , and 1D +\f3\fs18 "cubic" +\f4\fs20 interpolation.\ \pard\pardeftab397\li720\fi-446\ri720\sb180\sa60\partightenfactor0 \f3\fs18 \cf2 \'96\'a0(string)mapColor(numeric\'a0value)\ @@ -6067,13 +6073,7 @@ If \f3\fs18 defineSpatialMap() \f4\fs20 ) to translate the coordinates of \f3\fs18 point -\f4\fs20 into a corresponding map value. The parameter -\f3\fs18 map -\f4\fs20 may specify the map either as a -\f3\fs18 SpatialMap -\f4\fs20 object, or by its -\f3\fs18 string -\f4\fs20 name; in either case, the map must have been added to the subpopulation. The length of +\f4\fs20 into a corresponding map value. The length of \f3\fs18 point \f4\fs20 must be equal to the spatiality of the spatial map; in other words, for a spatial map with spatiality \f3\fs18 "xz" @@ -6230,7 +6230,9 @@ Note that unlike \f3\fs18 point \f4\fs20 , returns a nearby point sampled from a kernel weighted by the spatial map\'92s values. Only points within the maximum distance of the kernel, \f3\fs18 maxDistance -\f4\fs20 , will be chosen, and the probability that a given point is chosen will be proportional to the density of the kernel at that point multiplied by the value of the map at that point (interpolated, if interpolation is enabled for the map). Negative values of the map will be treated as zero.\ +\f4\fs20 , will be chosen, and the probability that a given point is chosen will be proportional to the density of the kernel at that point multiplied by the value of the map at that point (interpolated, if interpolation is enabled for the map). Negative values of the map will be treated as zero. The point returned will be within spatial bounds, respecting periodic boundaries if in effect (so there is no need to call +\f3\fs18 pointPeriodic() +\f4\fs20 on the result).\ The kernel is specified with a kernel shape, \f3\fs18 functionType \f4\fs20 , followed by zero or more ellipsis arguments; see @@ -6289,7 +6291,7 @@ The kernel specification is similar to that for the \f3\fs18 "e" \f4\fs20 , in which case the ellipsis should supply a \f3\fs18 numeric$ -\f4\fs20 lambda (shape) parameter for a negative exponential function; +\f4\fs20 lambda (rate) parameter for a negative exponential function; \f3\fs18 "n" \f4\fs20 , in which case the ellipsis should supply a \f3\fs18 numeric$ @@ -6310,7 +6312,7 @@ The kernel specification is similar to that for the \f4\fs20 class documentation for discussions of these kernel types.\ Distance metrics specified to this method, such as \f3\fs18 maxDistance -\f4\fs20 and the additional kernel shape parameters, are measured in the distance scale of the spatial map \'96 the same distance scale in which the spatial bounds of the map are specified. The operation is performed upon the grid values of the spatial map; distances are internally translated into the scale of the value grid. Clipping at the edge of the spatial map is done; in particular, the weights of edge and corner grid values are adjusted for their partial (one-half and one-quarter) coverage.\ +\f4\fs20 and the additional kernel shape parameters, are measured in the distance scale of the spatial map \'96 the same distance scale in which the spatial bounds of the map are specified. The operation is performed upon the grid values of the spatial map; distances are internally translated into the scale of the value grid. For non-periodic boundaries, clipping at the edge of the spatial map is done; in a 2D map with no periodic boundaries, for example, the weights of edge and corner grid values are adjusted for their partial (one-half and one-quarter) coverage. For periodic boundaries, the smoothing operation will automatically wrap around based upon the assumption that the grid values at the two connected edges of the periodic boundary have identical values (which they should, since by definition they represent the same position in space).\ The density scale of the kernel has no effect and will be normalized; this is the reason that \f3\fs18 smooth() \f4\fs20 , unlike diff --git a/core/spatial_map.cpp b/core/spatial_map.cpp index 453d855b..65415fb2 100644 --- a/core/spatial_map.cpp +++ b/core/spatial_map.cpp @@ -49,6 +49,9 @@ SpatialMap::SpatialMap(std::string p_name, std::string p_spatiality_string, Subp required_dimensionality_ = 1; spatiality_ = 1; spatiality_type_ = 1; + p_subpop->species_.SpatialPeriodicity(&periodic_a_, nullptr, nullptr); + periodic_b_ = false; + periodic_c_ = false; bounds_a0_ = p_subpop->bounds_x0_; bounds_a1_ = p_subpop->bounds_x1_; bounds_b0_ = 0; @@ -60,6 +63,9 @@ SpatialMap::SpatialMap(std::string p_name, std::string p_spatiality_string, Subp required_dimensionality_ = 2; spatiality_ = 1; spatiality_type_ = 2; + p_subpop->species_.SpatialPeriodicity(nullptr, &periodic_a_, nullptr); + periodic_b_ = false; + periodic_c_ = false; bounds_a0_ = p_subpop->bounds_y0_; bounds_a1_ = p_subpop->bounds_y1_; bounds_b0_ = 0; @@ -71,6 +77,9 @@ SpatialMap::SpatialMap(std::string p_name, std::string p_spatiality_string, Subp required_dimensionality_ = 3; spatiality_ = 1; spatiality_type_ = 3; + p_subpop->species_.SpatialPeriodicity(nullptr, nullptr, &periodic_a_); + periodic_b_ = false; + periodic_c_ = false; bounds_a0_ = p_subpop->bounds_z0_; bounds_a1_ = p_subpop->bounds_z1_; bounds_b0_ = 0; @@ -82,6 +91,8 @@ SpatialMap::SpatialMap(std::string p_name, std::string p_spatiality_string, Subp required_dimensionality_ = 2; spatiality_ = 2; spatiality_type_ = 4; + p_subpop->species_.SpatialPeriodicity(&periodic_a_, &periodic_b_, nullptr); + periodic_c_ = false; bounds_a0_ = p_subpop->bounds_x0_; bounds_a1_ = p_subpop->bounds_x1_; bounds_b0_ = p_subpop->bounds_y0_; @@ -93,6 +104,8 @@ SpatialMap::SpatialMap(std::string p_name, std::string p_spatiality_string, Subp required_dimensionality_ = 3; spatiality_ = 2; spatiality_type_ = 5; + p_subpop->species_.SpatialPeriodicity(&periodic_a_, nullptr, &periodic_b_); + periodic_c_ = false; bounds_a0_ = p_subpop->bounds_x0_; bounds_a1_ = p_subpop->bounds_x1_; bounds_b0_ = p_subpop->bounds_z0_; @@ -104,6 +117,8 @@ SpatialMap::SpatialMap(std::string p_name, std::string p_spatiality_string, Subp required_dimensionality_ = 3; spatiality_ = 2; spatiality_type_ = 6; + p_subpop->species_.SpatialPeriodicity(nullptr, &periodic_a_, &periodic_b_); + periodic_c_ = false; bounds_a0_ = p_subpop->bounds_y0_; bounds_a1_ = p_subpop->bounds_y1_; bounds_b0_ = p_subpop->bounds_z0_; @@ -115,6 +130,7 @@ SpatialMap::SpatialMap(std::string p_name, std::string p_spatiality_string, Subp required_dimensionality_ = 3; spatiality_ = 3; spatiality_type_ = 7; + p_subpop->species_.SpatialPeriodicity(&periodic_a_, &periodic_b_, &periodic_c_); bounds_a0_ = p_subpop->bounds_x0_; bounds_a1_ = p_subpop->bounds_x1_; bounds_b0_ = p_subpop->bounds_y0_; @@ -130,7 +146,7 @@ SpatialMap::SpatialMap(std::string p_name, std::string p_spatiality_string, Subp } SpatialMap::SpatialMap(std::string p_name, SpatialMap &p_original) : - name_(p_name), tag_value_(SLIM_TAG_UNSET_VALUE), spatiality_string_(p_original.spatiality_string_), spatiality_(p_original.spatiality_), spatiality_type_(p_original.spatiality_type_), required_dimensionality_(p_original.required_dimensionality_), bounds_a0_(p_original.bounds_a0_), bounds_a1_(p_original.bounds_a1_), bounds_b0_(p_original.bounds_b0_), bounds_b1_(p_original.bounds_b1_), bounds_c0_(p_original.bounds_c0_), bounds_c1_(p_original.bounds_c1_), interpolate_(p_original.interpolate_), values_min_(p_original.values_min_), values_max_(p_original.values_max_), n_colors_(p_original.n_colors_), colors_min_(p_original.colors_min_), colors_max_(p_original.colors_max_) + name_(p_name), tag_value_(SLIM_TAG_UNSET_VALUE), spatiality_string_(p_original.spatiality_string_), spatiality_(p_original.spatiality_), spatiality_type_(p_original.spatiality_type_), periodic_a_(p_original.periodic_a_), periodic_b_(p_original.periodic_b_), periodic_c_(p_original.periodic_c_), required_dimensionality_(p_original.required_dimensionality_), bounds_a0_(p_original.bounds_a0_), bounds_a1_(p_original.bounds_a1_), bounds_b0_(p_original.bounds_b0_), bounds_b1_(p_original.bounds_b1_), bounds_c0_(p_original.bounds_c0_), bounds_c1_(p_original.bounds_c1_), interpolate_(p_original.interpolate_), values_min_(p_original.values_min_), values_max_(p_original.values_max_), n_colors_(p_original.n_colors_), colors_min_(p_original.colors_min_), colors_max_(p_original.colors_max_) { // Note that this does not copy the information from EidosDictionaryRetained, and it leaves tag unset // This is intentional (that is very instance-specific state that should arguably not be copied) @@ -401,23 +417,29 @@ bool SpatialMap::IsCompatibleWithSubpopulation(Subpopulation *p_subpop) { // This checks that spatiality/dimensionality and bounds are compatible between the spatial map and a given subpopulation int spatial_dimensionality = p_subpop->species_.SpatialDimensionality(); + bool subpop_periodic_x, subpop_periodic_y, subpop_periodic_z; + + p_subpop->species_.SpatialPeriodicity(&subpop_periodic_x, &subpop_periodic_y, &subpop_periodic_z); if (spatiality_type_ == 1) { // "x" if ((required_dimensionality_ > spatial_dimensionality) || (bounds_a0_ != p_subpop->bounds_x0_) || - (bounds_a1_ != p_subpop->bounds_x1_)) + (bounds_a1_ != p_subpop->bounds_x1_) || + (periodic_a_ != subpop_periodic_x)) return false; } else if (spatiality_type_ == 2) { // "y" if ((required_dimensionality_ > spatial_dimensionality) || (bounds_a0_ != p_subpop->bounds_y0_) || - (bounds_a1_ != p_subpop->bounds_y1_)) + (bounds_a1_ != p_subpop->bounds_y1_) || + (periodic_a_ != subpop_periodic_y)) return false; } else if (spatiality_type_ == 3) { // "z" if ((required_dimensionality_ > spatial_dimensionality) || (bounds_a0_ != p_subpop->bounds_z0_) || - (bounds_a1_ != p_subpop->bounds_z1_)) + (bounds_a1_ != p_subpop->bounds_z1_) || + (periodic_a_ != subpop_periodic_z)) return false; } else if (spatiality_type_ == 4) { // "xy" @@ -425,7 +447,9 @@ bool SpatialMap::IsCompatibleWithSubpopulation(Subpopulation *p_subpop) (bounds_a0_ != p_subpop->bounds_x0_) || (bounds_a1_ != p_subpop->bounds_x1_) || (bounds_b0_ != p_subpop->bounds_y0_) || - (bounds_b1_ != p_subpop->bounds_y1_)) + (bounds_b1_ != p_subpop->bounds_y1_) || + (periodic_a_ != subpop_periodic_x) || + (periodic_b_ != subpop_periodic_y)) return false; } else if (spatiality_type_ == 5) { // "xz" @@ -433,7 +457,9 @@ bool SpatialMap::IsCompatibleWithSubpopulation(Subpopulation *p_subpop) (bounds_a0_ != p_subpop->bounds_x0_) || (bounds_a1_ != p_subpop->bounds_x1_) || (bounds_b0_ != p_subpop->bounds_z0_) || - (bounds_b1_ != p_subpop->bounds_z1_)) + (bounds_b1_ != p_subpop->bounds_z1_) || + (periodic_a_ != subpop_periodic_x) || + (periodic_b_ != subpop_periodic_z)) return false; } else if (spatiality_type_ == 6) { // "yz" @@ -441,7 +467,9 @@ bool SpatialMap::IsCompatibleWithSubpopulation(Subpopulation *p_subpop) (bounds_a0_ != p_subpop->bounds_y0_) || (bounds_a1_ != p_subpop->bounds_y1_) || (bounds_b0_ != p_subpop->bounds_z0_) || - (bounds_b1_ != p_subpop->bounds_z1_)) + (bounds_b1_ != p_subpop->bounds_z1_) || + (periodic_a_ != subpop_periodic_y) || + (periodic_b_ != subpop_periodic_z)) return false; } else if (spatiality_type_ == 7) { // "xyz" @@ -451,7 +479,10 @@ bool SpatialMap::IsCompatibleWithSubpopulation(Subpopulation *p_subpop) (bounds_b0_ != p_subpop->bounds_y0_) || (bounds_b1_ != p_subpop->bounds_y1_) || (bounds_c0_ != p_subpop->bounds_z0_) || - (bounds_c1_ != p_subpop->bounds_z1_)) + (bounds_c1_ != p_subpop->bounds_z1_) || + (periodic_a_ != subpop_periodic_x) || + (periodic_b_ != subpop_periodic_y) || + (periodic_c_ != subpop_periodic_z)) return false; } @@ -460,9 +491,11 @@ bool SpatialMap::IsCompatibleWithSubpopulation(Subpopulation *p_subpop) bool SpatialMap::IsCompatibleWithMap(SpatialMap *p_map) { - // This checks that spatiality/dimensionality and bounds are compatible between the spatial map and a given spatial map + // This checks that spatiality/dimensionality/periodicity and bounds are compatible between the spatial map and a given spatial map if ((spatiality_ != p_map->spatiality_) || (spatiality_type_ != p_map->spatiality_type_)) return false; + if ((periodic_a_ != p_map->periodic_a_) || (periodic_b_ != p_map->periodic_b_) || (periodic_c_ != p_map->periodic_c_)) + return false; if ((bounds_a0_ != p_map->bounds_a0_) || (bounds_a1_ != p_map->bounds_a1_) || (grid_size_[0] != p_map->grid_size_[0])) return false; @@ -530,6 +563,7 @@ bool SpatialMap::IsCompatibleWithValue(EidosValue *p_value) double SpatialMap::ValueAtPoint_S1(double *p_point) { // This looks up the value at point, which is in coordinates that have been normalized and clamped to [0,1] + // Note that this does NOT handle periodicity; it is assumed that the point has already been brought in bounds assert (spatiality_ == 1); double x_fraction = p_point[0]; @@ -558,6 +592,7 @@ double SpatialMap::ValueAtPoint_S1(double *p_point) double SpatialMap::ValueAtPoint_S2(double *p_point) { // This looks up the value at point, which is in coordinates that have been normalized and clamped to [0,1] + // Note that this does NOT handle periodicity; it is assumed that the point has already been brought in bounds assert (spatiality_ == 2); double x_fraction = p_point[0]; @@ -596,6 +631,7 @@ double SpatialMap::ValueAtPoint_S2(double *p_point) double SpatialMap::ValueAtPoint_S3(double *p_point) { // This looks up the value at point, which is in coordinates that have been normalized and clamped to [0,1] + // Note that this does NOT handle periodicity; it is assumed that the point has already been brought in bounds assert (spatiality_ == 3); double x_fraction = p_point[0]; @@ -751,7 +787,7 @@ void SpatialMap::Convolve_S1(SpatialKernel &kernel) // FIXME: TO BE PARALLELIZED for (int64_t a = 0; a < dim_a; ++a) { - double coverage = ((a == 0) || (a == dim_a - 1)) ? 0.5 : 1.0; + double coverage = (!periodic_a_ && ((a == 0) || (a == dim_a - 1))) ? 0.5 : 1.0; // calculate the kernel's effect at point (a) double kernel_total = 0.0; @@ -761,9 +797,18 @@ void SpatialMap::Convolve_S1(SpatialKernel &kernel) { int64_t conv_a = a + kernel_a + kernel_a_offset; - // clip to bounds + // clip/wrap to bounds if ((conv_a < 0) || (conv_a >= dim_a)) - continue; + { + if (!periodic_a_) + continue; + + // periodicity: assume the two edges have identical values, skip over the edge value on the opposite side + while (conv_a < 0) + conv_a += (dim_a - 1); // move -1 to dim - 2 + while (conv_a >= dim_a) + conv_a -= (dim_a - 1); // move dim to 1 + } // this point is within bounds; add it in to the totals double kernel_value = kernel_values[kernel_a] * coverage; @@ -810,11 +855,11 @@ void SpatialMap::Convolve_S2(SpatialKernel &kernel) // FIXME: TO BE PARALLELIZED for (int64_t b = 0; b < dim_b; ++b) { - double coverage_b = ((b == 0) || (b == dim_b - 1)) ? 0.5 : 1.0; + double coverage_b = (!periodic_b_ && ((b == 0) || (b == dim_b - 1))) ? 0.5 : 1.0; for (int64_t a = 0; a < dim_a; ++a) { - double coverage_a = ((a == 0) || (a == dim_a - 1)) ? 0.5 : 1.0; + double coverage_a = (!periodic_a_ && ((a == 0) || (a == dim_a - 1))) ? 0.5 : 1.0; double coverage = coverage_a * coverage_b; // handles partial coverage at the edges of the spatial map // calculate the kernel's effect at point (a,b) @@ -827,7 +872,16 @@ void SpatialMap::Convolve_S2(SpatialKernel &kernel) // handle bounds: either clip or wrap if ((conv_a < 0) || (conv_a >= dim_a)) - continue; + { + if (!periodic_a_) + continue; + + // periodicity: assume the two edges have identical values, skip over the edge value on the opposite side + while (conv_a < 0) + conv_a += (dim_a - 1); // move -1 to dim - 2 + while (conv_a >= dim_a) + conv_a -= (dim_a - 1); // move dim to 1 + } for (int64_t kernel_b = 0; kernel_b < kernel_dim_b; kernel_b++) { @@ -835,7 +889,16 @@ void SpatialMap::Convolve_S2(SpatialKernel &kernel) // handle bounds: either clip or wrap if ((conv_b < 0) || (conv_b >= dim_b)) - continue; + { + if (!periodic_b_) + continue; + + // periodicity: assume the two edges have identical values, skip over the edge value on the opposite side + while (conv_b < 0) + conv_b += (dim_b - 1); // move -1 to dim - 2 + while (conv_b >= dim_b) + conv_b -= (dim_b - 1); // move dim to 1 + } // this point is within bounds; add it in to the totals double kernel_value = kernel_values[kernel_a + kernel_b * kernel_dim_a] * coverage; @@ -886,15 +949,15 @@ void SpatialMap::Convolve_S3(SpatialKernel &kernel) // FIXME: TO BE PARALLELIZED for (int64_t c = 0; c < dim_c; ++c) { - double coverage_c = ((c == 0) || (c == dim_c - 1)) ? 0.5 : 1.0; + double coverage_c = (!periodic_c_ && ((c == 0) || (c == dim_c - 1))) ? 0.5 : 1.0; for (int64_t b = 0; b < dim_b; ++b) { - double coverage_b = ((b == 0) || (b == dim_b - 1)) ? 0.5 : 1.0; + double coverage_b = (!periodic_b_ && ((b == 0) || (b == dim_b - 1))) ? 0.5 : 1.0; for (int64_t a = 0; a < dim_a; ++a) { - double coverage_a = ((a == 0) || (a == dim_a - 1)) ? 0.5 : 1.0; + double coverage_a = (!periodic_a_ && ((a == 0) || (a == dim_a - 1))) ? 0.5 : 1.0; double coverage = coverage_a * coverage_b * coverage_c; // handles partial coverage at the edges of the spatial map // calculate the kernel's effect at point (a,b,c) @@ -907,7 +970,16 @@ void SpatialMap::Convolve_S3(SpatialKernel &kernel) // handle bounds: either clip or wrap if ((conv_a < 0) || (conv_a >= dim_a)) - continue; + { + if (!periodic_a_) + continue; + + // periodicity: assume the two edges have identical values, skip over the edge value on the opposite side + while (conv_a < 0) + conv_a += (dim_a - 1); // move -1 to dim - 2 + while (conv_a >= dim_a) + conv_a -= (dim_a - 1); // move dim to 1 + } for (int64_t kernel_b = 0; kernel_b < kernel_dim_b; kernel_b++) { @@ -915,7 +987,16 @@ void SpatialMap::Convolve_S3(SpatialKernel &kernel) // handle bounds: either clip or wrap if ((conv_b < 0) || (conv_b >= dim_b)) - continue; + { + if (!periodic_b_) + continue; + + // periodicity: assume the two edges have identical values, skip over the edge value on the opposite side + while (conv_b < 0) + conv_b += (dim_b - 1); // move -1 to dim - 2 + while (conv_b >= dim_b) + conv_b -= (dim_b - 1); // move dim to 1 + } for (int64_t kernel_c = 0; kernel_c < kernel_dim_c; kernel_c++) { @@ -923,7 +1004,16 @@ void SpatialMap::Convolve_S3(SpatialKernel &kernel) // handle bounds: either clip or wrap if ((conv_c < 0) || (conv_c >= dim_c)) - continue; + { + if (!periodic_c_) + continue; + + // periodicity: assume the two edges have identical values, skip over the edge value on the opposite side + while (conv_c < 0) + conv_c += (dim_c - 1); // move -1 to dim - 2 + while (conv_c >= dim_c) + conv_c -= (dim_c - 1); // move dim to 1 + } // this point is within bounds; add it in to the totals double kernel_value = kernel_values[kernel_a + kernel_b * kernel_dim_a + kernel_c * kernel_dim_a * kernel_dim_b] * coverage; @@ -1572,12 +1662,18 @@ EidosValue_SP SpatialMap::ExecuteMethod_interpolate(EidosGlobalStringID p_method else // method == 3 { // This is cubic/bicubic interpolation; we use the GSL to do this for us + // Require all/nothing for periodicity + if (((spatiality_ == 2) && (periodic_a_ != periodic_b_)) || + ((spatiality_ == 3) && ((periodic_a_ != periodic_b_) || (periodic_a_ != periodic_c_)))) + EIDOS_TERMINATION << "ERROR (SpatialMap::ExecuteMethod_interpolate): interpolate() currently requires the spatial map to be either entirely non-periodic, or entirely periodic, for 'cubic' interpolation." << EidosTerminate(nullptr); + + bool periodic = periodic_a_; // now guaranteed to apply to all dimensions + switch (spatiality_) { case 1: { // cubic interpolation - // FIXME should correctly handle periodic boundaries int64_t dim_a = (factor * (grid_size_[0] - 1)) + 1; double *new_values = (double *)malloc(dim_a * sizeof(double)); double *x = (double *)malloc(grid_size_[0] * sizeof(double)); @@ -1594,7 +1690,8 @@ EidosValue_SP SpatialMap::ExecuteMethod_interpolate(EidosGlobalStringID p_method } gsl_interp_accel *acc = gsl_interp_accel_alloc(); - gsl_spline *spline = gsl_spline_alloc(gsl_interp_cspline, grid_size_[0]); + auto interpolation_type = periodic ? gsl_interp_cspline_periodic : gsl_interp_cspline; + gsl_spline *spline = gsl_spline_alloc(interpolation_type, grid_size_[0]); double *new_values_ptr = new_values; double scale = 1.0 / factor; @@ -1615,7 +1712,14 @@ EidosValue_SP SpatialMap::ExecuteMethod_interpolate(EidosGlobalStringID p_method case 2: { // bicubic interpolation - // FIXME should correctly handle periodic boundaries + // 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); @@ -1666,7 +1770,7 @@ EidosValue_SP SpatialMap::ExecuteMethod_interpolate(EidosGlobalStringID p_method case 3: { // tricubic interpolation - not supported by the GSL - EIDOS_TERMINATION << "ERROR (SpatialMap::ExecuteMethod_interpolate): cubic interpolation is not supported for 3D spatial maps." << EidosTerminate(nullptr); + EIDOS_TERMINATION << "ERROR (SpatialMap::ExecuteMethod_interpolate): cubic interpolation is not currently supported for 3D spatial maps; please open a feature request if you need this." << EidosTerminate(nullptr); break; } default: break; @@ -1848,6 +1952,7 @@ EidosValue_SP SpatialMap::ExecuteMethod_mapValue(EidosGlobalStringID p_method_id { // We need to use the correct spatial bounds for each coordinate, which depends upon our exact spatiality // There is doubtless a way to make this code smarter, but brute force is sometimes best... + // Note that we clamp coordinates here, whether they are periodic or not; the caller should use pointPeriodic() double map_value = 0; switch (spatiality_) @@ -1970,6 +2075,13 @@ EidosValue_SP SpatialMap::ExecuteMethod_sampleImprovedNearbyPoint(EidosGlobalStr if (point_count % spatiality_ != 0) EIDOS_TERMINATION << "ERROR (SpatialMap::ExecuteMethod_sampleImprovedNearbyPoint): sampleImprovedNearbyPoint() requires the length of point to be a multiple of the spatial map's spatiality (i.e., to contain complete points)." << EidosTerminate(nullptr); + // Require all/nothing for periodicity + if (((spatiality_ == 2) && (periodic_a_ != periodic_b_)) || + ((spatiality_ == 3) && ((periodic_a_ != periodic_b_) || (periodic_a_ != periodic_c_)))) + EIDOS_TERMINATION << "ERROR (SpatialMap::ExecuteMethod_sampleImprovedNearbyPoint): sampleImprovedNearbyPoint() currently requires the spatial map to be either entirely non-periodic, or entirely periodic." << EidosTerminate(nullptr); + + bool periodic = periodic_a_; // now guaranteed to apply to all dimensions + const EidosValue_Float_vector *point_vec = (point_count == 1) ? nullptr : point_value->FloatVector(); double point_singleton = (point_count == 1) ? point_value->FloatAtIndex(0, nullptr) : 0.0; const double *point_buf = (point_count == 1) ? &point_singleton : point_vec->data(); @@ -1994,21 +2106,36 @@ EidosValue_SP SpatialMap::ExecuteMethod_sampleImprovedNearbyPoint(EidosGlobalStr // displace the point by a draw from the kernel, looping until the displaced point is in bounds double displaced_point[1]; - do + if (periodic) { + // displace the point by a draw from the kernel, then enforce periodic boundaries kernel.DrawDisplacement_S1(displaced_point); displaced_point[0] += point[0]; + + while (displaced_point[0] < 0.0) + displaced_point[0] += bounds_a1_; + while (displaced_point[0] > bounds_a1_) + displaced_point[0] -= bounds_a1_; + } + else + { + // displace the point by a draw from the kernel, looping until the displaced point is in bounds + do + { + kernel.DrawDisplacement_S1(displaced_point); + displaced_point[0] += point[0]; + } + while ((displaced_point[0] < 0.0) || (displaced_point[0] > 1.0)); } - while ((displaced_point[0] < 0.0) || (displaced_point[0] > 1.0)); // Metropolis-Hastings: move to the new point if it is better, otherwise stay with probability equal to ratio of map values double original_map_value = ValueAtPoint_S1(point); double map_value = ValueAtPoint_S1(displaced_point); if ((map_value > original_map_value) || (map_value > original_map_value * Eidos_rng_uniform(rng))) - *(result_ptr++) = displaced_point[0]; + *(result_ptr++) = displaced_point[0] * (bounds_a1_ - bounds_a0_) + bounds_a0_; else - *(result_ptr++) = point[0]; + *(result_ptr++) = point[0] * (bounds_a1_ - bounds_a0_) + bounds_a0_; } } else if (spatiality_ == 2) @@ -2026,14 +2153,35 @@ EidosValue_SP SpatialMap::ExecuteMethod_sampleImprovedNearbyPoint(EidosGlobalStr // displace the point by a draw from the kernel, looping until the displaced point is in bounds double displaced_point[2]; - do + if (periodic) { + // displace the point by a draw from the kernel, then enforce periodic boundaries kernel.DrawDisplacement_S2(displaced_point); displaced_point[0] += point[0]; displaced_point[1] += point[1]; + + while (displaced_point[0] < 0.0) + displaced_point[0] += bounds_a1_; + while (displaced_point[0] > bounds_a1_) + displaced_point[0] -= bounds_a1_; + + while (displaced_point[1] < 0.0) + displaced_point[1] += bounds_b1_; + while (displaced_point[1] > bounds_b1_) + displaced_point[1] -= bounds_b1_; + } + else + { + // displace the point by a draw from the kernel, looping until the displaced point is in bounds + do + { + kernel.DrawDisplacement_S2(displaced_point); + displaced_point[0] += point[0]; + displaced_point[1] += point[1]; + } + while ((displaced_point[0] < 0.0) || (displaced_point[0] > 1.0) || + (displaced_point[1] < 0.0) || (displaced_point[1] > 1.0)); } - while ((displaced_point[0] < 0.0) || (displaced_point[0] > 1.0) || - (displaced_point[1] < 0.0) || (displaced_point[1] > 1.0)); // Metropolis-Hastings: move to the new point if it is better, otherwise stay with probability equal to ratio of map values double original_map_value = ValueAtPoint_S2(point); @@ -2041,13 +2189,13 @@ EidosValue_SP SpatialMap::ExecuteMethod_sampleImprovedNearbyPoint(EidosGlobalStr if ((map_value > original_map_value) || (map_value > original_map_value * Eidos_rng_uniform(rng))) { - *(result_ptr++) = displaced_point[0]; - *(result_ptr++) = displaced_point[1]; + *(result_ptr++) = displaced_point[0] * (bounds_a1_ - bounds_a0_) + bounds_a0_; + *(result_ptr++) = displaced_point[1] * (bounds_b1_ - bounds_b0_) + bounds_b0_; } else { - *(result_ptr++) = point[0]; - *(result_ptr++) = point[1]; + *(result_ptr++) = point[0] * (bounds_a1_ - bounds_a0_) + bounds_a0_; + *(result_ptr++) = point[1] * (bounds_b1_ - bounds_b0_) + bounds_b0_; } } } @@ -2068,16 +2216,43 @@ EidosValue_SP SpatialMap::ExecuteMethod_sampleImprovedNearbyPoint(EidosGlobalStr // displace the point by a draw from the kernel, looping until the displaced point is in bounds double displaced_point[3]; - do + if (periodic) { + // displace the point by a draw from the kernel, then enforce periodic boundaries kernel.DrawDisplacement_S3(displaced_point); displaced_point[0] += point[0]; displaced_point[1] += point[1]; displaced_point[2] += point[2]; + + while (displaced_point[0] < 0.0) + displaced_point[0] += bounds_a1_; + while (displaced_point[0] > bounds_a1_) + displaced_point[0] -= bounds_a1_; + + while (displaced_point[1] < 0.0) + displaced_point[1] += bounds_b1_; + while (displaced_point[1] > bounds_b1_) + displaced_point[1] -= bounds_b1_; + + while (displaced_point[2] < 0.0) + displaced_point[2] += bounds_c1_; + while (displaced_point[2] > bounds_c1_) + displaced_point[2] -= bounds_c1_; + } + else + { + // displace the point by a draw from the kernel, looping until the displaced point is in bounds + do + { + kernel.DrawDisplacement_S3(displaced_point); + displaced_point[0] += point[0]; + displaced_point[1] += point[1]; + displaced_point[2] += point[2]; + } + while ((displaced_point[0] < 0.0) || (displaced_point[0] > 1.0) || + (displaced_point[1] < 0.0) || (displaced_point[1] > 1.0) || + (displaced_point[2] < 0.0) || (displaced_point[2] > 1.0)); } - while ((displaced_point[0] < 0.0) || (displaced_point[0] > 1.0) || - (displaced_point[1] < 0.0) || (displaced_point[1] > 1.0) || - (displaced_point[2] < 0.0) || (displaced_point[2] > 1.0)); // Metropolis-Hastings: move to the new point if it is better, otherwise stay with probability equal to ratio of map values double original_map_value = ValueAtPoint_S3(point); @@ -2085,15 +2260,15 @@ EidosValue_SP SpatialMap::ExecuteMethod_sampleImprovedNearbyPoint(EidosGlobalStr if ((map_value > original_map_value) || (map_value > original_map_value * Eidos_rng_uniform(rng))) { - *(result_ptr++) = displaced_point[0]; - *(result_ptr++) = displaced_point[1]; - *(result_ptr++) = displaced_point[2]; + *(result_ptr++) = displaced_point[0] * (bounds_a1_ - bounds_a0_) + bounds_a0_; + *(result_ptr++) = displaced_point[1] * (bounds_b1_ - bounds_b0_) + bounds_b0_; + *(result_ptr++) = displaced_point[2] * (bounds_c1_ - bounds_c0_) + bounds_c0_; } else { - *(result_ptr++) = point[0]; - *(result_ptr++) = point[1]; - *(result_ptr++) = point[2]; + *(result_ptr++) = point[0] * (bounds_a1_ - bounds_a0_) + bounds_a0_; + *(result_ptr++) = point[1] * (bounds_b1_ - bounds_b0_) + bounds_b0_; + *(result_ptr++) = point[2] * (bounds_c1_ - bounds_c0_) + bounds_c0_; } } } @@ -2118,6 +2293,13 @@ EidosValue_SP SpatialMap::ExecuteMethod_sampleNearbyPoint(EidosGlobalStringID p_ if (point_count % spatiality_ != 0) EIDOS_TERMINATION << "ERROR (SpatialMap::ExecuteMethod_sampleNearbyPoint): sampleNearbyPoint() requires the length of point to be a multiple of the spatial map's spatiality (i.e., to contain complete points)." << EidosTerminate(nullptr); + // Require all/nothing for periodicity + if (((spatiality_ == 2) && (periodic_a_ != periodic_b_)) || + ((spatiality_ == 3) && ((periodic_a_ != periodic_b_) || (periodic_a_ != periodic_c_)))) + EIDOS_TERMINATION << "ERROR (SpatialMap::ExecuteMethod_sampleNearbyPoint): sampleNearbyPoint() currently requires the spatial map to be either entirely non-periodic, or entirely periodic." << EidosTerminate(nullptr); + + bool periodic = periodic_a_; // now guaranteed to apply to all dimensions + const EidosValue_Float_vector *point_vec = (point_count == 1) ? nullptr : point_value->FloatVector(); double point_singleton = (point_count == 1) ? point_value->FloatAtIndex(0, nullptr) : 0.0; const double *point_buf = (point_count == 1) ? &point_singleton : point_vec->data(); @@ -2144,19 +2326,34 @@ EidosValue_SP SpatialMap::ExecuteMethod_sampleNearbyPoint(EidosGlobalStringID p_ // rejection sample to draw a displaced point from the product of the kernel times the map do { - // displace the point by a draw from the kernel, looping until the displaced point is in bounds - do + if (periodic) { + // displace the point by a draw from the kernel, then enforce periodic boundaries kernel.DrawDisplacement_S1(displaced_point); displaced_point[0] += point_a; + + while (displaced_point[0] < 0.0) + displaced_point[0] += bounds_a1_; + while (displaced_point[0] > bounds_a1_) + displaced_point[0] -= bounds_a1_; + } + else + { + // displace the point by a draw from the kernel, looping until the displaced point is in bounds + do + { + kernel.DrawDisplacement_S1(displaced_point); + displaced_point[0] += point_a; + } + while ((displaced_point[0] < 0.0) || (displaced_point[0] > 1.0)); } - while ((displaced_point[0] < 0.0) || (displaced_point[0] > 1.0)); map_value = ValueAtPoint_S1(displaced_point); } while (values_max_ * Eidos_rng_uniform(rng) > map_value); - *(result_ptr++) = displaced_point[0]; + // scale point coordinates back to user space + *(result_ptr++) = displaced_point[0] * (bounds_a1_ - bounds_a0_) + bounds_a0_; } } else if (spatiality_ == 2) @@ -2176,22 +2373,43 @@ EidosValue_SP SpatialMap::ExecuteMethod_sampleNearbyPoint(EidosGlobalStringID p_ // rejection sample to draw a displaced point from the product of the kernel times the map do { - // displace the point by a draw from the kernel, looping until the displaced point is in bounds - do + if (periodic) { + // displace the point by a draw from the kernel, then enforce periodic boundaries kernel.DrawDisplacement_S2(displaced_point); displaced_point[0] += point_a; displaced_point[1] += point_b; + + while (displaced_point[0] < 0.0) + displaced_point[0] += bounds_a1_; + while (displaced_point[0] > bounds_a1_) + displaced_point[0] -= bounds_a1_; + + while (displaced_point[1] < 0.0) + displaced_point[1] += bounds_b1_; + while (displaced_point[1] > bounds_b1_) + displaced_point[1] -= bounds_b1_; + } + else + { + // displace the point by a draw from the kernel, looping until the displaced point is in bounds + do + { + kernel.DrawDisplacement_S2(displaced_point); + displaced_point[0] += point_a; + displaced_point[1] += point_b; + } + while ((displaced_point[0] < 0.0) || (displaced_point[0] > 1.0) || + (displaced_point[1] < 0.0) || (displaced_point[1] > 1.0)); } - while ((displaced_point[0] < 0.0) || (displaced_point[0] > 1.0) || - (displaced_point[1] < 0.0) || (displaced_point[1] > 1.0)); map_value = ValueAtPoint_S2(displaced_point); } while (values_max_ * Eidos_rng_uniform(rng) > map_value); - *(result_ptr++) = displaced_point[0]; - *(result_ptr++) = displaced_point[1]; + // scale point coordinates back to user space + *(result_ptr++) = displaced_point[0] * (bounds_a1_ - bounds_a0_) + bounds_a0_; + *(result_ptr++) = displaced_point[1] * (bounds_b1_ - bounds_b0_) + bounds_b0_; } } else // (spatiality_ == 3) @@ -2213,25 +2431,52 @@ EidosValue_SP SpatialMap::ExecuteMethod_sampleNearbyPoint(EidosGlobalStringID p_ // rejection sample to draw a displaced point from the product of the kernel times the map do { - // displace the point by a draw from the kernel, looping until the displaced point is in bounds - do + if (periodic) { + // displace the point by a draw from the kernel, then enforce periodic boundaries kernel.DrawDisplacement_S3(displaced_point); displaced_point[0] += point_a; displaced_point[1] += point_b; displaced_point[2] += point_c; + + while (displaced_point[0] < 0.0) + displaced_point[0] += bounds_a1_; + while (displaced_point[0] > bounds_a1_) + displaced_point[0] -= bounds_a1_; + + while (displaced_point[1] < 0.0) + displaced_point[1] += bounds_b1_; + while (displaced_point[1] > bounds_b1_) + displaced_point[1] -= bounds_b1_; + + while (displaced_point[2] < 0.0) + displaced_point[2] += bounds_c1_; + while (displaced_point[2] > bounds_c1_) + displaced_point[2] -= bounds_c1_; + } + else + { + // displace the point by a draw from the kernel, looping until the displaced point is in bounds + do + { + kernel.DrawDisplacement_S3(displaced_point); + displaced_point[0] += point_a; + displaced_point[1] += point_b; + displaced_point[2] += point_c; + } + while ((displaced_point[0] < 0.0) || (displaced_point[0] > 1.0) || + (displaced_point[1] < 0.0) || (displaced_point[1] > 1.0) || + (displaced_point[2] < 0.0) || (displaced_point[2] > 1.0)); } - while ((displaced_point[0] < 0.0) || (displaced_point[0] > 1.0) || - (displaced_point[1] < 0.0) || (displaced_point[1] > 1.0) || - (displaced_point[2] < 0.0) || (displaced_point[2] > 1.0)); map_value = ValueAtPoint_S3(displaced_point); } while (values_max_ * Eidos_rng_uniform(rng) > map_value); - *(result_ptr++) = displaced_point[0]; - *(result_ptr++) = displaced_point[1]; - *(result_ptr++) = displaced_point[2]; + // scale point coordinates back to user space + *(result_ptr++) = displaced_point[0] * (bounds_a1_ - bounds_a0_) + bounds_a0_; + *(result_ptr++) = displaced_point[1] * (bounds_b1_ - bounds_b0_) + bounds_b0_; + *(result_ptr++) = displaced_point[2] * (bounds_c1_ - bounds_c0_) + bounds_c0_; } } diff --git a/core/spatial_map.h b/core/spatial_map.h index ba08a535..83703f83 100644 --- a/core/spatial_map.h +++ b/core/spatial_map.h @@ -63,6 +63,7 @@ class SpatialMap : public EidosDictionaryRetained std::string spatiality_string_; // "x", "y", "z", "xy", "xz", "yz", or "xyz": the spatial dimensions for the map int spatiality_; // 1, 2, or 3 for 1D, 2D, or 3D: the number of spatial dimensions int spatiality_type_; // 1=="x", 2=="y", 3=="z", 4=="xy", 5=="yz", 6=="xz", 7=="xyz" + bool periodic_a_, periodic_b_, periodic_c_; // periodic boundary flags for spatiality dimensions a/b/c int required_dimensionality_; // 1, 2, or 3 for the dimensionality we require; enough to encompass spatiality_type_