From ca76d7692af1e184a27930d83fd5cae044292839 Mon Sep 17 00:00:00 2001 From: Ben Haller Date: Sun, 10 Sep 2023 13:08:32 -0400 Subject: [PATCH] add smoothValues() method and SpatialKernel class --- QtSLiM/QtSLiMTablesDrawer.cpp | 20 +- QtSLiM/help/SLiMHelpClasses.html | 12 +- SLiM.xcodeproj/project.pbxproj | 12 + SLiMgui/PopulationView.mm | 112 ++++++++++ SLiMgui/SLiMHelpClasses.rtf | 60 ++++- SLiMgui/SLiMWindowController.mm | 20 +- VERSIONS | 1 + core/interaction_type.cpp | 168 ++++---------- core/interaction_type.h | 16 +- core/slim_globals.cpp | 1 + core/slim_globals.h | 2 + core/spatial_kernel.cpp | 368 +++++++++++++++++++++++++++++++ core/spatial_kernel.h | 90 ++++++++ core/spatial_map.cpp | 296 +++++++++++++++++++++++-- core/spatial_map.h | 8 +- 15 files changed, 987 insertions(+), 199 deletions(-) create mode 100644 core/spatial_kernel.cpp create mode 100644 core/spatial_kernel.h diff --git a/QtSLiM/QtSLiMTablesDrawer.cpp b/QtSLiM/QtSLiMTablesDrawer.cpp index a133f5100..654403be1 100644 --- a/QtSLiM/QtSLiMTablesDrawer.cpp +++ b/QtSLiM/QtSLiMTablesDrawer.cpp @@ -906,11 +906,11 @@ QVariant QtSLiMInteractionTypeTableModel::data(const QModelIndex &p_index, int r { switch (interactionType->if_type_) { - case IFType::kFixed: return QVariant(QString("fixed")); - case IFType::kLinear: return QVariant(QString("linear")); - case IFType::kExponential: return QVariant(QString("exp")); - case IFType::kNormal: return QVariant(QString("normal")); - case IFType::kCauchy: return QVariant(QString("Cauchy")); + case SpatialKernelType::kFixed: return QVariant(QString("fixed")); + case SpatialKernelType::kLinear: return QVariant(QString("linear")); + case SpatialKernelType::kExponential: return QVariant(QString("exp")); + case SpatialKernelType::kNormal: return QVariant(QString("normal")); + case SpatialKernelType::kCauchy: return QVariant(QString("Cauchy")); } } else if (p_index.column() == 3) @@ -923,16 +923,16 @@ QVariant QtSLiMInteractionTypeTableModel::data(const QModelIndex &p_index, int r // append second parameters where applicable switch (interactionType->if_type_) { - case IFType::kFixed: - case IFType::kLinear: + case SpatialKernelType::kFixed: + case SpatialKernelType::kLinear: break; - case IFType::kExponential: + case SpatialKernelType::kExponential: paramString += QString(", β=%1").arg(interactionType->if_param2_, 0, 'f', 3); break; - case IFType::kNormal: + case SpatialKernelType::kNormal: paramString += QString(", σ=%1").arg(interactionType->if_param2_, 0, 'f', 3); break; - case IFType::kCauchy: + case SpatialKernelType::kCauchy: paramString += QString(", γ=%1").arg(interactionType->if_param2_, 0, 'f', 3); break; } diff --git a/QtSLiM/help/SLiMHelpClasses.html b/QtSLiM/help/SLiMHelpClasses.html index 34e20f887..0008bd597 100644 --- a/QtSLiM/help/SLiMHelpClasses.html +++ b/QtSLiM/help/SLiMHelpClasses.html @@ -701,10 +701,10 @@

5.14.2  SpatialMap methods

– (void)changeValues(numeric values)

Changes the grid values used for the target spatial map.  The values parameter should be a vector, matrix, or array of numeric values as described in the documentation for defineSpatialMap().  Other characteristics of the spatial map, such as its color mapping (if defined), its spatial bounds, and its spatiality, will remain unchanged.  The grid resolution of the spatial map is allowed to change with this method.  This method is useful for changing the values of a spatial map over time, such as to implement changes to the landscape’s characteristics due to seasonality, climate change, processes such as fire or urbanization, and so forth.  As with the original map values provided to defineSpatialMap(), it is often useful to read map values from a PNG image file using the Eidos class Image.

-

– (void)interpolateValues(integer$ factor)

-

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, 100].

-

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 kN−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.

-

Linear (or bilinear or trilinear) interpolation will be used to interpolate the values for the new grid points; this is done even if interpolation is not normally enabled for the spatial map (i.e., if the map’s interpolate property is F).

+

– (void)interpolateValues(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, 201], somewhat arbitrarily

+

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) 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.  Other interpolation methods may be added in future.  The choice of interpolation method used here is independent of the map’s interpolate property.

– (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])

@@ -715,6 +715,10 @@

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.

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.

+

– (void)smoothValues(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", or "c", and additional parameters in the ellipsis that depend upon the kernel type and further specify its shape.  The kernel types and their additional parameters are exactly the same as for the setInteractionType() method of InteractionType; see the documentation of that method for details.

+

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.

+

The density scale of the kernel is unimportant and will be normalized; for example, a kernel specified to have a maximum value of 1, versus a kernel specified to have a maximum value of 10, will produce the same effect.  This normalization prevents the kernel from increasing or decreasing the average spatial map value (apart from possible edge effects).

5.15  Class Species

5.15.1  Species properties

avatar => (string$)

diff --git a/SLiM.xcodeproj/project.pbxproj b/SLiM.xcodeproj/project.pbxproj index 41b7891cb..36fc59a48 100644 --- a/SLiM.xcodeproj/project.pbxproj +++ b/SLiM.xcodeproj/project.pbxproj @@ -338,6 +338,10 @@ 985F3F0E24BA31D900E712E0 /* eidos_test_functions_statistics.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 985F3F0A24BA31D900E712E0 /* eidos_test_functions_statistics.cpp */; }; 98606AEE1DED0DCD00821CFF /* mutation_run.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 98606AEC1DED0DCD00821CFF /* mutation_run.cpp */; }; 98606AEF1DED0DCD00821CFF /* mutation_run.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 98606AEC1DED0DCD00821CFF /* mutation_run.cpp */; }; + 986070EA2AACECD600FD6143 /* spatial_kernel.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 986070E82AACECD600FD6143 /* spatial_kernel.cpp */; }; + 986070EB2AACECD600FD6143 /* spatial_kernel.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 986070E82AACECD600FD6143 /* spatial_kernel.cpp */; }; + 986070EC2AACECD600FD6143 /* spatial_kernel.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 986070E82AACECD600FD6143 /* spatial_kernel.cpp */; }; + 986070ED2AACECD600FD6143 /* spatial_kernel.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 986070E82AACECD600FD6143 /* spatial_kernel.cpp */; }; 986926D41AA1337A0000E138 /* graph_submenu_H.pdf in Resources */ = {isa = PBXBuildFile; fileRef = 986926D21AA1337A0000E138 /* graph_submenu_H.pdf */; }; 986926D51AA1337A0000E138 /* graph_submenu.pdf in Resources */ = {isa = PBXBuildFile; fileRef = 986926D31AA1337A0000E138 /* graph_submenu.pdf */; }; 986926D91AA140550000E138 /* GraphView.mm in Sources */ = {isa = PBXBuildFile; fileRef = 986926D81AA140550000E138 /* GraphView.mm */; }; @@ -1612,6 +1616,8 @@ 986021081A3BB504001BDCFE /* LICENSE */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = text; path = LICENSE; sourceTree = ""; }; 98606AEC1DED0DCD00821CFF /* mutation_run.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = mutation_run.cpp; sourceTree = ""; }; 98606AED1DED0DCD00821CFF /* mutation_run.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = mutation_run.h; sourceTree = ""; }; + 986070E82AACECD600FD6143 /* spatial_kernel.cpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.cpp; path = spatial_kernel.cpp; sourceTree = ""; }; + 986070E92AACECD600FD6143 /* spatial_kernel.h */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.c.h; path = spatial_kernel.h; sourceTree = ""; }; 986764A92089066A00E81B2E /* tables.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; name = tables.h; path = treerec/tskit/tables.h; sourceTree = ""; }; 986926D21AA1337A0000E138 /* graph_submenu_H.pdf */ = {isa = PBXFileReference; lastKnownFileType = image.pdf; path = graph_submenu_H.pdf; sourceTree = ""; }; 986926D31AA1337A0000E138 /* graph_submenu.pdf */ = {isa = PBXFileReference; lastKnownFileType = image.pdf; path = graph_submenu.pdf; sourceTree = ""; }; @@ -2186,6 +2192,8 @@ 98E9A69A1A3CD542000AD4FC /* polymorphism.cpp */, 985D1D8A2808B84F00461CFA /* sparse_vector.h */, 985D1D892808B84F00461CFA /* sparse_vector.cpp */, + 986070E92AACECD600FD6143 /* spatial_kernel.h */, + 986070E82AACECD600FD6143 /* spatial_kernel.cpp */, 98DB3D6E1E6122AE00E2C200 /* interaction_type.h */, 98DB3D6D1E6122AE00E2C200 /* interaction_type.cpp */, 98DEB47D2AA632AA00ABE60F /* spatial_map.h */, @@ -3466,6 +3474,7 @@ 981DC36C28E26F8B000ABE91 /* eidos_functions_other.cpp in Sources */, 98332AC21FDBA53F00274FF0 /* xerbla.c in Sources */, 985724D51AD481070047C223 /* eidos_test.cpp in Sources */, + 986070EA2AACECD600FD6143 /* spatial_kernel.cpp in Sources */, 98332AB41FDBA1E100274FF0 /* blas.c in Sources */, 9809DFA02550F32500C4E82D /* log_file.cpp in Sources */, 989A5BE92525304100E7192D /* eidos_class_Dictionary.cpp in Sources */, @@ -3779,6 +3788,7 @@ 98CF5251294A3FC900557BBA /* eidos_class_Dictionary.cpp in Sources */, 98CF5252294A3FC900557BBA /* elementary.c in Sources */, 98CF5253294A3FC900557BBA /* SLiMHaplotypeManager.mm in Sources */, + 986070EC2AACECD600FD6143 /* spatial_kernel.cpp in Sources */, 98CF5254294A3FC900557BBA /* GraphView_FitnessOverTime.mm in Sources */, 98CF5255294A3FC900557BBA /* SLiMWindowController.mm in Sources */, 98CF5256294A3FC900557BBA /* PopulationView.mm in Sources */, @@ -4087,6 +4097,7 @@ 989A5BEA2525304100E7192D /* eidos_class_Dictionary.cpp in Sources */, 98C821561C7A983700548839 /* elementary.c in Sources */, 98B40DFA1FB2DA8C0046F695 /* SLiMHaplotypeManager.mm in Sources */, + 986070EB2AACECD600FD6143 /* spatial_kernel.cpp in Sources */, 986926E81AA40AFF0000E138 /* GraphView_FitnessOverTime.mm in Sources */, 98D4C2041A701D5A00FFB083 /* SLiMWindowController.mm in Sources */, 98D4C2071A704EA700FFB083 /* PopulationView.mm in Sources */, @@ -4367,6 +4378,7 @@ 98D7ECE228CE58FC00DEAAC4 /* eidos_token.cpp in Sources */, 98D7ECE328CE58FC00DEAAC4 /* core.c in Sources */, 98D7ECE428CE58FC00DEAAC4 /* minmax.c in Sources */, + 986070ED2AACECD600FD6143 /* spatial_kernel.cpp in Sources */, 98D7ECE528CE58FC00DEAAC4 /* trees.c in Sources */, 98D7ECE628CE58FC00DEAAC4 /* xerbla.c in Sources */, 98D7ECE728CE58FC00DEAAC4 /* eidos_test.cpp in Sources */, diff --git a/SLiMgui/PopulationView.mm b/SLiMgui/PopulationView.mm index a7feb4c74..e4c3cc8fb 100644 --- a/SLiMgui/PopulationView.mm +++ b/SLiMgui/PopulationView.mm @@ -892,6 +892,118 @@ - (void)_drawBackgroundSpatialMap:(SpatialMap *)background_map inBounds:(NSRect) glDisableClientState(GL_VERTEX_ARRAY); glDisableClientState(GL_COLOR_ARRAY); + +#if 0 + // Experimental feature: draw boxes showing where the grid nodes are, since that is rather confusing! + NSRect individualArea = NSMakeRect(bounds.origin.x, bounds.origin.y, bounds.size.width - 1, bounds.size.height - 1); + int64_t xsize = background_map->grid_size_[0]; + int64_t ysize = background_map->grid_size_[1]; + double *values = background_map->values_; + + if ((xsize <= 51) && (ysize <= 51)) + { + // Set up to draw rects + displayListIndex = 0; + + vertices = glArrayVertices; + glEnableClientState(GL_VERTEX_ARRAY); + glVertexPointer(2, GL_FLOAT, 0, glArrayVertices); + + colors = glArrayColors; + glEnableClientState(GL_COLOR_ARRAY); + glColorPointer(4, GL_FLOAT, 0, glArrayColors); + + // first pass we draw squares to make outlines, second pass we draw the interiors in color + for (int pass = 0; pass <= 1; ++pass) + { + for (int x = 0; x < xsize; ++x) + { + for (int y = 0; y < ysize; ++y) + { + float position_x = x / (float)(xsize - 1); // 0 to 1 + float position_y = y / (float)(ysize - 1); // 0 to 1 + + float centerX = (float)(individualArea.origin.x + round(position_x * individualArea.size.width) + 0.5); + float centerY = (float)(individualArea.origin.y + individualArea.size.height - round(position_y * individualArea.size.height) + 0.5); + const float margin = ((pass == 0) ? 5.5f : 3.5f); + float left = centerX - margin; + float top = centerY - margin; + float right = centerX + margin; + float bottom = centerY + margin; + + if (left < individualArea.origin.x) + left = (float)individualArea.origin.x; + if (top < individualArea.origin.y) + top = (float)individualArea.origin.y; + if (right > individualArea.origin.x + individualArea.size.width) + right = (float)(individualArea.origin.x + individualArea.size.width); + if (bottom > individualArea.origin.y + individualArea.size.height) + bottom = (float)(individualArea.origin.y + individualArea.size.height); + + *(vertices++) = left; + *(vertices++) = top; + *(vertices++) = left; + *(vertices++) = bottom; + *(vertices++) = right; + *(vertices++) = bottom; + *(vertices++) = right; + *(vertices++) = top; + + if (pass == 0) + { + for (int j = 0; j < 4; ++j) + { + *(colors++) = 1.0; + *(colors++) = 0.25; + *(colors++) = 0.25; + *(colors++) = 1.0; + } + } + else + { + // look up the map's color at this grid point + float rgb[3]; + double value = values[x + y * xsize]; + + background_map->ColorForValue(value, rgb); + + for (int j = 0; j < 4; ++j) + { + *(colors++) = rgb[0]; + *(colors++) = rgb[1]; + *(colors++) = rgb[2]; + *(colors++) = 1.0; + } + } + + displayListIndex++; + + // If we've filled our buffers, get ready to draw more + if (displayListIndex == kMaxGLRects) + { + // Draw our arrays + glDrawArrays(GL_QUADS, 0, 4 * displayListIndex); + + // And get ready to draw more + vertices = glArrayVertices; + colors = glArrayColors; + displayListIndex = 0; + } + } + } + } + + // Draw any leftovers + if (displayListIndex) + { + // Draw our arrays + glDrawArrays(GL_QUADS, 0, 4 * displayListIndex); + } + + glDisableClientState(GL_VERTEX_ARRAY); + glDisableClientState(GL_COLOR_ARRAY); + } +#endif } - (void)chooseDefaultBackgroundSettings:(PopulationViewBackgroundSettings *)background map:(SpatialMap **)returnMap forSubpopulation:(Subpopulation *)subpop diff --git a/SLiMgui/SLiMHelpClasses.rtf b/SLiMgui/SLiMHelpClasses.rtf index 696e8353a..bf894d4da 100644 --- a/SLiMgui/SLiMHelpClasses.rtf +++ b/SLiMgui/SLiMHelpClasses.rtf @@ -5718,7 +5718,7 @@ If the model is being profiled, or is executing forward to a tick number entered \f4\fs20 .\ \pard\pardeftab397\li720\fi-446\ri720\sb180\sa60\partightenfactor0 -\f3\fs18 \cf2 \'96\'a0(void)interpolateValues(integer$\'a0factor)\ +\f3\fs18 \cf2 \'96\'a0(void)interpolateValues(integer$\'a0factor, [string$\'a0method\'a0=\'a0"linear"])\ \pard\pardeftab720\li547\ri720\sb60\sa60\partightenfactor0 \f4\fs20 \cf2 Increases the resolution of the spatial map by @@ -5728,8 +5728,8 @@ If the model is being profiled, or is executing forward to a tick number entered \f4\fs20 must be an integer in [ \f3\fs18 2 \f4\fs20 , -\f3\fs18 100 -\f4\fs20 ].\ +\f3\fs18 201 +\f4\fs20 ], somewhat arbitrarily\ For a 1D spatial map, \f3\fs18 factor-1 \f4\fs20 new values will be inserted between every pair of values in the original value grid. A @@ -5743,19 +5743,29 @@ For a 1D spatial map, \f4\fs20 of \f1\i k \f4\i0 would be -\f1\i kN -\f4\i0 \uc0\u8722 1, not +\f1\i k +\f4\i0 ( +\f1\i N +\f4\i0 \uc0\u8722 1)+1, not \f1\i kN \f4\i0 , because new values are inserted only \f1\i between \f4\i0 existing values. For 2D and 3D spatial maps, essentially the same process is conducted along each axis of the map\'92s spatiality, increasing the resolution of the map by \f3\fs18 factor \f4\fs20 in every dimension.\ -Linear (or bilinear or trilinear) interpolation will be used to interpolate the values for the new grid points; this is done even if interpolation is not normally enabled for the spatial map (i.e., if the map\'92s +If +\f3\fs18 method +\f4\fs20 is +\f3\fs18 "linear" +\f4\fs20 (the default), linear (or bilinear or trilinear) interpolation will be used to interpolate the values for the new grid points. Alternatively, if +\f3\fs18 method +\f4\fs20 is +\f3\fs18 "nearest" +\f4\fs20 , the nearest value in the old grid will be used for new grid points; with this method, it is recommended that +\f3\fs18 factor +\f4\fs20 be odd, not even, to avoid artifacts due to rounding of coordinates midway between the original grid positions. Other interpolation methods may be added in future. The choice of interpolation method used here is independent of the map\'92s \f3\fs18 interpolate -\f4\fs20 property is -\f3\fs18 F -\f4\fs20 ).\ +\f4\fs20 property.\ \pard\pardeftab397\li720\fi-446\ri720\sb180\sa60\partightenfactor0 \f3\fs18 \cf2 \'96\'a0(string)mapColor(numeric\'a0value)\ @@ -5928,6 +5938,38 @@ The \f4\fs20 method of \f3\fs18 Subpopulation \f4\fs20 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.\ +\pard\pardeftab397\li720\fi-446\ri720\sb180\sa60\partightenfactor0 + +\f3\fs18 \cf2 \'96\'a0(void)smoothValues(float$\'a0maxDistance, string$\'a0functionType, ...)\ +\pard\pardeftab720\li547\ri720\sb60\sa60\partightenfactor0 + +\f4\fs20 \cf2 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 +\f3\fs18 maxDistance +\f4\fs20 (beyond which the kernel cuts off to a value of zero), a kernel type +\f3\fs18 functionType +\f4\fs20 that should be +\f3\fs18 "f" +\f4\fs20 , +\f3\fs18 "l" +\f4\fs20 , +\f3\fs18 "e" +\f4\fs20 , +\f3\fs18 "n" +\f4\fs20 , or +\f3\fs18 "c" +\f4\fs20 , and additional parameters in the ellipsis that depend upon the kernel type and further specify its shape. The kernel types and their additional parameters are exactly the same as for the +\f3\fs18 setInteractionType() +\f4\fs20 method of +\f3\fs18 InteractionType +\f4\fs20 ; see the documentation of that method for details.\ +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.\ +The density scale of the kernel is unimportant and will be normalized; for example, a kernel specified to have a maximum value of +\f3\fs18 1 +\f4\fs20 , versus a kernel specified to have a maximum value of +\f3\fs18 10 +\f4\fs20 , will produce the same effect. This normalization prevents the kernel from increasing or decreasing the average spatial map value (apart from possible edge effects).\ \pard\pardeftab720\ri720\sb360\sa60\partightenfactor0 \f0\b\fs22 \cf0 5.15 Class Species\ diff --git a/SLiMgui/SLiMWindowController.mm b/SLiMgui/SLiMWindowController.mm index a93599b31..30dd47550 100644 --- a/SLiMgui/SLiMWindowController.mm +++ b/SLiMgui/SLiMWindowController.mm @@ -4941,11 +4941,11 @@ - (id)tableView:(NSTableView *)aTableView objectValueForTableColumn:(NSTableColu { switch (interactionType->if_type_) { - case IFType::kFixed: return @"fixed"; - case IFType::kLinear: return @"linear"; - case IFType::kExponential: return @"exp"; - case IFType::kNormal: return @"normal"; - case IFType::kCauchy: return @"Cauchy"; + case SpatialKernelType::kFixed: return @"fixed"; + case SpatialKernelType::kLinear: return @"linear"; + case SpatialKernelType::kExponential: return @"exp"; + case SpatialKernelType::kNormal: return @"normal"; + case SpatialKernelType::kCauchy: return @"Cauchy"; } } else if (aTableColumn == interactionTypeIFParamsColumn) @@ -4958,16 +4958,16 @@ - (id)tableView:(NSTableView *)aTableView objectValueForTableColumn:(NSTableColu // append second parameters where applicable switch (interactionType->if_type_) { - case IFType::kFixed: - case IFType::kLinear: + case SpatialKernelType::kFixed: + case SpatialKernelType::kLinear: break; - case IFType::kExponential: + case SpatialKernelType::kExponential: [paramString appendFormat:@", β=%.3f", interactionType->if_param2_]; break; - case IFType::kNormal: + case SpatialKernelType::kNormal: [paramString appendFormat:@", σ=%.3f", interactionType->if_param2_]; break; - case IFType::kCauchy: + case SpatialKernelType::kCauchy: [paramString appendFormat:@", γ=%.3f", interactionType->if_param2_]; break; } diff --git a/VERSIONS b/VERSIONS index 12277b6e0..1e9e170ad 100644 --- a/VERSIONS +++ b/VERSIONS @@ -147,6 +147,7 @@ development head (in the master branch): add a changeValues() method on SpatialMap to change the grid values of an existing spatial map add various properties on SpatialMap: gridDimensions, interpolate, name, spatialBounds, spatiality, tag add interpolateValues() to SpatialMap, to increase the resolution of a spatial map by interpolating new values + add smoothValues() to SpatialMap, to smooth/blur a map; internally, add the SpatialKernel class (not user-visible) version 4.0.1 (Eidos version 3.0.1): diff --git a/core/interaction_type.cpp b/core/interaction_type.cpp index 0a290cfee..551a71ab9 100755 --- a/core/interaction_type.cpp +++ b/core/interaction_type.cpp @@ -33,22 +33,6 @@ #include -// stream output for enumerations -std::ostream& operator<<(std::ostream& p_out, IFType p_if_type) -{ - switch (p_if_type) - { - case IFType::kFixed: p_out << gStr_f; break; - case IFType::kLinear: p_out << gStr_l; break; - case IFType::kExponential: p_out << gStr_e; break; - case IFType::kNormal: p_out << gEidosStr_n; break; - case IFType::kCauchy: p_out << gEidosStr_c; break; - } - - return p_out; -} - - #pragma mark - #pragma mark InteractionType #pragma mark - @@ -93,7 +77,7 @@ void InteractionType::_WarmUp(void) InteractionType::InteractionType(Community &p_community, slim_objectid_t p_interaction_type_id, std::string p_spatiality_string, bool p_reciprocal, double p_max_distance, IndividualSex p_receiver_sex, IndividualSex p_exerter_sex) : self_symbol_(EidosStringRegistry::GlobalStringIDForString(SLiMEidosScript::IDStringWithPrefix('i', p_interaction_type_id)), EidosValue_SP(new (gEidosValuePool->AllocateChunk()) EidosValue_Object_singleton(this, gSLiM_InteractionType_Class))), - spatiality_string_(p_spatiality_string), reciprocal_(p_reciprocal), max_distance_(p_max_distance), max_distance_sq_(p_max_distance * p_max_distance), receiver_sex_(p_receiver_sex), exerter_sex_(p_exerter_sex), if_type_(IFType::kFixed), if_param1_(1.0), if_param2_(0.0), + spatiality_string_(p_spatiality_string), reciprocal_(p_reciprocal), max_distance_(p_max_distance), max_distance_sq_(p_max_distance * p_max_distance), receiver_sex_(p_receiver_sex), exerter_sex_(p_exerter_sex), if_type_(SpatialKernelType::kFixed), if_param1_(1.0), if_param2_(0.0), community_(p_community), interaction_type_id_(p_interaction_type_id) { // Figure out our spatiality, which is the number of spatial dimensions we actively use for distances @@ -763,23 +747,24 @@ double InteractionType::CalculateStrengthNoCallbacks(double p_distance) // It is the caller's responsibility to do that filtering, for performance reasons! // The caller is also responsible for guaranteeing that this is not a self-interaction, // and that it is not ruled out by sex-selectivity. + // SEE ALSO: Kernel::DensityForDistance(), which is parallel to this switch (if_type_) { - case IFType::kFixed: + case SpatialKernelType::kFixed: return (if_param1_); // fmax - case IFType::kLinear: + case SpatialKernelType::kLinear: return (if_param1_ * (1.0 - p_distance / max_distance_)); // fmax * (1 − d/dmax) - case IFType::kExponential: + case SpatialKernelType::kExponential: return (if_param1_ * exp(-if_param2_ * p_distance)); // fmax * exp(−λd) - case IFType::kNormal: + case SpatialKernelType::kNormal: return (if_param1_ * exp(-(p_distance * p_distance) / n_2param2sq_)); // fmax * exp(−d^2/2σ^2) - case IFType::kCauchy: + case SpatialKernelType::kCauchy: { double temp = p_distance / if_param2_; return (if_param1_ / (1.0 + temp * temp)); // fmax / (1+(d/λ)^2) } } - EIDOS_TERMINATION << "ERROR (InteractionType::CalculateStrengthNoCallbacks): (internal error) unexpected if_type_ value." << EidosTerminate(); + EIDOS_TERMINATION << "ERROR (InteractionType::CalculateStrengthNoCallbacks): (internal error) unexpected SpatialKernelType." << EidosTerminate(); } double InteractionType::CalculateStrengthWithCallbacks(double p_distance, Individual *p_receiver, Individual *p_exerter, std::vector &p_interaction_callbacks) @@ -2356,7 +2341,7 @@ void InteractionType::BuildSV_Distances_SS_3(SLiM_kdNode *root, double *nd, slim } } -// add neighbor strengths of type "f" (IFType::kFixed : fixed) to the sparse vector in 2D +// add neighbor strengths of type "f" (SpatialKernelType::kFixed : fixed) to the sparse vector in 2D void InteractionType::BuildSV_Strengths_f_2(SLiM_kdNode *root, double *nd, slim_popsize_t p_focal_individual_index, SparseVector *p_sparse_vector, int p_phase) { double d = dist_sq2(root, nd); @@ -2385,7 +2370,7 @@ void InteractionType::BuildSV_Strengths_f_2(SLiM_kdNode *root, double *nd, slim_ } } -// add neighbor strengths of type "l" (IFType::kLinear : linear) to the sparse vector in 2D +// add neighbor strengths of type "l" (SpatialKernelType::kLinear : linear) to the sparse vector in 2D void InteractionType::BuildSV_Strengths_l_2(SLiM_kdNode *root, double *nd, slim_popsize_t p_focal_individual_index, SparseVector *p_sparse_vector, int p_phase) { double d = dist_sq2(root, nd); @@ -2414,7 +2399,7 @@ void InteractionType::BuildSV_Strengths_l_2(SLiM_kdNode *root, double *nd, slim_ } } -// add neighbor strengths of type "e" (IFType::kExponential : exponential) to the sparse vector in 2D +// add neighbor strengths of type "e" (SpatialKernelType::kExponential : exponential) to the sparse vector in 2D void InteractionType::BuildSV_Strengths_e_2(SLiM_kdNode *root, double *nd, slim_popsize_t p_focal_individual_index, SparseVector *p_sparse_vector, int p_phase) { double d = dist_sq2(root, nd); @@ -2443,7 +2428,7 @@ void InteractionType::BuildSV_Strengths_e_2(SLiM_kdNode *root, double *nd, slim_ } } -// add neighbor strengths of type "n" (IFType::kNormal : normal/Gaussian) to the sparse vector in 2D +// add neighbor strengths of type "n" (SpatialKernelType::kNormal : normal/Gaussian) to the sparse vector in 2D void InteractionType::BuildSV_Strengths_n_2(SLiM_kdNode *root, double *nd, slim_popsize_t p_focal_individual_index, SparseVector *p_sparse_vector, int p_phase) { double d = dist_sq2(root, nd); @@ -2473,7 +2458,7 @@ void InteractionType::BuildSV_Strengths_n_2(SLiM_kdNode *root, double *nd, slim_ } -// add neighbor strengths of type "c" (IFType::kCauchy : Cauchy) to the sparse vector in 2D +// add neighbor strengths of type "c" (SpatialKernelType::kCauchy : Cauchy) to the sparse vector in 2D void InteractionType::BuildSV_Strengths_c_2(SLiM_kdNode *root, double *nd, slim_popsize_t p_focal_individual_index, SparseVector *p_sparse_vector, int p_phase) { double d = dist_sq2(root, nd); @@ -2736,13 +2721,13 @@ void InteractionType::FillSparseVectorForReceiverStrengths(SparseVector *sv, Ind switch (if_type_) { - case IFType::kFixed: BuildSV_Strengths_f_2(exerter_subpop_data.kd_root_, receiver_position, excluded_index, sv, 0); break; - case IFType::kLinear: BuildSV_Strengths_l_2(exerter_subpop_data.kd_root_, receiver_position, excluded_index, sv, 0); break; - case IFType::kExponential: BuildSV_Strengths_e_2(exerter_subpop_data.kd_root_, receiver_position, excluded_index, sv, 0); break; - case IFType::kNormal: BuildSV_Strengths_n_2(exerter_subpop_data.kd_root_, receiver_position, excluded_index, sv, 0); break; - case IFType::kCauchy: BuildSV_Strengths_c_2(exerter_subpop_data.kd_root_, receiver_position, excluded_index, sv, 0); break; + case SpatialKernelType::kFixed: BuildSV_Strengths_f_2(exerter_subpop_data.kd_root_, receiver_position, excluded_index, sv, 0); break; + case SpatialKernelType::kLinear: BuildSV_Strengths_l_2(exerter_subpop_data.kd_root_, receiver_position, excluded_index, sv, 0); break; + case SpatialKernelType::kExponential: BuildSV_Strengths_e_2(exerter_subpop_data.kd_root_, receiver_position, excluded_index, sv, 0); break; + case SpatialKernelType::kNormal: BuildSV_Strengths_n_2(exerter_subpop_data.kd_root_, receiver_position, excluded_index, sv, 0); break; + case SpatialKernelType::kCauchy: BuildSV_Strengths_c_2(exerter_subpop_data.kd_root_, receiver_position, excluded_index, sv, 0); break; default: - EIDOS_TERMINATION << "ERROR (InteractionType::FillSparseVectorForReceiverStrengths): (internal error) unoptimized IFType value." << EidosTerminate(); + EIDOS_TERMINATION << "ERROR (InteractionType::FillSparseVectorForReceiverStrengths): (internal error) unoptimized SpatialKernelType value." << EidosTerminate(); } sv->Finished(); @@ -2788,13 +2773,13 @@ void InteractionType::FillSparseVectorForReceiverStrengths(SparseVector *sv, Ind // CalculateStrengthNoCallbacks() is basically inlined here, moved outside the loop; see that function for comments switch (if_type_) { - case IFType::kFixed: + case SpatialKernelType::kFixed: { for (uint32_t col_iter = 0; col_iter < nnz; ++col_iter) values[col_iter] = (sv_value_t)if_param1_; break; } - case IFType::kLinear: + case SpatialKernelType::kLinear: { for (uint32_t col_iter = 0; col_iter < nnz; ++col_iter) { @@ -2804,7 +2789,7 @@ void InteractionType::FillSparseVectorForReceiverStrengths(SparseVector *sv, Ind } break; } - case IFType::kExponential: + case SpatialKernelType::kExponential: { for (uint32_t col_iter = 0; col_iter < nnz; ++col_iter) { @@ -2814,7 +2799,7 @@ void InteractionType::FillSparseVectorForReceiverStrengths(SparseVector *sv, Ind } break; } - case IFType::kNormal: + case SpatialKernelType::kNormal: { for (uint32_t col_iter = 0; col_iter < nnz; ++col_iter) { @@ -2824,7 +2809,7 @@ void InteractionType::FillSparseVectorForReceiverStrengths(SparseVector *sv, Ind } break; } - case IFType::kCauchy: + case SpatialKernelType::kCauchy: { for (uint32_t col_iter = 0; col_iter < nnz; ++col_iter) { @@ -2845,7 +2830,7 @@ void InteractionType::FillSparseVectorForReceiverStrengths(SparseVector *sv, Ind values[col_iter] = (sv_value_t)CalculateStrengthNoCallbacks(distance); } - EIDOS_TERMINATION << "ERROR (InteractionType::FillSparseVectorForReceiverStrengths): (internal error) unimplemented IFType case." << EidosTerminate(); + EIDOS_TERMINATION << "ERROR (InteractionType::FillSparseVectorForReceiverStrengths): (internal error) unimplemented SpatialKernelType case." << EidosTerminate(); } } } @@ -3450,7 +3435,7 @@ void InteractionType::SetProperty(EidosGlobalStringID p_property_id, const Eidos if (max_distance_ < 0.0) EIDOS_TERMINATION << "ERROR (InteractionType::SetProperty): the maximum interaction distance must be greater than or equal to zero." << EidosTerminate(); - if ((if_type_ == IFType::kLinear) && (std::isinf(max_distance_) || (max_distance_ <= 0.0))) + if ((if_type_ == SpatialKernelType::kLinear) && (std::isinf(max_distance_) || (max_distance_ <= 0.0))) EIDOS_TERMINATION << "ERROR (InteractionType::SetProperty): the maximum interaction distance must be finite and greater than zero when interaction type 'l' has been chosen." << EidosTerminate(); // tweak a flag to make SLiMgui update @@ -4150,7 +4135,7 @@ EidosValue_SP InteractionType::ExecuteMethod_drawByStrength(EidosGlobalStringID EIDOS_TERMINATION << "ERROR (InteractionType::ExecuteMethod_drawByStrength): drawByStrength() requires count >= 0." << EidosTerminate(); bool has_interaction_callbacks = (exerter_subpop_data.evaluation_interaction_callbacks_.size() != 0); - bool optimize_fixed_interaction_strengths = (!has_interaction_callbacks && (if_type_ == IFType::kFixed)); + bool optimize_fixed_interaction_strengths = (!has_interaction_callbacks && (if_type_ == SpatialKernelType::kFixed)); if (!returnDict) { @@ -4189,7 +4174,7 @@ EidosValue_SP InteractionType::ExecuteMethod_drawByStrength(EidosGlobalStringID if (exerter_index_in_subpop != receiver_index) if ((exerter_sex_ == IndividualSex::kUnspecified) || (exerter_sex_ == exerter->sex_)) - strength = ApplyInteractionCallbacks(receiver, exerter, if_param1_, NAN, callbacks); // hard-coding interaction function "f" (IFType::kFixed), which is required + strength = ApplyInteractionCallbacks(receiver, exerter, if_param1_, NAN, callbacks); // hard-coding interaction function "f" (SpatialKernelType::kFixed), which is required total_interaction_strength += strength; cached_strength.emplace_back(strength); @@ -4674,7 +4659,7 @@ EidosValue_SP InteractionType::ExecuteMethod_localPopulationDensity(EidosGlobalS EidosValue *clipped_integrals = clipped_integrals_SP.get(); // Decide whether we can use our optimized case below - bool optimize_fixed_interaction_strengths = (!has_interaction_callbacks && (if_type_ == IFType::kFixed)); + bool optimize_fixed_interaction_strengths = (!has_interaction_callbacks && (if_type_ == SpatialKernelType::kFixed)); if (receivers_count == 1) { @@ -5607,96 +5592,19 @@ EidosValue_SP InteractionType::ExecuteMethod_neighborCountOfPoint(EidosGlobalStr EidosValue_SP InteractionType::ExecuteMethod_setInteractionFunction(EidosGlobalStringID p_method_id, const std::vector &p_arguments, EidosInterpreter &p_interpreter) { #pragma unused (p_method_id, p_arguments, p_interpreter) - EidosValue_String *functionType_value = (EidosValue_String *)p_arguments[0].get(); - - const std::string &if_type_string = functionType_value->StringRefAtIndex(0, nullptr); - IFType if_type; - int expected_if_param_count = 0; - std::vector if_parameters; - std::vector if_strings; if (AnyEvaluated()) EIDOS_TERMINATION << "ERROR (InteractionType::ExecuteMethod_setInteractionFunction): setInteractionFunction() cannot be called while the interaction is being evaluated; call unevaluate() first, or call setInteractionFunction() prior to evaluation of the interaction." << EidosTerminate(); - if (if_type_string.compare(gStr_f) == 0) - { - if_type = IFType::kFixed; - expected_if_param_count = 1; - } - else if (if_type_string.compare(gStr_l) == 0) - { - if_type = IFType::kLinear; - expected_if_param_count = 1; - - if (std::isinf(max_distance_) || (max_distance_ <= 0.0)) - EIDOS_TERMINATION << "ERROR (InteractionType::ExecuteMethod_setInteractionFunction): interaction type 'l' cannot be set in setInteractionFunction() unless a finite maximum interaction distance greater than zero has been set." << EidosTerminate(); - } - else if (if_type_string.compare(gStr_e) == 0) - { - if_type = IFType::kExponential; - expected_if_param_count = 2; - } - else if (if_type_string.compare(gEidosStr_n) == 0) - { - if_type = IFType::kNormal; - expected_if_param_count = 2; - } - else if (if_type_string.compare(gEidosStr_c) == 0) - { - if_type = IFType::kCauchy; - expected_if_param_count = 2; - } - else - EIDOS_TERMINATION << "ERROR (InteractionType::ExecuteMethod_setInteractionFunction): setInteractionFunction() functionType \"" << if_type_string << "\" must be \"f\", \"l\", \"e\", \"n\", or \"c\"." << EidosTerminate(); - - if ((spatiality_ == 0) && (if_type != IFType::kFixed)) - EIDOS_TERMINATION << "ERROR (InteractionType::ExecuteMethod_setInteractionFunction): setInteractionFunction() requires functionType 'f' for non-spatial interactions." << EidosTerminate(); - - if ((int)p_arguments.size() != 1 + expected_if_param_count) - EIDOS_TERMINATION << "ERROR (InteractionType::ExecuteMethod_setInteractionFunction): setInteractionFunction() functionType \"" << if_type << "\" requires exactly " << expected_if_param_count << " DFE parameter" << (expected_if_param_count == 1 ? "" : "s") << "." << EidosTerminate(); - - for (int if_param_index = 0; if_param_index < expected_if_param_count; ++if_param_index) - { - EidosValue *if_param_value = p_arguments[1 + if_param_index].get(); - EidosValueType if_param_type = if_param_value->Type(); - - if ((if_param_type != EidosValueType::kValueFloat) && (if_param_type != EidosValueType::kValueInt)) - EIDOS_TERMINATION << "ERROR (InteractionType::ExecuteMethod_setInteractionFunction): setInteractionFunction() requires that the parameters for this interaction function be of type numeric (integer or float)." << EidosTerminate(); - - if_parameters.emplace_back(if_param_value->FloatAtIndex(0, nullptr)); - } - - // Bounds-check the IF parameters in the cases where there is a hard bound - switch (if_type) - { - case IFType::kFixed: - // no limits on fixed IFs; doesn't make much sense to use 0.0, but it's not illegal - break; - case IFType::kLinear: - // no limits on linear IFs; doesn't make much sense to use 0.0, but it's not illegal - break; - case IFType::kExponential: - // no limits on exponential IFs; a shape of 0.0 doesn't make much sense, but it's not illegal - break; - case IFType::kNormal: - // no limits on the maximum strength (although 0.0 doesn't make much sense); sd must be >= 0 - if (if_parameters[1] < 0.0) - EIDOS_TERMINATION << "ERROR (InteractionType::ExecuteMethod_setInteractionFunction): an interaction function of type \"n\" must have a standard deviation parameter >= 0." << EidosTerminate(); - break; - case IFType::kCauchy: - // no limits on the maximum strength (although 0.0 doesn't make much sense); scale must be > 0 - if (if_parameters[1] <= 0.0) - EIDOS_TERMINATION << "ERROR (InteractionType::ExecuteMethod_setInteractionFunction): an interaction function of type \"c\" must have a scale parameter > 0." << EidosTerminate(); - break; - } + // SpatialKernel parses and bounds-checks our arguments for us + SpatialKernel kernel(spatiality_, max_distance_, p_arguments, 0); // Everything seems to be in order, so replace our IF info with the new info - if_type_ = if_type; - if_param1_ = ((if_parameters.size() >= 1) ? if_parameters[0] : 0.0); - if_param2_ = ((if_parameters.size() >= 2) ? if_parameters[1] : 0.0); - - if (if_type_ == IFType::kNormal) - n_2param2sq_ = 2.0 * if_param2_ * if_param2_; + // FIXME we could consider actually keeping an internal SpatialKernel instance permanently + if_type_ = kernel.kernel_type_; + if_param1_ = kernel.kernel_param1_; + if_param2_ = kernel.kernel_param2_; + n_2param2sq_ = kernel.n_2param2sq_; // mark that interaction types changed, so they get redisplayed in SLiMgui community_.interaction_types_changed_ = true; @@ -5848,7 +5756,7 @@ EidosValue_SP InteractionType::ExecuteMethod_strength(EidosGlobalStringID p_meth Individual *exerter = exerter_subpop->parent_individuals_[exerter_index]; if ((exerter_sex_ == IndividualSex::kUnspecified) || (exerter_sex_ == exerter->sex_)) - strength = ApplyInteractionCallbacks(receiver, exerter, if_param1_, NAN, callbacks); // hard-coding interaction function "f" (IFType::kFixed), which is required + strength = ApplyInteractionCallbacks(receiver, exerter, if_param1_, NAN, callbacks); // hard-coding interaction function "f" (SpatialKernelType::kFixed), which is required } result_vec->set_float_no_check(strength, exerter_index); @@ -5877,7 +5785,7 @@ EidosValue_SP InteractionType::ExecuteMethod_strength(EidosGlobalStringID p_meth if (exerter_index_in_subpop != receiver_index) if ((exerter_sex_ == IndividualSex::kUnspecified) || (exerter_sex_ == exerter->sex_)) - strength = ApplyInteractionCallbacks(receiver, exerter, if_param1_, NAN, callbacks); // hard-coding interaction function "f" (IFType::kFixed), which is required + strength = ApplyInteractionCallbacks(receiver, exerter, if_param1_, NAN, callbacks); // hard-coding interaction function "f" (SpatialKernelType::kFixed), which is required result_vec->set_float_no_check(strength, exerter_index); } diff --git a/core/interaction_type.h b/core/interaction_type.h index ad2e3ad66..0a81903f0 100644 --- a/core/interaction_type.h +++ b/core/interaction_type.h @@ -40,6 +40,7 @@ #include "slim_eidos_block.h" #include "sparse_vector.h" #include "subpopulation.h" +#include "spatial_kernel.h" class Species; @@ -51,19 +52,6 @@ class InteractionType_Class; extern EidosClass *gSLiM_InteractionType_Class; -// This enumeration represents a type of interaction function (IF) that an -// interaction type can use to convert distances to interaction strengths -enum class IFType : char { - kFixed = 0, - kLinear, - kExponential, - kNormal, - kCauchy -}; - -std::ostream& operator<<(std::ostream& p_out, IFType p_if_type); - - // This class uses an internal implementation of kd-trees for fast nearest-neighbor finding. We use the same data structure to // save computed distances and interaction strengths. A value of NaN is used as a placeholder to indicate that a given value // has not yet been calculated, and we fill the data structure in lazily. We keep one such data structure per evaluated @@ -149,7 +137,7 @@ class InteractionType : public EidosDictionaryUnretained slim_usertag_t tag_value_ = SLIM_TAG_UNSET_VALUE; // a user-defined tag value - IFType if_type_; // the interaction function (IF) to use + SpatialKernelType if_type_; // the interaction function (IF) to use double if_param1_, if_param2_; // the parameters for that IF (not all of which may be used) double n_2param2sq_; // for type "n", precalculated == 2.0 * if_param2_ * if_param2_ diff --git a/core/slim_globals.cpp b/core/slim_globals.cpp index 22ab94fe4..03a282beb 100644 --- a/core/slim_globals.cpp +++ b/core/slim_globals.cpp @@ -1368,6 +1368,7 @@ const std::string &gStr_interpolateValues = EidosRegisteredString("interpolateVa const std::string &gStr_mapColor = EidosRegisteredString("mapColor", gID_mapColor); const std::string &gStr_mapImage = EidosRegisteredString("mapImage", gID_mapImage); const std::string &gStr_mapValue = EidosRegisteredString("mapValue", gID_mapValue); +const std::string &gStr_smoothValues = EidosRegisteredString("smoothValues", gID_smoothValues); const std::string &gStr_outputMSSample = EidosRegisteredString("outputMSSample", gID_outputMSSample); const std::string &gStr_outputVCFSample = EidosRegisteredString("outputVCFSample", gID_outputVCFSample); const std::string &gStr_outputSample = EidosRegisteredString("outputSample", gID_outputSample); diff --git a/core/slim_globals.h b/core/slim_globals.h index 6a6b3432f..f4916387f 100644 --- a/core/slim_globals.h +++ b/core/slim_globals.h @@ -956,6 +956,7 @@ extern const std::string &gStr_interpolateValues; extern const std::string &gStr_mapColor; extern const std::string &gStr_mapImage; extern const std::string &gStr_mapValue; +extern const std::string &gStr_smoothValues; extern const std::string &gStr_outputMSSample; extern const std::string &gStr_outputVCFSample; extern const std::string &gStr_outputSample; @@ -1339,6 +1340,7 @@ enum _SLiMGlobalStringID : int { gID_mapColor, gID_mapImage, gID_mapValue, + gID_smoothValues, gID_outputMSSample, gID_outputVCFSample, gID_outputSample, diff --git a/core/spatial_kernel.cpp b/core/spatial_kernel.cpp new file mode 100644 index 000000000..a0d0fedc5 --- /dev/null +++ b/core/spatial_kernel.cpp @@ -0,0 +1,368 @@ +// +// spatial_kernel.cpp +// SLiM +// +// Created by Ben Haller on 9/9/23. +// Copyright (c) 2023 Philipp Messer. All rights reserved. +// A product of the Messer Lab, http://messerlab.org/slim/ +// + +// This file is part of SLiM. +// +// SLiM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +// +// SLiM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with SLiM. If not, see . + + +#include "spatial_kernel.h" +#include "eidos_value.h" + + +// stream output for enumerations +std::ostream& operator<<(std::ostream& p_out, SpatialKernelType p_kernel_type) +{ + switch (p_kernel_type) + { + case SpatialKernelType::kFixed: p_out << gStr_f; break; + case SpatialKernelType::kLinear: p_out << gStr_l; break; + case SpatialKernelType::kExponential: p_out << gStr_e; break; + case SpatialKernelType::kNormal: p_out << gEidosStr_n; break; + case SpatialKernelType::kCauchy: p_out << gEidosStr_c; break; + } + + return p_out; +} + + +#pragma mark - +#pragma mark SpatialKernel +#pragma mark - + +SpatialKernel::SpatialKernel(int p_dimensionality, double p_maxDistance, const std::vector &p_arguments, int p_first_kernel_arg) : dimensionality_(p_dimensionality), max_distance_(p_maxDistance) +{ + // This constructs a kernel from the arguments given, beginning at argument p_first_kernel_arg. + // For example, take the smoothValues() method of SpatialKernel: + // + // - (void)smoothValues(float$ maxDistance, string$ functionType, ...) + // + // It parses out maxDistance and passes it to us; it then forwards its remaining + // arguments, with p_first_kernel_arg == 1, to define the shape of the kernel it wants. + // The ellipsis arguments are as they are for setInteractionFunction(); this class is + // basically a grid-sampled version of the same style of kernel that InteractionType + // uses, and indeed, InteractionType now uses SpatialKernel for some of its work. + // + // The grid sampling is based upon the spatial scale established by a given SpatialMap; + // the max distance and other kernel parameters are in terms of that scale. + + if ((p_dimensionality < 0) || (p_dimensionality > 3)) + EIDOS_TERMINATION << "ERROR (SpatialKernel::SpatialKernel): spatial kernel dimensionality must be 0, 1, 2, or 3." << EidosTerminate(); + if (max_distance_ <= 0) + EIDOS_TERMINATION << "ERROR (SpatialKernel::SpatialKernel): spatial kernel maxDistance must be greater than zero." << EidosTerminate(); + + // Parse the arguments that define our kernel shape + if (p_arguments[p_first_kernel_arg]->Type() != EidosValueType::kValueString) + EIDOS_TERMINATION << "ERROR (SpatialKernel::SpatialKernel): (internal error) functionType is not a string." << EidosTerminate(); + + EidosValue_String *functionType_value = (EidosValue_String *)p_arguments[p_first_kernel_arg].get(); + + const std::string &k_type_string = functionType_value->StringRefAtIndex(0, nullptr); + SpatialKernelType k_type; + int expected_k_param_count = 0; + std::vector k_parameters; + + if (k_type_string.compare(gStr_f) == 0) + { + k_type = SpatialKernelType::kFixed; + expected_k_param_count = 1; + } + else if (k_type_string.compare(gStr_l) == 0) + { + k_type = SpatialKernelType::kLinear; + expected_k_param_count = 1; + + if (std::isinf(max_distance_) || (max_distance_ <= 0.0)) + EIDOS_TERMINATION << "ERROR (SpatialKernel::SpatialKernel): spatial kernel type 'l' cannot be used unless a finite maximum interaction distance greater than zero has been set." << EidosTerminate(); + } + else if (k_type_string.compare(gStr_e) == 0) + { + k_type = SpatialKernelType::kExponential; + expected_k_param_count = 2; + } + else if (k_type_string.compare(gEidosStr_n) == 0) + { + k_type = SpatialKernelType::kNormal; + expected_k_param_count = 2; + } + else if (k_type_string.compare(gEidosStr_c) == 0) + { + k_type = SpatialKernelType::kCauchy; + expected_k_param_count = 2; + } + else + EIDOS_TERMINATION << "ERROR (SpatialKernel::SpatialKernel): spatial kernel functionType \"" << k_type_string << "\" must be \"f\", \"l\", \"e\", \"n\", or \"c\"." << EidosTerminate(); + + if ((dimensionality_ == 0) && (k_type != SpatialKernelType::kFixed)) + EIDOS_TERMINATION << "ERROR (SpatialKernel::SpatialKernel): spatial kernel functionType 'f' is required for non-spatial interactions." << EidosTerminate(); + + if ((int)p_arguments.size() - p_first_kernel_arg != 1 + expected_k_param_count) + EIDOS_TERMINATION << "ERROR (SpatialKernel::SpatialKernel): spatial kernel functionType \"" << k_type << "\" requires exactly " << expected_k_param_count << " kernel configuration parameter" << (expected_k_param_count == 1 ? "" : "s") << "." << EidosTerminate(); + + for (int k_param_index = 0; k_param_index < expected_k_param_count; ++k_param_index) + { + EidosValue *k_param_value = p_arguments[1 + k_param_index + p_first_kernel_arg].get(); + EidosValueType k_param_type = k_param_value->Type(); + + if ((k_param_type != EidosValueType::kValueFloat) && (k_param_type != EidosValueType::kValueInt)) + EIDOS_TERMINATION << "ERROR (SpatialKernel::SpatialKernel): the parameters for this spatial kernel type must be numeric (integer or float)." << EidosTerminate(); + + k_parameters.emplace_back(k_param_value->FloatAtIndex(0, nullptr)); + } + + // Bounds-check the IF parameters in the cases where there is a hard bound + switch (k_type) + { + case SpatialKernelType::kFixed: + // no limits on fixed IFs; doesn't make much sense to use 0.0, but it's not illegal + break; + case SpatialKernelType::kLinear: + // no limits on linear IFs; doesn't make much sense to use 0.0, but it's not illegal + break; + case SpatialKernelType::kExponential: + // no limits on exponential IFs; a shape of 0.0 doesn't make much sense, but it's not illegal + break; + case SpatialKernelType::kNormal: + // no limits on the maximum strength (although 0.0 doesn't make much sense); sd must be >= 0 + if (k_parameters[1] < 0.0) + EIDOS_TERMINATION << "ERROR (SpatialKernel::SpatialKernel): spatial kernel type \"n\" must have a standard deviation parameter >= 0." << EidosTerminate(); + break; + case SpatialKernelType::kCauchy: + // no limits on the maximum strength (although 0.0 doesn't make much sense); scale must be > 0 + if (k_parameters[1] <= 0.0) + EIDOS_TERMINATION << "ERROR (SpatialKernel::SpatialKernel): spatial kernel type \"c\" must have a scale parameter > 0." << EidosTerminate(); + break; + } + + // Everything seems to be in order, so replace our kernel info with the new info + kernel_type_ = k_type; + kernel_param1_ = ((k_parameters.size() >= 1) ? k_parameters[0] : 0.0); + kernel_param2_ = ((k_parameters.size() >= 2) ? k_parameters[1] : 0.0); + n_2param2sq_ = (kernel_type_ == SpatialKernelType::kNormal) ? (2.0 * kernel_param2_ * kernel_param2_) : 0.0; +} + +SpatialKernel::~SpatialKernel(void) +{ + if (values_) + free(values_); +} + +void SpatialKernel::CalculateGridValues(SpatialMap &p_map) +{ + if ((dimensionality_ < 1) || (dimensionality_ > 3)) + EIDOS_TERMINATION << "ERROR (SpatialKernel::CalculateGridValues): grid values can only be calculated for kernels with dimensionality of 1, 2, or 3." << EidosTerminate(nullptr); + if ((max_distance_ <= 0.0) || (!std::isfinite(max_distance_))) + EIDOS_TERMINATION << "ERROR (SpatialKernel::CalculateGridValues): grid values can only be calculated for kernels with a maxDistance that is positive and finite." << EidosTerminate(nullptr); + + // Derive our spatial scale from the given spatial map, which provides a correspondence between + // spatial bounds and pixel sizes; after this, we do not use the SpatialMap, so these scales + // could instead be passed in. + double spatial_size_a = (p_map.bounds_a1_ - p_map.bounds_a0_); + double spatial_size_b = (dimensionality_ >= 2) ? (p_map.bounds_b1_ - p_map.bounds_b0_) : 0.0; + double spatial_size_c = (dimensionality_ >= 3) ? (p_map.bounds_c1_ - p_map.bounds_c0_) : 0.0; + + pixels_to_spatial_a_ = (spatial_size_a / (p_map.grid_size_[0] - 1)); + pixels_to_spatial_b_ = (dimensionality_ >= 2) ? (spatial_size_b / (p_map.grid_size_[1] - 1)) : 0.0; + pixels_to_spatial_c_ = (dimensionality_ >= 3) ? (spatial_size_c / (p_map.grid_size_[2] - 1)) : 0.0; + + dim[0] = 0; + dim[1] = 0; + dim[2] = 0; + + int64_t values_size = 1; + + double pixelsize_d = (max_distance_ * 2) / pixels_to_spatial_a_; // convert spatial to pixels + int64_t pixelsize = (int64_t)ceil(pixelsize_d); + if (pixelsize % 2 == 0) // round up to an odd integer + pixelsize++; + dim[0] = pixelsize; + values_size *= pixelsize; + + if (dimensionality_ >= 2) + { + pixelsize_d = (max_distance_ * 2) / pixels_to_spatial_b_; // convert spatial to pixels + pixelsize = (int64_t)ceil(pixelsize_d); + if (pixelsize % 2 == 0) // round up to an odd integer + pixelsize++; + dim[1] = pixelsize; + values_size *= pixelsize; + } + + if (dimensionality_ >= 3) + { + pixelsize_d = (max_distance_ * 2) / pixels_to_spatial_c_; // convert spatial to pixels + pixelsize = (int64_t)ceil(pixelsize_d); + if (pixelsize % 2 == 0) // round up to an odd integer + pixelsize++; + dim[2] = pixelsize; + values_size *= pixelsize; + } + + // Allocate our values buffer + values_ = (double *)malloc(values_size * sizeof(double)); + + if (!values_) + EIDOS_TERMINATION << "ERROR (SpatialKernel::CalculateGridValues): allocation failed; you may need to raise the memory limit for SLiM." << EidosTerminate(nullptr); + + // Set our values + switch (dimensionality_) + { + case 1: + { + int64_t kernel_offset_a = dim[0] / 2; // rounds down + + for (int64_t a = 0; a < dim[0]; ++a) + { + double distance = (a - kernel_offset_a) * pixels_to_spatial_a_; + double density = (distance > max_distance_) ? 0.0 : DensityForDistance(distance); + + values_[a] = density; + } + break; + } + case 2: + { + int64_t kernel_offset_a = dim[0] / 2; // rounds down + int64_t kernel_offset_b = dim[1] / 2; // rounds down + + for (int64_t a = 0; a < dim[0]; ++a) + { + double dist_a = (a - kernel_offset_a) * pixels_to_spatial_a_; + double dist_a_sq = dist_a * dist_a; + + for (int64_t b = 0; b < dim[1]; ++b) + { + double dist_b = (b - kernel_offset_b) * pixels_to_spatial_b_; + double dist_b_sq = dist_b * dist_b; + double distance = sqrt(dist_a_sq + dist_b_sq); + + double density = (distance > max_distance_) ? 0.0 : DensityForDistance(distance); + + values_[a + b * dim[0]] = density; + } + } + break; + } + case 3: + { + int64_t kernel_offset_a = dim[0] / 2; // rounds down + int64_t kernel_offset_b = dim[1] / 2; // rounds down + int64_t kernel_offset_c = dim[2] / 2; // rounds down + + for (int64_t a = 0; a < dim[0]; ++a) + { + double dist_a = (a - kernel_offset_a) * pixels_to_spatial_a_; + double dist_a_sq = dist_a * dist_a; + + for (int64_t b = 0; b < dim[1]; ++b) + { + double dist_b = (b - kernel_offset_b) * pixels_to_spatial_b_; + double dist_b_sq = dist_b * dist_b; + + for (int64_t c = 0; c < dim[2]; ++c) + { + double dist_c = (c - kernel_offset_c) * pixels_to_spatial_c_; + double dist_c_sq = dist_c * dist_c; + double distance = sqrt(dist_a_sq + dist_b_sq + dist_c_sq); + double density = (distance > max_distance_) ? 0.0 : DensityForDistance(distance); + + values_[a] = density; + } + } + } + break; + } + default: break; + } +} + +double SpatialKernel::DensityForDistance(double p_distance) +{ + // SEE ALSO: InteractionType::CalculateStrengthNoCallbacks(), which is parallel to this + switch (kernel_type_) + { + case SpatialKernelType::kFixed: + return (kernel_param1_); // fmax + case SpatialKernelType::kLinear: + return (kernel_param1_ * (1.0 - p_distance / max_distance_)); // fmax * (1 − d/dmax) + case SpatialKernelType::kExponential: + return (kernel_param1_ * exp(-kernel_param2_ * p_distance)); // fmax * exp(−λd) + case SpatialKernelType::kNormal: + return (kernel_param1_ * exp(-(p_distance * p_distance) / n_2param2sq_)); // fmax * exp(−d^2/2σ^2) + case SpatialKernelType::kCauchy: + { + double temp = p_distance / kernel_param2_; + return (kernel_param1_ / (1.0 + temp * temp)); // fmax / (1+(d/λ)^2) + } + } + EIDOS_TERMINATION << "ERROR (SpatialKernel::DensityForDistance): (internal error) unexpected SpatialKernelType value." << EidosTerminate(); +} + +std::ostream& operator<<(std::ostream& p_out, SpatialKernel &p_kernel) +{ + std::cout << "Kernel with dimensionality_ == " << p_kernel.dimensionality_ << ":" << std::endl; + std::cout << " max_distance_ == " << p_kernel.max_distance_ << std::endl; + std::cout << " kernel_type_ == \"" << p_kernel.kernel_type_ << "\"" << std::endl; + std::cout << " kernel_param1_ == " << p_kernel.kernel_param1_ << std::endl; + std::cout << " kernel_param2_ == " << p_kernel.kernel_param2_ << std::endl; + std::cout << " n_2param2sq_ == " << p_kernel.n_2param2sq_ << std::endl; + std::cout << " dim[3] == {" << p_kernel.dim[0] << ", " << p_kernel.dim[1] << ", " << p_kernel.dim[2] << "}" << std::endl; + + if (p_kernel.values_) + { + std::cout << " pixels_to_spatial_a_ == " << p_kernel.pixels_to_spatial_a_ << std::endl; + std::cout << " pixels_to_spatial_b_ == " << p_kernel.pixels_to_spatial_b_ << std::endl; + std::cout << " pixels_to_spatial_c_ == " << p_kernel.pixels_to_spatial_c_ << std::endl; + } + + std::cout << " values =="; + + switch (p_kernel.dimensionality_) + { + case 1: + break; + case 2: + for (int b = 0; b < p_kernel.dim[1]; b++) + { + std::cout << std::endl << " "; + + for (int a = 0; a < p_kernel.dim[0]; a++) + { + std::ostringstream os; + + os.precision(3); + os << std::fixed; + os << p_kernel.values_[a + b * p_kernel.dim[0]]; + + std::cout << os.str() << " "; + } + } + break; + case 3: + break; + default: break; + } + std::cout << std::endl; + + return p_out; +} + + + + + + + + diff --git a/core/spatial_kernel.h b/core/spatial_kernel.h new file mode 100644 index 000000000..e404903cd --- /dev/null +++ b/core/spatial_kernel.h @@ -0,0 +1,90 @@ +// +// spatial_kernel.h +// SLiM +// +// Created by Ben Haller on 9/9/23. +// Copyright (c) 2023 Philipp Messer. All rights reserved. +// A product of the Messer Lab, http://messerlab.org/slim/ +// + +// This file is part of SLiM. +// +// SLiM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +// +// SLiM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along with SLiM. If not, see . + +/* + + The class SpatialKernel represents a spatial kernel of some shape. It is used by both InteractionType and SpatialMap + to represent the kernels they use internally. It is not visible in Eidos, at least for now. + + */ + +#ifndef __SLiM__spatial_kernel__ +#define __SLiM__spatial_kernel__ + +#include +#include + +#include "eidos_value.h" +#include "spatial_map.h" + + +// This enumeration represents a type of interaction function (IF) that an +// interaction type can use to convert distances to interaction strengths +enum class SpatialKernelType : char { + kFixed = 0, + kLinear, + kExponential, + kNormal, + kCauchy +}; + +std::ostream& operator<<(std::ostream& p_out, SpatialKernelType p_kernel_type); + + +#pragma mark - +#pragma mark SpatialKernel +#pragma mark - + +class SpatialKernel +{ +public: + // core kernel definition + int dimensionality_; // 1, 2, or 3: how many dimensions the kernel data is + double max_distance_; // the maximum spatial distance out to which the kernel stretches + double pixels_to_spatial_a_; // multiply by this to convert pixels to spatial scale for a + double pixels_to_spatial_b_; // multiply by this to convert pixels to spatial scale for b + double pixels_to_spatial_c_; // multiply by this to convert pixels to spatial scale for c + + SpatialKernelType kernel_type_; // the kernel type to use + double kernel_param1_, kernel_param2_; // parameters for that kernel type (not all of which may be used) + double n_2param2sq_; // for type "n", precalc 2.0 * kernel_param2_ * kernel_param2_ + + // discrete grid values; set up only if CalculateGridValues() is called + double *values_ = nullptr; // raw kernel pixel data, malloced + int64_t dim[3] = {0, 0, 0}; // pixel dimensions of values_ for 1, 2, or 3 axes + +public: + SpatialKernel(const SpatialKernel&) = delete; // no copying + SpatialKernel& operator=(const SpatialKernel&) = delete; // no copying + SpatialKernel(void) = delete; // no null construction + SpatialKernel(int p_dimensionality, double p_maxDistance, const std::vector &p_arguments, int p_first_kernel_arg); + ~SpatialKernel(void); + + void CalculateGridValues(SpatialMap &p_map); + double DensityForDistance(double p_distance); + + friend void SpatialMap::Convolve_S1(SpatialKernel &kernel); + friend void SpatialMap::Convolve_S2(SpatialKernel &kernel); + friend void SpatialMap::Convolve_S3(SpatialKernel &kernel); +}; + +std::ostream& operator<<(std::ostream& p_out, SpatialKernel &p_kernel); + + +#endif /* __SLiM__spatial_kernel__ */ diff --git a/core/spatial_map.cpp b/core/spatial_map.cpp index 84b4537f6..7754dbee1 100644 --- a/core/spatial_map.cpp +++ b/core/spatial_map.cpp @@ -19,6 +19,7 @@ #include "spatial_map.h" +#include "spatial_kernel.h" #include "subpopulation.h" #include "eidos_class_Image.h" @@ -28,7 +29,7 @@ #pragma mark - -#pragma mark _SpatialMap +#pragma mark SpatialMap #pragma mark - SpatialMap::SpatialMap(std::string p_name, std::string p_spatiality_string, Subpopulation *p_subpop, EidosValue *p_values, bool p_interpolate, EidosValue *p_value_range, EidosValue *p_colors) : @@ -284,14 +285,10 @@ void SpatialMap::TakeValuesFromEidosValue(EidosValue *p_values, std::string p_co if (values_dimcount != spatiality_) EIDOS_TERMINATION << "ERROR (" << p_code_name << "): " << p_code_name << " the dimensionality of the supplied vector/matrix/array does not match the spatiality defined for the map." << EidosTerminate(); + int dimension_index; values_size_ = 1; - // There is no longer a gridSize parameter (as of SLiM 3.5), so values must be a matrix/array that matches the spatiality of the map - grid_size_[0] = 0; - grid_size_[1] = 0; - grid_size_[2] = 0; - - for (int dimension_index = 0; dimension_index < spatiality_; ++dimension_index) + for (dimension_index = 0; dimension_index < spatiality_; ++dimension_index) { int64_t dimension_size = (values_dimcount == 1) ? p_values->Count() : values_dim[dimension_index]; // treat a vector as a 1D matrix @@ -301,6 +298,8 @@ void SpatialMap::TakeValuesFromEidosValue(EidosValue *p_values, std::string p_co grid_size_[dimension_index] = dimension_size; values_size_ *= dimension_size; } + for ( ; dimension_index < 3; ++dimension_index) + grid_size_[dimension_index] = 0; // Matrices and arrays use dim[0] as the number of rows, and dim[1] as the number of cols; spatial maps do the opposite, // following standard image conventions (by row, not by column); we therefore need to swap grid_size_[0] and grid_size_[1] @@ -369,14 +368,10 @@ void SpatialMap::TakeOverMallocedValues(double *p_values, int64_t p_dimcount, in if (p_dimcount != spatiality_) EIDOS_TERMINATION << "ERROR (SpatialMap::TakeOverMallocedValues): (internal error) the dimensionality of the supplied values does not match the spatiality defined for the map." << EidosTerminate(); + int dimension_index; values_size_ = 1; - // There is no longer a gridSize parameter (as of SLiM 3.5), so values must be a matrix/array that matches the spatiality of the map - grid_size_[0] = 0; - grid_size_[1] = 0; - grid_size_[2] = 0; - - for (int dimension_index = 0; dimension_index < spatiality_; ++dimension_index) + for (dimension_index = 0; dimension_index < spatiality_; ++dimension_index) { int64_t dimension_size = p_dimensions[dimension_index]; @@ -386,6 +381,8 @@ void SpatialMap::TakeOverMallocedValues(double *p_values, int64_t p_dimcount, in grid_size_[dimension_index] = dimension_size; values_size_ *= dimension_size; } + for ( ; dimension_index < 3; ++dimension_index) + grid_size_[dimension_index] = 0; // Take over the passed buffer free(values_); @@ -655,6 +652,210 @@ void SpatialMap::ColorForValue(double p_value, float *p_rgb_ptr) } } +void SpatialMap::Convolve_S1(SpatialKernel &kernel) +{ + if (spatiality_ != 1) + EIDOS_TERMINATION << "ERROR (SpatialMap::Convolve_S1): (internal error) map spatiality 1 required." << EidosTerminate(); + if (kernel.dimensionality_ != 1) + EIDOS_TERMINATION << "ERROR (SpatialMap::Convolve_S1): (internal error) kernel dimensionality 1 required." << EidosTerminate(); + + int64_t kernel_dim_a = kernel.dim[0]; + + if ((kernel_dim_a < 1) || (kernel_dim_a % 2 == 0)) + EIDOS_TERMINATION << "ERROR (SpatialMap::Convolve_S1): (internal error) kernel dimensions must be odd." << EidosTerminate(); + + int64_t dim_a = grid_size_[0]; + double *new_values = (double *)malloc(dim_a * sizeof(double)); + + if (!new_values) + EIDOS_TERMINATION << "ERROR (SpatialMap::Convolve_S1): allocation failed; you may need to raise the memory limit for SLiM." << EidosTerminate(nullptr); + + // this assumes the kernel's dimensions are symmetrical around its center, and relies on rounding (which is guaranteed) + int64_t kernel_a_offset = -(kernel_dim_a / 2); + double *kernel_values = kernel.values_; + double *new_values_ptr = new_values; + + for (int64_t a = 0; a < dim_a; ++a) + { + // calculate the kernel's effect at point (a) + double kernel_total = 0.0; + double conv_total = 0.0; + + for (int64_t kernel_a = 0; kernel_a < kernel_dim_a; kernel_a++) + { + int64_t conv_a = a + kernel_a + kernel_a_offset; + + // clip to bounds + if ((conv_a < 0) || (conv_a >= dim_a)) + continue; + + // this point is within bounds; add it in to the totals + double kernel_value = kernel_values[kernel_a]; + double pixel_value = values_[conv_a]; + + // we keep a total of the kernel values that were within bounds, for this point + kernel_total += kernel_value; + + // and we keep a total of the convolution - kernel values times pixel values + conv_total += kernel_value * pixel_value; + } + + *(new_values_ptr++) = ((kernel_total > 0) ? (conv_total / kernel_total) : 0); + } + + TakeOverMallocedValues(new_values, 1, grid_size_); // takes new_values from us +} + +void SpatialMap::Convolve_S2(SpatialKernel &kernel) +{ + if (spatiality_ != 2) + EIDOS_TERMINATION << "ERROR (SpatialMap::Convolve_S2): (internal error) map spatiality 2 required." << EidosTerminate(); + if (kernel.dimensionality_ != 2) + EIDOS_TERMINATION << "ERROR (SpatialMap::Convolve_S2): (internal error) kernel dimensionality 2 required." << EidosTerminate(); + + int64_t kernel_dim_a = kernel.dim[0]; + int64_t kernel_dim_b = kernel.dim[1]; + + if ((kernel_dim_a < 1) || (kernel_dim_a % 2 == 0) || + (kernel_dim_b < 1) || (kernel_dim_b % 2 == 0)) + EIDOS_TERMINATION << "ERROR (SpatialMap::Convolve_S2): (internal error) kernel dimensions must be odd." << EidosTerminate(); + + int64_t dim_a = grid_size_[0], dim_b = grid_size_[1]; + double *new_values = (double *)malloc(dim_a * dim_b * sizeof(double)); + + if (!new_values) + EIDOS_TERMINATION << "ERROR (SpatialMap::Convolve_S2): allocation failed; you may need to raise the memory limit for SLiM." << EidosTerminate(nullptr); + + // this assumes the kernel's dimensions are symmetrical around its center, and relies on rounding (which is guaranteed) + int64_t kernel_a_offset = -(kernel_dim_a / 2), kernel_b_offset = -(kernel_dim_b / 2); + double *kernel_values = kernel.values_; + double *new_values_ptr = new_values; + + for (int64_t b = 0; b < dim_b; ++b) + { + for (int64_t a = 0; a < dim_a; ++a) + { + // calculate the kernel's effect at point (a,b) + double kernel_total = 0.0; + double conv_total = 0.0; + + for (int64_t kernel_a = 0; kernel_a < kernel_dim_a; kernel_a++) + { + int64_t conv_a = a + kernel_a + kernel_a_offset; + + // handle bounds: either clip or wrap + if ((conv_a < 0) || (conv_a >= dim_a)) + continue; + + for (int64_t kernel_b = 0; kernel_b < kernel_dim_b; kernel_b++) + { + int64_t conv_b = b + kernel_b + kernel_b_offset; + + // handle bounds: either clip or wrap + if ((conv_b < 0) || (conv_b >= dim_b)) + continue; + + // this point is within bounds; add it in to the totals + double kernel_value = kernel_values[kernel_a + kernel_b * kernel_dim_a]; + double pixel_value = values_[conv_a + conv_b * dim_a]; + + // we keep a total of the kernel values that were within bounds, for this point + kernel_total += kernel_value; + + // and we keep a total of the convolution - kernel values times pixel values + conv_total += kernel_value * pixel_value; + } + } + + *(new_values_ptr++) = ((kernel_total > 0) ? (conv_total / kernel_total) : 0); + } + } + + TakeOverMallocedValues(new_values, 2, grid_size_); // takes new_values from us +} + +void SpatialMap::Convolve_S3(SpatialKernel &kernel) +{ + if (spatiality_ != 3) + EIDOS_TERMINATION << "ERROR (SpatialMap::Convolve_S3): (internal error) map spatiality 3 required." << EidosTerminate(); + if (kernel.dimensionality_ != 3) + EIDOS_TERMINATION << "ERROR (SpatialMap::Convolve_S3): (internal error) kernel dimensionality 3 required." << EidosTerminate(); + + int64_t kernel_dim_a = kernel.dim[0]; + int64_t kernel_dim_b = kernel.dim[1]; + int64_t kernel_dim_c = kernel.dim[2]; + + if ((kernel_dim_a < 1) || (kernel_dim_a % 2 == 0) || + (kernel_dim_b < 1) || (kernel_dim_b % 2 == 0) || + (kernel_dim_c < 1) || (kernel_dim_c % 2 == 0)) + EIDOS_TERMINATION << "ERROR (SpatialMap::Convolve_S3): (internal error) kernel dimensions must be odd." << EidosTerminate(); + + int64_t dim_a = grid_size_[0], dim_b = grid_size_[1], dim_c = grid_size_[2]; + double *new_values = (double *)malloc(dim_a * dim_b * dim_c * sizeof(double)); + + if (!new_values) + EIDOS_TERMINATION << "ERROR (SpatialMap::Convolve_S3): allocation failed; you may need to raise the memory limit for SLiM." << EidosTerminate(nullptr); + + // this assumes the kernel's dimensions are symmetrical around its center, and relies on rounding (which is guaranteed) + int64_t kernel_a_offset = -(kernel_dim_a / 2), kernel_b_offset = -(kernel_dim_b / 2), kernel_c_offset = -(kernel_dim_c / 2); + double *kernel_values = kernel.values_; + double *new_values_ptr = new_values; + + for (int64_t c = 0; c < dim_c; ++c) + { + for (int64_t b = 0; b < dim_b; ++b) + { + for (int64_t a = 0; a < dim_a; ++a) + { + // calculate the kernel's effect at point (a,b,c) + double kernel_total = 0.0; + double conv_total = 0.0; + + for (int64_t kernel_a = 0; kernel_a < kernel_dim_a; kernel_a++) + { + int64_t conv_a = a + kernel_a + kernel_a_offset; + + // handle bounds: either clip or wrap + if ((conv_a < 0) || (conv_a >= dim_a)) + continue; + + for (int64_t kernel_b = 0; kernel_b < kernel_dim_b; kernel_b++) + { + int64_t conv_b = b + kernel_b + kernel_b_offset; + + // handle bounds: either clip or wrap + if ((conv_b < 0) || (conv_b >= dim_b)) + continue; + + for (int64_t kernel_c = 0; kernel_c < kernel_dim_c; kernel_c++) + { + int64_t conv_c = c + kernel_c + kernel_c_offset; + + // handle bounds: either clip or wrap + if ((conv_c < 0) || (conv_c >= dim_c)) + continue; + + // 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]; + double pixel_value = values_[conv_a + conv_b * dim_a + conv_c * dim_a * dim_b]; + + // we keep a total of the kernel values that were within bounds, for this point + kernel_total += kernel_value; + + // and we keep a total of the convolution - kernel values times pixel values + conv_total += kernel_value * pixel_value; + } + } + } + + *(new_values_ptr++) = ((kernel_total > 0) ? (conv_total / kernel_total) : 0); + } + } + } + + TakeOverMallocedValues(new_values, 3, grid_size_); // takes new_values from us +} + // // Eidos support @@ -765,6 +966,7 @@ EidosValue_SP SpatialMap::ExecuteInstanceMethod(EidosGlobalStringID p_method_id, case gID_mapColor: return ExecuteMethod_mapColor(p_method_id, p_arguments, p_interpreter); case gID_mapImage: return ExecuteMethod_mapImage(p_method_id, p_arguments, p_interpreter); case gID_mapValue: return ExecuteMethod_mapValue(p_method_id, p_arguments, p_interpreter); + case gID_smoothValues: return ExecuteMethod_smoothValues(p_method_id, p_arguments, p_interpreter); default: return super::ExecuteInstanceMethod(p_method_id, p_arguments, p_interpreter); } } @@ -786,7 +988,7 @@ EidosValue_SP SpatialMap::ExecuteMethod_changeValues(EidosGlobalStringID p_metho return gStaticEidosValueVOID; } -// ********************* - (void)interpolateValues(integer$ factor) +// ********************* - (void)interpolateValues(integer$ factor, [string$ method = "linear"]) // EidosValue_SP SpatialMap::ExecuteMethod_interpolateValues(EidosGlobalStringID p_method_id, const std::vector &p_arguments, EidosInterpreter &p_interpreter) { @@ -794,19 +996,31 @@ EidosValue_SP SpatialMap::ExecuteMethod_interpolateValues(EidosGlobalStringID p_ EidosValue *factor_value = p_arguments[0].get(); int64_t factor = factor_value->IntAtIndex(0, nullptr); - if ((factor < 2) || (factor > 100)) - EIDOS_TERMINATION << "ERROR (SpatialMap::ExecuteMethod_interpolateValues): interpolateValues() requires factor to be in [2, 100]." << EidosTerminate(); + // the upper limit here is arbitrary, but the goal is to prevent users from blowing up their memory usage unintentionally + if ((factor < 2) || (factor > 201)) + EIDOS_TERMINATION << "ERROR (SpatialMap::ExecuteMethod_interpolateValues): interpolateValues() requires factor to be in [2, 201]." << EidosTerminate(); + + EidosValue_String *method_value = (EidosValue_String *)p_arguments[1].get(); + const std::string &method_string = method_value->StringRefAtIndex(0, nullptr); + int method; + + if (method_string == "nearest") + method = 0; + else if (method_string == "linear") + method = 1; + else + EIDOS_TERMINATION << "ERROR (SpatialMap::ExecuteMethod_interpolateValues): interpolateValues() requires method to be 'nearest' or 'linear'." << EidosTerminate(); // Temporarily force interpolation on bool old_interpolate = interpolate_; - interpolate_ = true; + interpolate_ = (method ? true : false); // Generate the new values and set them switch (spatiality_) { case 1: { - int64_t dim_a = factor * grid_size_[0] - 1; + int64_t dim_a = (factor * (grid_size_[0] - 1)) + 1; double *new_values = (double *)malloc(dim_a * sizeof(double)); double *new_values_ptr = new_values; double point_vec[1]; @@ -826,7 +1040,7 @@ EidosValue_SP SpatialMap::ExecuteMethod_interpolateValues(EidosGlobalStringID p_ } case 2: { - int64_t dim_a = factor * grid_size_[0] - 1, dim_b = factor * grid_size_[1] - 1; + 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 *new_values_ptr = new_values; double point_vec[2]; @@ -852,7 +1066,7 @@ EidosValue_SP SpatialMap::ExecuteMethod_interpolateValues(EidosGlobalStringID p_ } case 3: { - int64_t dim_a = factor * grid_size_[0] - 1, dim_b = factor * grid_size_[1] - 1, dim_c = factor * grid_size_[2] - 1; + int64_t dim_a = (factor * (grid_size_[0] - 1)) + 1, dim_b = (factor * (grid_size_[1] - 1)) + 1, dim_c = (factor * (grid_size_[2] - 1)) + 1; double *new_values = (double *)malloc(dim_a * dim_b * dim_c * sizeof(double)); double *new_values_ptr = new_values; double point_vec[3]; @@ -1124,6 +1338,45 @@ EidosValue_SP SpatialMap::ExecuteMethod_mapValue(EidosGlobalStringID p_method_id return EidosValue_SP(float_singleton_result); } +// ********************* - (void)smoothValues(float$ maxDistance, string$ functionType, ...) +// +EidosValue_SP SpatialMap::ExecuteMethod_smoothValues(EidosGlobalStringID p_method_id, const std::vector &p_arguments, EidosInterpreter &p_interpreter) +{ +#pragma unused (p_method_id, p_arguments, p_interpreter) + // Our arguments go to Kernel::Kernel(), which creates the kernel object that we use + EidosValue *maxDistance_value = p_arguments[0].get(); + double max_distance = maxDistance_value->FloatAtIndex(0, nullptr); + + SpatialKernel kernel(spatiality_, max_distance, p_arguments, 1); // uses our arguments starting at index 1 + + // Ask the kernel to create a discrete grid of values, at our spatial scale (we define the + // relationship between spatial bounds and pixels, used by the kernel to make its grid) + kernel.CalculateGridValues(*this); + + //std::cout << kernel << std::endl; + + // Generate the new spatial map values and set them into ourselves + switch (spatiality_) + { + case 1: + Convolve_S1(kernel); break; + case 2: + Convolve_S2(kernel); break; + case 3: + Convolve_S3(kernel); break; + + default: break; + } + + // Reassess our min and max if we're using the default grayscale color map; + // otherwise they are user-supplied and should not be modified + if (n_colors_ == 0) + SetAutomaticColorMinMax(); + + return gStaticEidosValueVOID; +} + + // // Object instantiation @@ -1196,10 +1449,11 @@ const std::vector *SpatialMap_Class::Methods(void) con methods = new std::vector(*super::Methods()); methods->emplace_back((EidosInstanceMethodSignature *)(new EidosInstanceMethodSignature(gStr_changeValues, kEidosValueMaskVOID))->AddNumeric("values")); - methods->emplace_back((EidosInstanceMethodSignature *)(new EidosInstanceMethodSignature(gStr_interpolateValues, kEidosValueMaskVOID))->AddInt_S("factor")); + methods->emplace_back((EidosInstanceMethodSignature *)(new EidosInstanceMethodSignature(gStr_interpolateValues, kEidosValueMaskVOID))->AddInt_S("factor")->AddString_OS("method", EidosValue_String_SP(new (gEidosValuePool->AllocateChunk()) EidosValue_String_singleton("linear")))); methods->emplace_back((EidosInstanceMethodSignature *)(new EidosInstanceMethodSignature(gStr_mapColor, kEidosValueMaskString))->AddNumeric("value")); methods->emplace_back((EidosInstanceMethodSignature *)(new EidosInstanceMethodSignature(gStr_mapImage, kEidosValueMaskObject | kEidosValueMaskSingleton, gEidosImage_Class))->AddInt_OSN(gEidosStr_width, gStaticEidosValueNULL)->AddInt_OSN(gEidosStr_height, gStaticEidosValueNULL)->AddLogical_OS("centers", gStaticEidosValue_LogicalF)->AddLogical_OS(gEidosStr_color, gStaticEidosValue_LogicalT)); methods->emplace_back((EidosInstanceMethodSignature *)(new EidosInstanceMethodSignature(gStr_mapValue, kEidosValueMaskFloat))->AddFloat("point")); + methods->emplace_back((EidosInstanceMethodSignature *)(new EidosInstanceMethodSignature(gStr_smoothValues, kEidosValueMaskVOID))->AddFloat_S("maxDistance")->AddString_S("functionType")->AddEllipsis()); std::sort(methods->begin(), methods->end(), CompareEidosCallSignatures); } diff --git a/core/spatial_map.h b/core/spatial_map.h index 606971bf8..ad239feb0 100644 --- a/core/spatial_map.h +++ b/core/spatial_map.h @@ -36,10 +36,11 @@ #include "eidos_class_Dictionary.h" class Subpopulation; +class SpatialKernel; #pragma mark - -#pragma mark _SpatialMap +#pragma mark SpatialMap #pragma mark - extern EidosClass *gSLiM_SpatialMap_Class; @@ -103,6 +104,10 @@ class SpatialMap : public EidosDictionaryRetained void ColorForValue(double p_value, double *p_rgb_ptr); void ColorForValue(double p_value, float *p_rgb_ptr); + void Convolve_S1(SpatialKernel &kernel); + void Convolve_S2(SpatialKernel &kernel); + void Convolve_S3(SpatialKernel &kernel); + // // Eidos support // @@ -118,6 +123,7 @@ class SpatialMap : public EidosDictionaryRetained EidosValue_SP ExecuteMethod_mapColor(EidosGlobalStringID p_method_id, const std::vector &p_arguments, EidosInterpreter &p_interpreter); EidosValue_SP ExecuteMethod_mapImage(EidosGlobalStringID p_method_id, const std::vector &p_arguments, EidosInterpreter &p_interpreter); EidosValue_SP ExecuteMethod_mapValue(EidosGlobalStringID p_method_id, const std::vector &p_arguments, EidosInterpreter &p_interpreter); + EidosValue_SP ExecuteMethod_smoothValues(EidosGlobalStringID p_method_id, const std::vector &p_arguments, EidosInterpreter &p_interpreter); }; class SpatialMap_Class : public EidosDictionaryRetained_Class