From d71d392f4499a13f2df3f07a0c6a2403ff1af5eb Mon Sep 17 00:00:00 2001 From: Ben Haller Date: Wed, 13 Sep 2023 19:37:05 -0400 Subject: [PATCH] add t-distribution to SpatialKernel, fix up kernel code --- QtSLiM/QtSLiMTablesDrawer.cpp | 4 + QtSLiM/help/SLiMHelpClasses.html | 10 +- SLiM.xcodeproj/project.pbxproj | 36 ++++ SLiMgui/SLiMHelpClasses.rtf | 62 +++++-- SLiMgui/SLiMWindowController.mm | 4 + VERSIONS | 1 + core/interaction_type.cpp | 43 ++++- core/interaction_type.h | 3 +- core/spatial_kernel.cpp | 273 +++++++++++++++++++++++-------- core/spatial_kernel.h | 17 +- eidos/eidos_globals.cpp | 1 + eidos/eidos_globals.h | 2 + gsl/gsl.pro | 2 + gsl/randist/chisq.c | 61 +++++++ gsl/randist/tdist.c | 80 +++++++++ 15 files changed, 498 insertions(+), 101 deletions(-) create mode 100644 gsl/randist/chisq.c create mode 100644 gsl/randist/tdist.c diff --git a/QtSLiM/QtSLiMTablesDrawer.cpp b/QtSLiM/QtSLiMTablesDrawer.cpp index 654403be1..6a84d6089 100644 --- a/QtSLiM/QtSLiMTablesDrawer.cpp +++ b/QtSLiM/QtSLiMTablesDrawer.cpp @@ -911,6 +911,7 @@ QVariant QtSLiMInteractionTypeTableModel::data(const QModelIndex &p_index, int r case SpatialKernelType::kExponential: return QVariant(QString("exp")); case SpatialKernelType::kNormal: return QVariant(QString("normal")); case SpatialKernelType::kCauchy: return QVariant(QString("Cauchy")); + case SpatialKernelType::kStudentsT: return QVariant(QString("Student's t")); } } else if (p_index.column() == 3) @@ -935,6 +936,9 @@ QVariant QtSLiMInteractionTypeTableModel::data(const QModelIndex &p_index, int r case SpatialKernelType::kCauchy: paramString += QString(", γ=%1").arg(interactionType->if_param2_, 0, 'f', 3); break; + case SpatialKernelType::kStudentsT: + paramString += QString(", ν=%1, σ=%2").arg(interactionType->if_param2_, 0, 'f', 3).arg(interactionType->if_param3_, 0, 'f', 3); + break; } return QVariant(paramString); diff --git a/QtSLiM/help/SLiMHelpClasses.html b/QtSLiM/help/SLiMHelpClasses.html index f8d9d10fc..08ac0d50d 100644 --- a/QtSLiM/help/SLiMHelpClasses.html +++ b/QtSLiM/help/SLiMHelpClasses.html @@ -500,7 +500,7 @@

Returns the number of individuals in exerterSubpop that are within the maximum interaction distance according to the distance metric of the InteractionType.  The subpopulation may be supplied either as an integer ID, or as a Subpopulation object.  The evaluate() method must have been previously called for exerterSubpop, and positions saved at evaluation time will be used.  If the InteractionType is non-spatial, this method may not be called.

This method is similar to nearestNeighborsOfPoint() (when passed a large count so as to guarantee that all neighbors are returned), but this method returns only a count of the individuals, not a vector containing the individuals.

– (void)setInteractionFunction(string$ functionType, ...)

-

Set the function used to translate spatial distances into interaction strengths for an interaction type.  The functionType may be "f", in which case the ellipsis ... should supply a numeric$ fixed interaction strength; "l", in which case the ellipsis should supply a numeric$ maximum strength for a linear function; "e", in which case the ellipsis should supply a numeric$ maximum strength and a numeric$ lambda (shape) parameter for a negative exponential function; ; "n", in which case the ellipsis should supply a numeric$ maximum strength and a numeric$ sigma (standard deviation) parameter for a Gaussian function; or "c", in which case the ellipsis should supply a numeric$ maximum strength and a numeric$ scale parameter for a Cauchy distribution function.  Non-spatial interactions must use function type "f", since no distance values are available in that case.

+

Set the function used to translate spatial distances into interaction strengths for an interaction type.  The functionType may be "f", in which case the ellipsis ... should supply a numeric$ fixed interaction strength; "l", in which case the ellipsis should supply a numeric$ maximum strength for a linear function; "e", in which case the ellipsis should supply a numeric$ maximum strength and a numeric$ lambda (shape) parameter for a negative exponential function; "n", in which case the ellipsis should supply a numeric$ maximum strength and a numeric$ sigma (standard deviation) parameter for a Gaussian function; "c", in which case the ellipsis should supply a numeric$ maximum strength and a numeric$ scale parameter for a Cauchy distribution function; or "t", in which case the ellipsis should supply a numeric$ maximum strength, a numeric$ degrees of freedom, and a numeric$ scale parameter for a t-distribution function.  Non-spatial interactions must use function type "f", since no distance values are available in that case.

The interaction function for an interaction type is normally a constant in simulations; in any case, it cannot be changed when an interaction has already been evaluated, so either it should be set prior to evaluation, or unevaluate() should be called first.

– (float)strength(object<Individual>$ receiver, [No<Individual> exerters = NULL])

Returns a vector containing the interaction strengths exerted upon receiver by the individuals in exerters.  If exerters is NULL (the default), then a vector of the interaction strengths exerted by all individuals in the subpopulation of receiver (including receiver itself) is returned; this case may be handled much more efficiently than if a vector of all individuals in the subpopulation is explicitly provided.  Otherwise, all individuals in exerters must belong to a single subpopulation (but not necessarily the same subpopulation as receiver).  The evaluate() method must have been previously called for the receiver and exerter subpopulations, and positions saved at evaluation time will be used.

@@ -735,11 +735,11 @@

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.

-

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 a flat kernel (type "f") or a Gaussian kernel (type "n") are supported.

-

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.

+

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", or "c", 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; or "c", in which case the ellipsis should supply a numeric$ scale parameter for a Cauchy distribution function.  See the InteractionType class documentation for discussions of these kernel types.

+

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 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)

diff --git a/SLiM.xcodeproj/project.pbxproj b/SLiM.xcodeproj/project.pbxproj index ec28f1dc9..8e2a6740b 100644 --- a/SLiM.xcodeproj/project.pbxproj +++ b/SLiM.xcodeproj/project.pbxproj @@ -1195,6 +1195,22 @@ 98D524661F2E6B0B005AD9A6 /* prettyprint.pdf in Resources */ = {isa = PBXBuildFile; fileRef = 98D524621F2E6AFB005AD9A6 /* prettyprint.pdf */; }; 98D524691F2EB4DD005AD9A6 /* EidosPrettyprinter.mm in Sources */ = {isa = PBXBuildFile; fileRef = 98D524681F2EB4DD005AD9A6 /* EidosPrettyprinter.mm */; }; 98D5246A1F2EB4DD005AD9A6 /* EidosPrettyprinter.mm in Sources */ = {isa = PBXBuildFile; fileRef = 98D524681F2EB4DD005AD9A6 /* EidosPrettyprinter.mm */; }; + 98D7D65C2AB24C40002AFE34 /* tdist.c in Sources */ = {isa = PBXBuildFile; fileRef = 98D7D65B2AB24C40002AFE34 /* tdist.c */; }; + 98D7D65D2AB24C40002AFE34 /* tdist.c in Sources */ = {isa = PBXBuildFile; fileRef = 98D7D65B2AB24C40002AFE34 /* tdist.c */; }; + 98D7D65E2AB24C40002AFE34 /* tdist.c in Sources */ = {isa = PBXBuildFile; fileRef = 98D7D65B2AB24C40002AFE34 /* tdist.c */; }; + 98D7D65F2AB24C40002AFE34 /* tdist.c in Sources */ = {isa = PBXBuildFile; fileRef = 98D7D65B2AB24C40002AFE34 /* tdist.c */; }; + 98D7D6602AB24C40002AFE34 /* tdist.c in Sources */ = {isa = PBXBuildFile; fileRef = 98D7D65B2AB24C40002AFE34 /* tdist.c */; }; + 98D7D6612AB24C40002AFE34 /* tdist.c in Sources */ = {isa = PBXBuildFile; fileRef = 98D7D65B2AB24C40002AFE34 /* tdist.c */; }; + 98D7D6622AB24C40002AFE34 /* tdist.c in Sources */ = {isa = PBXBuildFile; fileRef = 98D7D65B2AB24C40002AFE34 /* tdist.c */; }; + 98D7D6632AB24C40002AFE34 /* tdist.c in Sources */ = {isa = PBXBuildFile; fileRef = 98D7D65B2AB24C40002AFE34 /* tdist.c */; }; + 98D7D6652AB24CBC002AFE34 /* chisq.c in Sources */ = {isa = PBXBuildFile; fileRef = 98D7D6642AB24CBC002AFE34 /* chisq.c */; }; + 98D7D6662AB24CBC002AFE34 /* chisq.c in Sources */ = {isa = PBXBuildFile; fileRef = 98D7D6642AB24CBC002AFE34 /* chisq.c */; }; + 98D7D6672AB24CBC002AFE34 /* chisq.c in Sources */ = {isa = PBXBuildFile; fileRef = 98D7D6642AB24CBC002AFE34 /* chisq.c */; }; + 98D7D6682AB24CBC002AFE34 /* chisq.c in Sources */ = {isa = PBXBuildFile; fileRef = 98D7D6642AB24CBC002AFE34 /* chisq.c */; }; + 98D7D6692AB24CBC002AFE34 /* chisq.c in Sources */ = {isa = PBXBuildFile; fileRef = 98D7D6642AB24CBC002AFE34 /* chisq.c */; }; + 98D7D66A2AB24CBC002AFE34 /* chisq.c in Sources */ = {isa = PBXBuildFile; fileRef = 98D7D6642AB24CBC002AFE34 /* chisq.c */; }; + 98D7D66B2AB24CBC002AFE34 /* chisq.c in Sources */ = {isa = PBXBuildFile; fileRef = 98D7D6642AB24CBC002AFE34 /* chisq.c */; }; + 98D7D66C2AB24CBC002AFE34 /* chisq.c in Sources */ = {isa = PBXBuildFile; fileRef = 98D7D6642AB24CBC002AFE34 /* chisq.c */; }; 98D7EB8528CE557C00DEAAC4 /* coerce.c in Sources */ = {isa = PBXBuildFile; fileRef = 98C821201C7A980000548839 /* coerce.c */; }; 98D7EB8628CE557C00DEAAC4 /* compress.c in Sources */ = {isa = PBXBuildFile; fileRef = 98076603244934A800F6CBB4 /* compress.c */; }; 98D7EB8728CE557C00DEAAC4 /* eidos_value.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 985724A31AD435310047C223 /* eidos_value.cpp */; }; @@ -1969,6 +1985,8 @@ 98D524621F2E6AFB005AD9A6 /* prettyprint.pdf */ = {isa = PBXFileReference; lastKnownFileType = image.pdf; path = prettyprint.pdf; sourceTree = ""; }; 98D524671F2EB4DD005AD9A6 /* EidosPrettyprinter.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = EidosPrettyprinter.h; sourceTree = ""; }; 98D524681F2EB4DD005AD9A6 /* EidosPrettyprinter.mm */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.objcpp; path = EidosPrettyprinter.mm; sourceTree = ""; }; + 98D7D65B2AB24C40002AFE34 /* tdist.c */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.c; path = tdist.c; sourceTree = ""; }; + 98D7D6642AB24CBC002AFE34 /* chisq.c */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.c; path = chisq.c; sourceTree = ""; }; 98D7EBEE28CE557C00DEAAC4 /* eidos_multi */ = {isa = PBXFileReference; explicitFileType = "compiled.mach-o.executable"; includeInIndex = 0; path = eidos_multi; sourceTree = BUILT_PRODUCTS_DIR; }; 98D7ED2D28CE58FC00DEAAC4 /* slim_multi */ = {isa = PBXFileReference; explicitFileType = "compiled.mach-o.executable"; includeInIndex = 0; path = slim_multi; sourceTree = BUILT_PRODUCTS_DIR; }; 98DB3D6D1E6122AE00E2C200 /* interaction_type.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = interaction_type.cpp; sourceTree = ""; }; @@ -2718,6 +2736,7 @@ 98C820F71C7A980000548839 /* beta.c */, 98C820F81C7A980000548839 /* binomial_tpe.c */, 988880EB20744EE800E10172 /* cauchy.c */, + 98D7D6642AB24CBC002AFE34 /* chisq.c */, 98C821A71C7A9B1600548839 /* discrete.c */, 98C820F91C7A980000548839 /* exponential.c */, 980566E125A7C5B9008D3C7F /* fdist.c */, @@ -2731,6 +2750,7 @@ 982B50C52704048E006E91BC /* nbinomial.c */, 98C821001C7A980000548839 /* poisson.c */, 98C821B11C7A9B9F00548839 /* shuffle.c */, + 98D7D65B2AB24C40002AFE34 /* tdist.c */, 98C821011C7A980000548839 /* weibull.c */, ); path = randist; @@ -3453,6 +3473,7 @@ 98C821B51C7A9B9F00548839 /* shuffle.c in Sources */, 98CEFCF92AAFABAA00D2C9B4 /* bilinear.c in Sources */, 98C821861C7A983800548839 /* exponential.c in Sources */, + 98D7D6612AB24C40002AFE34 /* tdist.c in Sources */, 9876E60B1ED55B5000FF9762 /* erfc.c in Sources */, 981DC35F28E26F8B000ABE91 /* eidos_functions_matrices.cpp in Sources */, 9876E6101ED55C0C00FF9762 /* expint.c in Sources */, @@ -3475,6 +3496,7 @@ 982556B31BA8EF7C0054CB3F /* eidos_ast_node.cpp in Sources */, 9893C7DE1CDA2D650029AC94 /* eidos_beep.cpp in Sources */, 9807663A24493E0B00F6CBB4 /* adler32.c in Sources */, + 98D7D66A2AB24CBC002AFE34 /* chisq.c in Sources */, 9876E5FF1ED559A700FF9762 /* gauss.c in Sources */, 98C821991C7A983800548839 /* zeta.c in Sources */, 98C8218B1C7A983800548839 /* multinomial.c in Sources */, @@ -3647,6 +3669,7 @@ 98076622244934A800F6CBB4 /* trees.c in Sources */, 981DC36C28E26F8B000ABE91 /* eidos_functions_other.cpp in Sources */, 98332AC21FDBA53F00274FF0 /* xerbla.c in Sources */, + 98D7D6652AB24CBC002AFE34 /* chisq.c in Sources */, 985724D51AD481070047C223 /* eidos_test.cpp in Sources */, 98CEFD752AAFB12F00D2C9B4 /* cspline.c in Sources */, 986070EA2AACECD600FD6143 /* spatial_kernel.cpp in Sources */, @@ -3681,6 +3704,7 @@ 989790DA1AF3D0E100C6B14C /* eidos_class_TestElement.cpp in Sources */, 982663541A3BABD300A0CBBF /* main.cpp in Sources */, 981DC36828E26F8B000ABE91 /* eidos_functions_strings.cpp in Sources */, + 98D7D65C2AB24C40002AFE34 /* tdist.c in Sources */, 98DDAED6221765480038C133 /* slim_functions.cpp in Sources */, 9876E5FC1ED5599F00FF9762 /* gauss.c in Sources */, 988880EC20744EE900E10172 /* cauchy.c in Sources */, @@ -3807,6 +3831,7 @@ 981DC36228E26F8B000ABE91 /* eidos_functions_values.cpp in Sources */, 98CEFCF72AAFABAA00D2C9B4 /* bilinear.c in Sources */, 98090FA61B1B978900791DBF /* EidosValueWrapper.mm in Sources */, + 98D7D65F2AB24C40002AFE34 /* tdist.c in Sources */, 98CEFD6F2AAFABAA00D2C9B4 /* inline.c in Sources */, 98C8217B1C7A983800548839 /* zeta.c in Sources */, 9807662E24493E0B00F6CBB4 /* deflate.c in Sources */, @@ -3860,6 +3885,7 @@ 989790DC1AF3D0E100C6B14C /* eidos_class_TestElement.cpp in Sources */, 9876E60A1ED55B4F00FF9762 /* erfc.c in Sources */, 985724D11AD479010047C223 /* eidos_interpreter.cpp in Sources */, + 98D7D6682AB24CBC002AFE34 /* chisq.c in Sources */, 98C8216A1C7A983800548839 /* gauss.c in Sources */, 98C821681C7A983800548839 /* exponential.c in Sources */, 984824F3210B9F2D002402A5 /* dtrsv.c in Sources */, @@ -4041,6 +4067,7 @@ 98CF5281294A3FC900557BBA /* rowcol.c in Sources */, 98CF5282294A3FC900557BBA /* gausszig.c in Sources */, 98CF5283294A3FC900557BBA /* eidos_type_interpreter.cpp in Sources */, + 98D7D6672AB24CBC002AFE34 /* chisq.c in Sources */, 98CF5284294A3FC900557BBA /* kastore.c in Sources */, 98CF5285294A3FC900557BBA /* gzwrite.c in Sources */, 98CF5286294A3FC900557BBA /* chromosome.cpp in Sources */, @@ -4060,6 +4087,7 @@ 98CF5293294A3FC900557BBA /* slim_gui.mm in Sources */, 98CF5294294A3FC900557BBA /* mutation.cpp in Sources */, 98CF5295294A3FC900557BBA /* slim_globals.cpp in Sources */, + 98D7D65E2AB24C40002AFE34 /* tdist.c in Sources */, 98CF5296294A3FC900557BBA /* eidos_test_operators_arithmetic.cpp in Sources */, 98CF5297294A3FC900557BBA /* eidos_ast_node.cpp in Sources */, 98CF5298294A3FC900557BBA /* submatrix.c in Sources */, @@ -4142,6 +4170,7 @@ 98CF53D0294A714200557BBA /* eidos_functions_values.cpp in Sources */, 98CEFCF82AAFABAA00D2C9B4 /* bilinear.c in Sources */, 98CF53D1294A714200557BBA /* EidosValueWrapper.mm in Sources */, + 98D7D6602AB24C40002AFE34 /* tdist.c in Sources */, 98CEFD702AAFABAA00D2C9B4 /* inline.c in Sources */, 98CF53D2294A714200557BBA /* zeta.c in Sources */, 98CF53D3294A714200557BBA /* deflate.c in Sources */, @@ -4195,6 +4224,7 @@ 98CF53FB294A714200557BBA /* eidos_class_TestElement.cpp in Sources */, 98CF53FC294A714200557BBA /* erfc.c in Sources */, 98CF53FD294A714200557BBA /* eidos_interpreter.cpp in Sources */, + 98D7D6692AB24CBC002AFE34 /* chisq.c in Sources */, 98CF53FE294A714200557BBA /* gauss.c in Sources */, 98CF53FF294A714200557BBA /* exponential.c in Sources */, 98CF5400294A714200557BBA /* dtrsv.c in Sources */, @@ -4376,6 +4406,7 @@ 98332AED1FDBC29400274FF0 /* rowcol.c in Sources */, 98C8214D1C7A983700548839 /* gausszig.c in Sources */, 9893C7EA1CDFE9870029AC94 /* eidos_type_interpreter.cpp in Sources */, + 98D7D6662AB24CBC002AFE34 /* chisq.c in Sources */, 987D19A6209A53850030D28D /* kastore.c in Sources */, 98076621244934A800F6CBB4 /* gzwrite.c in Sources */, 98D4C1DD1A6F543600FFB083 /* chromosome.cpp in Sources */, @@ -4395,6 +4426,7 @@ 98EBCD1421F3CFC600B385CF /* slim_gui.mm in Sources */, 98D4C1D81A6F542600FFB083 /* mutation.cpp in Sources */, 98D4C1D31A6F541200FFB083 /* slim_globals.cpp in Sources */, + 98D7D65D2AB24C40002AFE34 /* tdist.c in Sources */, 985F3F0224BA307200E712E0 /* eidos_test_operators_arithmetic.cpp in Sources */, 98A513551B66B6CA005A753D /* eidos_ast_node.cpp in Sources */, 98332AF11FDBC36300274FF0 /* submatrix.c in Sources */, @@ -4435,6 +4467,7 @@ 98D7EB9528CE557C00DEAAC4 /* eidos_test.cpp in Sources */, 98CEFCFA2AAFABAA00D2C9B4 /* bilinear.c in Sources */, 98D7EB9628CE557C00DEAAC4 /* oper.c in Sources */, + 98D7D6622AB24C40002AFE34 /* tdist.c in Sources */, 98D7EB9728CE557C00DEAAC4 /* eidos_class_Image.cpp in Sources */, 98D7EB9828CE557C00DEAAC4 /* swap.c in Sources */, 98D7EB9928CE557C00DEAAC4 /* eidos_class_Dictionary.cpp in Sources */, @@ -4457,6 +4490,7 @@ 98D7EBA628CE557C00DEAAC4 /* binomial_tpe.c in Sources */, 98D7EBA728CE557C00DEAAC4 /* cauchy.c in Sources */, 98D7EBA828CE557C00DEAAC4 /* infnan.c in Sources */, + 98D7D66B2AB24CBC002AFE34 /* chisq.c in Sources */, 98D7EBA928CE557C00DEAAC4 /* eidos_functions.cpp in Sources */, 98D7EBAA28CE557C00DEAAC4 /* eidos_interpreter.cpp in Sources */, 98D7EBAB28CE557C00DEAAC4 /* mt.c in Sources */, @@ -4629,6 +4663,7 @@ 98D7ECE128CE58FC00DEAAC4 /* geometric.c in Sources */, 98D7ECE228CE58FC00DEAAC4 /* eidos_token.cpp in Sources */, 98D7ECE328CE58FC00DEAAC4 /* core.c in Sources */, + 98D7D66C2AB24CBC002AFE34 /* chisq.c in Sources */, 98D7ECE428CE58FC00DEAAC4 /* minmax.c in Sources */, 98CEFD7C2AAFB12F00D2C9B4 /* cspline.c in Sources */, 986070ED2AACECD600FD6143 /* spatial_kernel.cpp in Sources */, @@ -4663,6 +4698,7 @@ 98D7ECFE28CE58FC00DEAAC4 /* tdist.c in Sources */, 98D7ECFF28CE58FC00DEAAC4 /* polymorphism.cpp in Sources */, 98D7ED0028CE58FC00DEAAC4 /* eidos_class_TestElement.cpp in Sources */, + 98D7D6632AB24C40002AFE34 /* tdist.c in Sources */, 98D7ED0128CE58FC00DEAAC4 /* main.cpp in Sources */, 98D7ED0228CE58FC00DEAAC4 /* slim_functions.cpp in Sources */, 98D7ED0328CE58FC00DEAAC4 /* gauss.c in Sources */, diff --git a/SLiMgui/SLiMHelpClasses.rtf b/SLiMgui/SLiMHelpClasses.rtf index 71d99de02..b7fa93c73 100644 --- a/SLiMgui/SLiMHelpClasses.rtf +++ b/SLiMgui/SLiMHelpClasses.rtf @@ -4140,7 +4140,7 @@ This method is similar to \f5 \ \pard\pardeftab543\li547\ri720\sb60\sa60\partightenfactor0 -\f4\fs20 \cf0 Set the function used to translate spatial distances into interaction strengths for an interaction type. The +\f4\fs20 \cf2 Set the function used to translate spatial distances into interaction strengths for an interaction type. The \f3\fs18 functionType \f4\fs20 may be \f3\fs18 "f" @@ -4158,24 +4158,32 @@ This method is similar to \f3\fs18 numeric$ \f4\fs20 maximum strength and a \f3\fs18 numeric$ -\f4\fs20 lambda (shape) parameter for a negative exponential function; \cf2 \expnd0\expndtw0\kerning0 -; +\f4\fs20 lambda (shape) parameter for a negative exponential function; \f3\fs18 "n" \f4\fs20 , in which case the ellipsis should supply a \f3\fs18 numeric$ \f4\fs20 maximum strength and a \f3\fs18 numeric$ -\f4\fs20 sigma (standard deviation) parameter for a Gaussian function; or +\f4\fs20 sigma (standard deviation) parameter for a Gaussian function; \f3\fs18 "c" \f4\fs20 , in which case the ellipsis should supply a \f3\fs18 numeric$ \f4\fs20 maximum strength and a \f3\fs18 numeric$ -\f4\fs20 scale parameter for a Cauchy distribution function\cf0 \kerning1\expnd0\expndtw0 . Non-spatial interactions must use function type +\f4\fs20 scale parameter for a Cauchy distribution function; or +\f3\fs18 "t" +\f4\fs20 , in which case the ellipsis should supply a +\f3\fs18 numeric$ +\f4\fs20 maximum strength, a +\f3\fs18 numeric$ +\f4\fs20 degrees of freedom, and a +\f3\fs18 numeric$ +\f4\fs20 scale parameter for a +\f1\i t +\f4\i0 -distribution function. Non-spatial interactions must use function type \f3\fs18 "f" \f4\fs20 , since no distance values are available in that case.\ -\pard\pardeftab543\li547\ri720\sb60\sa60\partightenfactor0 -\cf2 The interaction function for an interaction type is normally a constant in simulations; in any case, it cannot be changed when an interaction has already been evaluated, so either it should be set prior to evaluation, or +The interaction function for an interaction type is normally a constant in simulations; in any case, it cannot be changed when an interaction has already been evaluated, so either it should be set prior to evaluation, or \f3\fs18 unevaluate() \f4\fs20 should be called first. \f5 \cf0 \ @@ -6080,12 +6088,22 @@ The kernel is specified with a kernel shape, \f3\fs18 functionType \f4\fs20 , followed by zero or more ellipsis arguments; see \f3\fs18 smooth() -\f4\fs20 for further information. For this method, at present only a flat kernel (type +\f4\fs20 for further information. For this method, at present only kernel types \f3\fs18 "f" -\f4\fs20 ) or a Gaussian kernel (type +\f4\fs20 , +\f3\fs18 "l" +\f4\fs20 , +\f3\fs18 "e" +\f4\fs20 , \f3\fs18 "n" -\f4\fs20 ) are supported.\ -This method can be used to find points in the vicinity of individuals that are favorable \'96 possessing more resources, or better environmental conditions, etc. It can also be used to guide the dispersal or foraging behavior of individuals.\ +\f4\fs20 , and +\f3\fs18 "t" +\f4\fs20 are supported, and type +\f3\fs18 "t" +\f4\fs20 is not presently supported for 3D kernels.\ +This method can be used to find points in the vicinity of individuals that are favorable \'96 possessing more resources, or better environmental conditions, etc. It can also be used to guide the dispersal or foraging behavior of individuals. See +\f3\fs18 sampleImprovedNearbyPoint() +\f4\fs20 for a variant that may be useful for directed movement across a landscape.\ \pard\pardeftab397\li720\fi-446\ri720\sb180\sa60\partightenfactor0 \f3\fs18 \cf2 \'96\'a0(object$)smooth(float$\'a0maxDistance, string$\'a0functionType, ...)\ @@ -6103,8 +6121,10 @@ This method can be used to find points in the vicinity of individuals that are f \f3\fs18 "e" \f4\fs20 , \f3\fs18 "n" -\f4\fs20 , or +\f4\fs20 , \f3\fs18 "c" +\f4\fs20 , or +\f3\fs18 "t" \f4\fs20 , and additional parameters in the ellipsis \f3\fs18 ... \f4\fs20 that depend upon the kernel type and further specify its shape. The target spatial map is returned, to allow easy chaining of operations.\ @@ -6126,11 +6146,19 @@ The kernel specification is similar to that for the \f3\fs18 "n" \f4\fs20 , in which case the ellipsis should supply a \f3\fs18 numeric$ -\f4\fs20 sigma (standard deviation) parameter for a Gaussian function; or +\f4\fs20 sigma (standard deviation) parameter for a Gaussian function; \f3\fs18 "c" \f4\fs20 , in which case the ellipsis should supply a \f3\fs18 numeric$ -\f4\fs20 scale parameter for a Cauchy distribution function. See the +\f4\fs20 scale parameter for a Cauchy distribution function; or +\f3\fs18 "t" +\f4\fs20 , in which case the ellipsis should supply a +\f3\fs18 numeric$ +\f4\fs20 degrees of freedom and a +\f3\fs18 numeric$ +\f4\fs20 scale parameter for a +\f1\i t +\f4\i0 -distribution function. See the \f3\fs18 InteractionType \f4\fs20 class documentation for discussions of these kernel types.\ Distance metrics specified to this method, such as @@ -6789,7 +6817,8 @@ Beginning with SLiM 3.3, the \f4\i0 include the nucleotides associated with any nucleotide-based mutations; the \f3\fs18 ancestralNucleotides \f4\fs20 flag governs only the ancestral sequence.\ -\kerning1\expnd0\expndtw0 Beginning with SLiM 3.5, the +\pard\pardeftab543\li547\ri720\sb60\sa60\partightenfactor0 +\cf2 \kerning1\expnd0\expndtw0 Beginning with SLiM 3.5, the \f3\fs18 pedigreeIDs \f4\fs20 parameter may be used to request that pedigree IDs be written out (and read in by \f3\fs18 readFromPopulationFile() @@ -6919,7 +6948,8 @@ As of SLiM 3.0, this method will read and restore the ages of individuals if tha As of SLiM 3.3, this method will restore the nucleotides of nucleotide-based mutations, and will restore the ancestral nucleotide sequence, if that information is present in the output file. Loading an output file that contains nucleotide information in a non-nucleotide-based model, and \f1\i vice versa \f4\i0 , will produce an error.\ -\kerning1\expnd0\expndtw0 As of SLiM 3.5, this method will read and restore the pedigree IDs of individuals and genomes if that information is present in the output file (as requested with +\pard\pardeftab543\li547\ri720\sb60\sa60\partightenfactor0 +\cf2 \kerning1\expnd0\expndtw0 As of SLiM 3.5, this method will read and restore the pedigree IDs of individuals and genomes if that information is present in the output file (as requested with \f3\fs18 outputFull(pedigreeIDs=T) \f4\fs20 ) \f1\i and diff --git a/SLiMgui/SLiMWindowController.mm b/SLiMgui/SLiMWindowController.mm index 30dd47550..1cd6c744a 100644 --- a/SLiMgui/SLiMWindowController.mm +++ b/SLiMgui/SLiMWindowController.mm @@ -4946,6 +4946,7 @@ - (id)tableView:(NSTableView *)aTableView objectValueForTableColumn:(NSTableColu case SpatialKernelType::kExponential: return @"exp"; case SpatialKernelType::kNormal: return @"normal"; case SpatialKernelType::kCauchy: return @"Cauchy"; + case SpatialKernelType::kStudentsT: return @"Student's t"; } } else if (aTableColumn == interactionTypeIFParamsColumn) @@ -4970,6 +4971,9 @@ - (id)tableView:(NSTableView *)aTableView objectValueForTableColumn:(NSTableColu case SpatialKernelType::kCauchy: [paramString appendFormat:@", γ=%.3f", interactionType->if_param2_]; break; + case SpatialKernelType::kStudentsT: + [paramString appendFormat:@", ν=%.3f, σ=%.3f", interactionType->if_param2_, interactionType->if_param3_]; + break; } return [paramString autorelease]; diff --git a/VERSIONS b/VERSIONS index 85b4881b4..b498afd0d 100644 --- a/VERSIONS +++ b/VERSIONS @@ -150,6 +150,7 @@ development head (in the master branch): add smooth() to SpatialMap, to smooth/blur a map; internally, add the SpatialKernel class (not user-visible) add add(), multiply(), subtract(), divide(), power(), exp(), range(), and gridValues() methods to SpatialMap add sampleNearbyPoint() and sampleImprovedNearbyPoint() to SpatialMap, for searching habitat and such + add a t-distribution option for SpatialKernel, to provide a fat-tailed distribution that behaves better than Cauchy version 4.0.1 (Eidos version 3.0.1): diff --git a/core/interaction_type.cpp b/core/interaction_type.cpp index 17e6c9e83..8cbc2b470 100755 --- a/core/interaction_type.cpp +++ b/core/interaction_type.cpp @@ -763,6 +763,8 @@ double InteractionType::CalculateStrengthNoCallbacks(double p_distance) double temp = p_distance / if_param2_; return (if_param1_ / (1.0 + temp * temp)); // fmax / (1+(d/λ)^2) } + case SpatialKernelType::kStudentsT: + return SpatialKernel::tdist(p_distance, if_param1_, if_param2_, if_param3_); // fmax / (1+(d/t)^2/n)^(−(ν+1)/2) } EIDOS_TERMINATION << "ERROR (InteractionType::CalculateStrengthNoCallbacks): (internal error) unexpected SpatialKernelType." << EidosTerminate(); } @@ -2457,7 +2459,6 @@ void InteractionType::BuildSV_Strengths_n_2(SLiM_kdNode *root, double *nd, slim_ } } - // 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) { @@ -2487,6 +2488,35 @@ void InteractionType::BuildSV_Strengths_c_2(SLiM_kdNode *root, double *nd, slim_ } } +// add neighbor strengths of type "t" (SpatialKernelType::kStudentsT : Student's t) to the sparse vector in 2D +void InteractionType::BuildSV_Strengths_t_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); +#ifndef __clang_analyzer__ + double dx = root->x[p_phase] - nd[p_phase]; +#else + double dx = 0.0; +#endif + double dx2 = dx * dx; + + if ((d <= max_distance_sq_) && (root->individual_index_ != p_focal_individual_index)) + { + d = sqrt(d); + p_sparse_vector->AddEntryStrength(root->individual_index_, (sv_value_t)SpatialKernel::tdist(d, if_param1_, if_param2_, if_param3_)); + } + + if (++p_phase >= 2) p_phase = 0; + if (dx > 0) { + if (root->left) BuildSV_Strengths_t_2(root->left, nd, p_focal_individual_index, p_sparse_vector, p_phase); + if (dx2 > max_distance_sq_) return; + if (root->right) BuildSV_Strengths_t_2(root->right, nd, p_focal_individual_index, p_sparse_vector, p_phase); + } else { + if (root->right) BuildSV_Strengths_t_2(root->right, nd, p_focal_individual_index, p_sparse_vector, p_phase); + if (dx2 > max_distance_sq_) return; + if (root->left) BuildSV_Strengths_t_2(root->left, nd, p_focal_individual_index, p_sparse_vector, p_phase); + } +} + void InteractionType::FillSparseVectorForReceiverPresences(SparseVector *sv, Individual *receiver, double *receiver_position, Subpopulation *exerter_subpop, InteractionsData &exerter_subpop_data) { #if DEBUG @@ -2726,6 +2756,7 @@ void InteractionType::FillSparseVectorForReceiverStrengths(SparseVector *sv, Ind 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; + case SpatialKernelType::kStudentsT: BuildSV_Strengths_t_2(exerter_subpop_data.kd_root_, receiver_position, excluded_index, sv, 0); break; default: EIDOS_TERMINATION << "ERROR (InteractionType::FillSparseVectorForReceiverStrengths): (internal error) unoptimized SpatialKernelType value." << EidosTerminate(); } @@ -2820,6 +2851,16 @@ void InteractionType::FillSparseVectorForReceiverStrengths(SparseVector *sv, Ind } break; } + case SpatialKernelType::kStudentsT: + { + for (uint32_t col_iter = 0; col_iter < nnz; ++col_iter) + { + sv_value_t distance = values[col_iter]; + + values[col_iter] = (sv_value_t)SpatialKernel::tdist(distance, if_param1_, if_param2_, if_param3_); + } + break; + } default: { // should never be hit, but this is the base case diff --git a/core/interaction_type.h b/core/interaction_type.h index 0a81903f0..21143cbe4 100644 --- a/core/interaction_type.h +++ b/core/interaction_type.h @@ -138,7 +138,7 @@ class InteractionType : public EidosDictionaryUnretained slim_usertag_t tag_value_ = SLIM_TAG_UNSET_VALUE; // a user-defined tag value 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 if_param1_, if_param2_, if_param3_; // 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_ std::map data_; // cached data for the interaction, for each "exerter" subpopulation @@ -197,6 +197,7 @@ class InteractionType : public EidosDictionaryUnretained void BuildSV_Strengths_e_2(SLiM_kdNode *root, double *nd, slim_popsize_t p_focal_individual_index, SparseVector *p_sparse_vector, int p_phase); void BuildSV_Strengths_n_2(SLiM_kdNode *root, double *nd, slim_popsize_t p_focal_individual_index, SparseVector *p_sparse_vector, int p_phase); void BuildSV_Strengths_c_2(SLiM_kdNode *root, double *nd, slim_popsize_t p_focal_individual_index, SparseVector *p_sparse_vector, int p_phase); + void BuildSV_Strengths_t_2(SLiM_kdNode *root, double *nd, slim_popsize_t p_focal_individual_index, SparseVector *p_sparse_vector, int p_phase); int CountNeighbors_1(SLiM_kdNode *root, double *nd, slim_popsize_t p_focal_individual_index); int CountNeighbors_2(SLiM_kdNode *root, double *nd, slim_popsize_t p_focal_individual_index, int p_phase); diff --git a/core/spatial_kernel.cpp b/core/spatial_kernel.cpp index 0c82c2a5c..653a00509 100644 --- a/core/spatial_kernel.cpp +++ b/core/spatial_kernel.cpp @@ -33,6 +33,7 @@ std::ostream& operator<<(std::ostream& p_out, SpatialKernelType p_kernel_type) case SpatialKernelType::kExponential: p_out << gStr_e; break; case SpatialKernelType::kNormal: p_out << gEidosStr_n; break; case SpatialKernelType::kCauchy: p_out << gEidosStr_c; break; + case SpatialKernelType::kStudentsT: p_out << gEidosStr_t; break; } return p_out; @@ -107,8 +108,13 @@ SpatialKernel::SpatialKernel(int p_dimensionality, double p_maxDistance, const s k_type = SpatialKernelType::kCauchy; expected_k_param_count = (p_expect_max_density ? 2 : 1); } + else if (k_type_string.compare(gEidosStr_t) == 0) + { + k_type = SpatialKernelType::kStudentsT; + expected_k_param_count = (p_expect_max_density ? 3 : 2); + } else - EIDOS_TERMINATION << "ERROR (SpatialKernel::SpatialKernel): spatial kernel functionType \"" << k_type_string << "\" must be \"f\", \"l\", \"e\", \"n\", or \"c\"." << EidosTerminate(); + EIDOS_TERMINATION << "ERROR (SpatialKernel::SpatialKernel): spatial kernel functionType \"" << k_type_string << "\" must be \"f\", \"l\", \"e\", \"n\", \"c\", or \"t\"." << EidosTerminate(); if ((dimensionality_ == 0) && (k_type != SpatialKernelType::kFixed)) EIDOS_TERMINATION << "ERROR (SpatialKernel::SpatialKernel): spatial kernel functionType 'f' is required for non-spatial interactions." << EidosTerminate(); @@ -154,12 +160,20 @@ SpatialKernel::SpatialKernel(int p_dimensionality, double p_maxDistance, const s if (k_parameters[1] <= 0.0) EIDOS_TERMINATION << "ERROR (SpatialKernel::SpatialKernel): spatial kernel type \"c\" must have a scale parameter > 0." << EidosTerminate(); break; + case SpatialKernelType::kStudentsT: + // nu can range from -inf to +inf but must be greater than the dimensionality minus one; scale (sigma) must be >= 0 + if (k_parameters[1] <= p_dimensionality - 1) + EIDOS_TERMINATION << "ERROR (SpatialKernel::SpatialKernel): spatial kernel type \"t\" must have a degrees of freedom parameter that is greater than the kernel dimensionality minus one." << EidosTerminate(); + if (k_parameters[2] < 0.0) + EIDOS_TERMINATION << "ERROR (SpatialKernel::SpatialKernel): spatial kernel type \"t\" 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); + kernel_param3_ = ((k_parameters.size() >= 3) ? k_parameters[2] : 0.0); n_2param2sq_ = (kernel_type_ == SpatialKernelType::kNormal) ? (2.0 * kernel_param2_ * kernel_param2_) : 0.0; } @@ -316,6 +330,8 @@ double SpatialKernel::DensityForDistance(double p_distance) double temp = p_distance / kernel_param2_; return (kernel_param1_ / (1.0 + temp * temp)); // fmax / (1+(d/λ)^2) } + case SpatialKernelType::kStudentsT: + return SpatialKernel::tdist(p_distance, kernel_param1_, kernel_param2_, kernel_param3_); // fmax / (1+(d/t)^2/n)^(−(ν+1)/2) } EIDOS_TERMINATION << "ERROR (SpatialKernel::DensityForDistance): (internal error) unexpected SpatialKernelType value." << EidosTerminate(); } @@ -326,25 +342,59 @@ void SpatialKernel::DrawDisplacement_S1(double *displacement) // Note that we could be going either plus or minus from the center gsl_rng *rng = EIDOS_GSL_RNG(omp_get_thread_num()); - if (kernel_type_ == SpatialKernelType::kFixed) - { - displacement[0] = Eidos_rng_uniform(rng) * 2 * max_distance_ - max_distance_; - } - else if (kernel_type_ == SpatialKernelType::kNormal) - { - double d; - - do { - d = gsl_ran_gaussian(rng, kernel_param2_); - } while (d > max_distance_); - - displacement[0] = d; - } - else + switch (kernel_type_) { - // Other distributions are of unclear utility, since draws may cluster at the max distance; this is - // particularly bad for the Cauchy, because the area under it out to infinity is infinite for D>1 - EIDOS_TERMINATION << "ERROR (SpatialKernel::DrawDisplacement_S1): 1D case not implemented." << EidosTerminate(); + case SpatialKernelType::kFixed: + { + displacement[0] = Eidos_rng_uniform(rng) * 2 * max_distance_ - max_distance_; + return; + } + case SpatialKernelType::kLinear: + { + double d = (1 - sqrt(Eidos_rng_uniform(rng))) * max_distance_; + + displacement[0] = d; + return; + } + case SpatialKernelType::kExponential: + { + double d; + + do { + d = gsl_ran_exponential(rng, kernel_param2_); + } while (d > max_distance_); + + displacement[0] = d; + return; + } + case SpatialKernelType::kNormal: + { + double d; + + do { + d = gsl_ran_gaussian(rng, kernel_param2_); + } while (d > max_distance_); + + displacement[0] = d; + return; + } + case SpatialKernelType::kStudentsT: + { + double d; + + do { + d = gsl_ran_tdist(rng, kernel_param2_) * kernel_param3_; + } while (d > max_distance_); + + displacement[0] = d; + return; + } + default: + { + // Other distributions are of unclear utility, since draws may cluster at the max distance; this is + // particularly bad for the Cauchy, because the area under it out to infinity is infinite for D>1 + EIDOS_TERMINATION << "ERROR (SpatialKernel::DrawDisplacement_S1): kernel type not supported." << EidosTerminate(); + } } } @@ -354,30 +404,73 @@ void SpatialKernel::DrawDisplacement_S2(double *displacement) // Note that we could be going in any direction from the center gsl_rng *rng = EIDOS_GSL_RNG(omp_get_thread_num()); - if (kernel_type_ == SpatialKernelType::kFixed) - { - double theta = Eidos_rng_uniform(rng) * 2 * M_PI; - double d = sqrt(Eidos_rng_uniform(rng)) * max_distance_; - displacement[0] = cos(theta) * d; - displacement[1] = sin(theta) * d; - } - else if (kernel_type_ == SpatialKernelType::kNormal) - { - double d1, d2; - - do { - d1 = gsl_ran_gaussian(rng, kernel_param2_); - d2 = gsl_ran_gaussian(rng, kernel_param2_); - } while (sqrt(d1*d1 + d2*d2) > max_distance_); - - displacement[0] = d1; - displacement[1] = d2; - } - else + switch (kernel_type_) { - // Other distributions are of unclear utility, since draws may cluster at the max distance; this is - // particularly bad for the Cauchy, because the area under it out to infinity is infinite for D>1 - EIDOS_TERMINATION << "ERROR (SpatialKernel::DrawDisplacement_S2): 2D case not implemented." << EidosTerminate(); + case SpatialKernelType::kFixed: + { + double theta = Eidos_rng_uniform(rng) * 2 * M_PI; + double d = sqrt(Eidos_rng_uniform(rng)) * max_distance_; + displacement[0] = cos(theta) * d; + displacement[1] = sin(theta) * d; + return; + } + case SpatialKernelType::kLinear: + { + double theta = Eidos_rng_uniform(rng) * 2 * M_PI; + double d = gsl_ran_beta(rng, 2.0, 2.0) * max_distance_; + displacement[0] = cos(theta) * d; + displacement[1] = sin(theta) * d; + return; + } + case SpatialKernelType::kExponential: + { + double d; + + do { + d = gsl_ran_gamma(rng, 2.0, 1.0 / kernel_param2_); + } while (d > max_distance_); + + double theta = Eidos_rng_uniform(rng) * 2 * M_PI; + + displacement[0] = cos(theta) * d; + displacement[1] = sin(theta) * d; + return; + } + case SpatialKernelType::kNormal: + { + double d1, d2; + + do { + d1 = gsl_ran_gaussian(rng, kernel_param2_); + d2 = gsl_ran_gaussian(rng, kernel_param2_); + } while (sqrt(d1*d1 + d2*d2) > max_distance_); + + displacement[0] = d1; + displacement[1] = d2; + return; + } + case SpatialKernelType::kStudentsT: + { + // df (nu) is kernel_param2_, scale is kernel_param3_ + double d; + + do { + double x = 0.5 + abs(Eidos_rng_uniform(rng) - 0.5); + d = sqrt(std::max(0.0, kernel_param2_ * (pow(2.0 - 2.0 * x, -2.0 / (kernel_param2_ - 1.0)) - 1.0))) * kernel_param3_; + } while (d > max_distance_); + + double theta = Eidos_rng_uniform(rng) * 2 * M_PI; + + displacement[0] = cos(theta) * d; + displacement[1] = sin(theta) * d; + return; + } + default: + { + // Other distributions are of unclear utility, since draws may cluster at the max distance; this is + // particularly bad for the Cauchy, because the area under it out to infinity is infinite for D>1 + EIDOS_TERMINATION << "ERROR (SpatialKernel::DrawDisplacement_S2): kernel type not supported." << EidosTerminate(); + } } } @@ -387,37 +480,73 @@ void SpatialKernel::DrawDisplacement_S3(double *displacement) // Note that we could be going in any direction from the center gsl_rng *rng = EIDOS_GSL_RNG(omp_get_thread_num()); - if (kernel_type_ == SpatialKernelType::kFixed) - { - double dx = gsl_ran_gaussian(rng, 1.0); - double dy = gsl_ran_gaussian(rng, 1.0); - double dz = gsl_ran_gaussian(rng, 1.0); - double sphere_dist = sqrt(dx*dx + dy*dy + dz*dz); - double d = pow(Eidos_rng_uniform(rng), 1/3.0) * max_distance_; - - displacement[0] = dx * d / sphere_dist; - displacement[1] = dy * d / sphere_dist; - displacement[2] = dz * d / sphere_dist; - } - else if (kernel_type_ == SpatialKernelType::kNormal) - { - double d1, d2, d3; - - do { - d1 = gsl_ran_gaussian(rng, kernel_param2_); - d2 = gsl_ran_gaussian(rng, kernel_param2_); - d3 = gsl_ran_gaussian(rng, kernel_param2_); - } while (sqrt(d1*d1 + d2*d2 + d3*d3) > max_distance_); - - displacement[0] = d1; - displacement[1] = d2; - displacement[2] = d3; - } - else + switch (kernel_type_) { - // Other distributions are of unclear utility, since draws may cluster at the max distance; this is - // particularly bad for the Cauchy, because the area under it out to infinity is infinite for D>1 - EIDOS_TERMINATION << "ERROR (SpatialKernel::DrawDisplacement_S3): 3D case not implemented." << EidosTerminate(); + case SpatialKernelType::kFixed: + { + double dx = gsl_ran_gaussian(rng, 1.0); + double dy = gsl_ran_gaussian(rng, 1.0); + double dz = gsl_ran_gaussian(rng, 1.0); + double sphere_dist = sqrt(dx*dx + dy*dy + dz*dz); + double d = pow(Eidos_rng_uniform(rng), 1/3.0) * max_distance_; + + displacement[0] = dx * d / sphere_dist; + displacement[1] = dy * d / sphere_dist; + displacement[2] = dz * d / sphere_dist; + return; + } + case SpatialKernelType::kLinear: + { + double dx = gsl_ran_gaussian(rng, 1.0); + double dy = gsl_ran_gaussian(rng, 1.0); + double dz = gsl_ran_gaussian(rng, 1.0); + double sphere_dist = sqrt(dx*dx + dy*dy + dz*dz); + double d = gsl_ran_beta(rng, 3.0, 2.0) * max_distance_; + + displacement[0] = dx * d / sphere_dist; + displacement[1] = dy * d / sphere_dist; + displacement[2] = dz * d / sphere_dist; + return; + } + case SpatialKernelType::kExponential: + { + double dx = gsl_ran_gaussian(rng, 1.0); + double dy = gsl_ran_gaussian(rng, 1.0); + double dz = gsl_ran_gaussian(rng, 1.0); + double sphere_dist = sqrt(dx*dx + dy*dy + dz*dz); + double d; + + do { + d = gsl_ran_gamma(rng, 3.0, 1.0 / kernel_param2_) * max_distance_; + } while (d > max_distance_); + + displacement[0] = dx * d / sphere_dist; + displacement[1] = dy * d / sphere_dist; + displacement[2] = dz * d / sphere_dist; + return; + } + case SpatialKernelType::kNormal: + { + double d1, d2, d3; + + do { + d1 = gsl_ran_gaussian(rng, kernel_param2_); + d2 = gsl_ran_gaussian(rng, kernel_param2_); + d3 = gsl_ran_gaussian(rng, kernel_param2_); + } while (sqrt(d1*d1 + d2*d2 + d3*d3) > max_distance_); + + displacement[0] = d1; + displacement[1] = d2; + displacement[2] = d3; + return; + } + case SpatialKernelType::kStudentsT: // punting for now + default: + { + // Other distributions are of unclear utility, since draws may cluster at the max distance; this is + // particularly bad for the Cauchy, because the area under it out to infinity is infinite for D>1 + EIDOS_TERMINATION << "ERROR (SpatialKernel::DrawDisplacement_S3): kernel type not supported." << EidosTerminate(); + } } } diff --git a/core/spatial_kernel.h b/core/spatial_kernel.h index b97fb23cd..5a7ca55a1 100644 --- a/core/spatial_kernel.h +++ b/core/spatial_kernel.h @@ -37,11 +37,12 @@ // 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 + kFixed = 0, // "f" + kLinear, // "l" + kExponential, // "e" + kNormal, // "n" + kCauchy, // "c" + kStudentsT, // "t" }; std::ostream& operator<<(std::ostream& p_out, SpatialKernelType p_kernel_type); @@ -62,13 +63,17 @@ class SpatialKernel 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 kernel_param1_, kernel_param2_, kernel_param3_; // 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 + // calculate t-distribution PDF values in our fashion, for which the function is normalized to a maximum value + // we don't use the GSL for this, because it does two gamma-function calculations that we don't need (they normalize away) + static inline double tdist(double x, double max, double nu, double tau) { double x_over_tau = x / tau; return max / pow(1.0 + x_over_tau * x_over_tau / nu, -(nu + 1.0) / 2.0); }; + public: SpatialKernel(const SpatialKernel&) = delete; // no copying SpatialKernel& operator=(const SpatialKernel&) = delete; // no copying diff --git a/eidos/eidos_globals.cpp b/eidos/eidos_globals.cpp index 83aaf7cc6..de376a674 100644 --- a/eidos/eidos_globals.cpp +++ b/eidos/eidos_globals.cpp @@ -3762,6 +3762,7 @@ const std::string &gEidosStr_end = EidosRegisteredString("end", gEidosID_end); const std::string &gEidosStr_weights = EidosRegisteredString("weights", gEidosID_weights); const std::string &gEidosStr_range = EidosRegisteredString("range", gEidosID_range); const std::string &gEidosStr_c = EidosRegisteredString("c", gEidosID_c); +const std::string &gEidosStr_t = EidosRegisteredString("t", gEidosID_t); const std::string &gEidosStr_n = EidosRegisteredString("n", gEidosID_n); const std::string &gEidosStr_s = EidosRegisteredString("s", gEidosID_s); const std::string &gEidosStr_x = EidosRegisteredString("x", gEidosID_x); diff --git a/eidos/eidos_globals.h b/eidos/eidos_globals.h index bd5386c8b..3d45ec87e 100644 --- a/eidos/eidos_globals.h +++ b/eidos/eidos_globals.h @@ -974,6 +974,7 @@ extern const std::string &gEidosStr_end; extern const std::string &gEidosStr_weights; extern const std::string &gEidosStr_range; extern const std::string &gEidosStr_c; +extern const std::string &gEidosStr_t; extern const std::string &gEidosStr_n; extern const std::string &gEidosStr_s; extern const std::string &gEidosStr_x; @@ -1100,6 +1101,7 @@ enum _EidosGlobalStringID : uint32_t gEidosID_weights, gEidosID_range, gEidosID_c, + gEidosID_t, gEidosID_n, gEidosID_s, gEidosID_x, diff --git a/gsl/gsl.pro b/gsl/gsl.pro index 536601883..a12e97a2b 100644 --- a/gsl/gsl.pro +++ b/gsl/gsl.pro @@ -75,6 +75,7 @@ SOURCES += \ randist/beta.c \ randist/binomial_tpe.c \ randist/cauchy.c \ + randist/chisq.c \ randist/discrete.c \ randist/exponential.c \ randist/fdist.c \ @@ -88,6 +89,7 @@ SOURCES += \ randist/nbinomial.c \ randist/poisson.c \ randist/shuffle.c \ + randist/tdist.c \ randist/weibull.c \ rng/inline.c \ rng/mt.c \ diff --git a/gsl/randist/chisq.c b/gsl/randist/chisq.c new file mode 100644 index 000000000..c540e7a18 --- /dev/null +++ b/gsl/randist/chisq.c @@ -0,0 +1,61 @@ +/* randist/chisq.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough + * + * This program 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. + * + * This program 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 this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include "config.h" +#include +#include "gsl_sf_gamma.h" +#include "gsl_rng.h" +#include "gsl_randist.h" + +/* The chisq distribution has the form + + p(x) dx = (1/(2*Gamma(nu/2))) (x/2)^(nu/2 - 1) exp(-x/2) dx + + for x = 0 ... +infty */ + +double +gsl_ran_chisq (const gsl_rng * r, const double nu) +{ + double chisq = 2 * gsl_ran_gamma (r, nu / 2, 1.0); + return chisq; +} + +double +gsl_ran_chisq_pdf (const double x, const double nu) +{ + if (x < 0) + { + return 0 ; + } + else + { + if(nu == 2.0) + { + return exp(-x/2.0) / 2.0; + } + else + { + double p; + double lngamma = gsl_sf_lngamma (nu / 2); + + p = exp ((nu / 2 - 1) * log (x/2) - x/2 - lngamma) / 2; + return p; + } + } +} diff --git a/gsl/randist/tdist.c b/gsl/randist/tdist.c new file mode 100644 index 000000000..ece3a33da --- /dev/null +++ b/gsl/randist/tdist.c @@ -0,0 +1,80 @@ +/* randist/tdist.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough + * + * This program 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. + * + * This program 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 this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include "config.h" +#include +#include "gsl_math.h" +#include "gsl_sf_gamma.h" +#include "gsl_rng.h" +#include "gsl_randist.h" + +/* The t-distribution has the form + + p(x) dx = (Gamma((nu + 1)/2)/(sqrt(pi nu) Gamma(nu/2)) + * (1 + (x^2)/nu)^-((nu + 1)/2) dx + + The method used here is the one described in Knuth */ + +double +gsl_ran_tdist (const gsl_rng * r, const double nu) +{ + if (nu <= 2) + { + double Y1 = gsl_ran_ugaussian (r); + double Y2 = gsl_ran_chisq (r, nu); + + double t = Y1 / sqrt (Y2 / nu); + + return t; + } + else + { + double Y1, Y2, Z, t; + do + { + Y1 = gsl_ran_ugaussian (r); + Y2 = gsl_ran_exponential (r, 1 / (nu/2 - 1)); + + Z = Y1 * Y1 / (nu - 2); + } + while (1 - Z < 0 || exp (-Y2 - Z) > (1 - Z)); + + /* Note that there is a typo in Knuth's formula, the line below + is taken from the original paper of Marsaglia, Mathematics of + Computation, 34 (1980), p 234-256 */ + + t = Y1 / sqrt ((1 - 2 / nu) * (1 - Z)); + return t; + } +} + +double +gsl_ran_tdist_pdf (const double x, const double nu) +{ + double p; + + double lg1 = gsl_sf_lngamma (nu / 2); + double lg2 = gsl_sf_lngamma ((nu + 1) / 2); + + p = ((exp (lg2 - lg1) / sqrt (M_PI * nu)) + * pow ((1 + x * x / nu), -(nu + 1) / 2)); + return p; +} + +