Allow configuring sampling and non-self-dispersal strategy

docs/simulate.ron

@@ -300,11 +300,60 @@
                     /* downscale factor along the x-axis of the habitat
                      * expressed either as a power-of-two integer
                      *  or as a base-two scientific notation string */
-                    x: (u32 = 2^{1..=15}) or ("1B{1..=15}"),
+                    x: (u32 = 2^{1..=16}) or ("1B{1..=16}"),
                     /* downscale factor along the y-axis of the habitat
                      * expressed either as a power-of-two integer
                      *  or as a base-two scientific notation string */
-                    y: (u32 = 2^{1..=15}) or ("1B{1..=15}"),
+                    y: (u32 = 2^{1..=16}) or ("1B{1..=16}"),
+                    /* number of random samples used to emperically determine
+                     *  the self-dispersal probability (and the non-self-
+                     *  dispersal target probability distribution,
+                     *  if precomputing is used, see non_self_dispersal below)
+                     * optional, default = 2^22 */
+                    samples: (0 < u32),
+                    /* strategy for sampling non-self-dispersal iff non-self-
+                     *  dispersal must be sampled separately, e.g. for the
+                     *  EventSkipping algorithm
+                     * Note: this option has no effect when normal dispersal
+                     *       sampling is used 
+                     * optional, default = Accurate */
+                    non_self_dispersal: (
+                        /* use accurate rejection sampling to sample non-self-
+                         *  dispersals
+                         * rejection sampling may be very slow if the
+                         *  probability of self-dispersal is very high */
+                      | Accurate
+                        /* choose between accurate but slow rejection sampling
+                         *  and the approximate but fast precomputing based
+                         *  on the probability of non-self-dispersal
+                         * refer to the Accurate and PrecomputeApproximation
+                         *  strategies above and below for more information on
+                         *  their merits and when to choose which */
+                      | UnlikelyApproximation(
+                            /* threshold for the non-self-dispersal probability
+                             *  for which rejection sampling is chosen
+                             * choose rejection sampling if
+                             *  P(non-self-dispersal) >= likelihood
+                             * choose precomputing if
+                             *  P(non-self-dispersal) < likelihood
+                             * likelihood = 0.0 is equivalent to Accurate
+                             * likelihood = 1.0 is equivalent to
+                             *  PrecomputeApproximation */
+                            likelihood: (0.0 <= f64 <= 1.0),
+                        )
+                        /* precompute the non-self-dispersal target distribution
+                         *  using the number of samples given above
+                         * precomputing the distribution is an approximation
+                         *  that only stores a discrete set of dispersal targets
+                         *  (which may exclude some valid targets) with
+                         *  approximate dispersal probabilities
+                         * precomputing only be chosen when this approximation
+                         *  can sufficiently represent the dispersal
+                         * sampling a precomputed dispersal distribution may be
+                         *  significantly faster than using rejection sampling
+                         *  when the probability of self-dispersal is high */
+                      | PrecomputeApproximation
+                    )
                 )
             )
         )

necsim/core/bond/src/off_by_one_u32.rs

@@ -1,4 +1,8 @@
-use core::{convert::TryFrom, fmt, num::NonZeroU64};
+use core::{
+    convert::TryFrom,
+    fmt,
+    num::{NonZeroU32, NonZeroU64},
+};
 
 use serde::{Deserialize, Deserializer, Serialize};
 
@@ -56,6 +60,11 @@ impl OffByOneU32 {
         (self.0 as u64) + 1_u64
     }
 
+    #[must_use]
+    pub const fn sub_one(self) -> u32 {
+        self.0
+    }
+
     #[must_use]
     pub const fn add_incl(self, other: u32) -> u32 {
         other.wrapping_add(self.0)
@@ -82,6 +91,12 @@ impl OffByOneU32 {
     }
 }
 
+impl From<NonZeroU32> for OffByOneU32 {
+    fn from(value: NonZeroU32) -> Self {
+        Self(value.get() - 1)
+    }
+}
+
 impl TryFrom<u64> for OffByOneU32 {
     type Error = OffByOneU32Error;
 

necsim/core/bond/src/off_by_one_u64.rs

@@ -63,6 +63,11 @@ impl OffByOneU64 {
         (self.0 as u128) + 1_u128
     }
 
+    #[must_use]
+    pub const fn sub_one(self) -> u64 {
+        self.0
+    }
+
     #[must_use]
     pub const fn add_incl(self, other: u64) -> u64 {
         other.wrapping_add(self.0)

necsim/core/src/cogs/rng.rs

@@ -1,14 +1,15 @@
 use core::{
     convert::AsMut,
     default::Default,
-    num::{NonZeroU128, NonZeroU32, NonZeroU64, NonZeroUsize},
+    num::{NonZeroU128, NonZeroUsize},
     ptr::copy_nonoverlapping,
 };
 
 use serde::{de::DeserializeOwned, Serialize};
 
 use necsim_core_bond::{
-    ClosedOpenUnitF64, ClosedUnitF64, NonNegativeF64, OpenClosedUnitF64, PositiveF64,
+    ClosedOpenUnitF64, ClosedUnitF64, NonNegativeF64, OffByOneU32, OffByOneU64, OpenClosedUnitF64,
+    PositiveF64,
 };
 
 use crate::{
@@ -111,8 +112,8 @@ pub trait RngSampler<M: MathsCore>: RngCore<M> {
 
     #[must_use]
     #[inline]
-    #[debug_ensures(ret < length.get(), "samples U(0, length - 1)")]
-    fn sample_index_u32(&mut self, length: NonZeroU32) -> u32 {
+    #[debug_ensures(u64::from(ret) < length.get(), "samples U(0, length - 1)")]
+    fn sample_index_u32(&mut self, length: OffByOneU32) -> u32 {
         // attributes on expressions are experimental
         // see https://github.com/rust-lang/rust/issues/15701
         #[allow(
@@ -121,15 +122,15 @@ pub trait RngSampler<M: MathsCore>: RngCore<M> {
             clippy::cast_sign_loss
         )]
         let index =
-            M::floor(self.sample_uniform_closed_open().get() * f64::from(length.get())) as u32;
+            M::floor(self.sample_uniform_closed_open().get() * (length.get() as f64)) as u32;
         // Safety in case of f64 rounding errors
-        index.min(length.get() - 1)
+        index.min(length.sub_one())
     }
 
     #[must_use]
     #[inline]
-    #[debug_ensures(ret < length.get(), "samples U(0, length - 1)")]
-    fn sample_index_u64(&mut self, length: NonZeroU64) -> u64 {
+    #[debug_ensures(u128::from(ret) < length.get(), "samples U(0, length - 1)")]
+    fn sample_index_u64(&mut self, length: OffByOneU64) -> u64 {
         // attributes on expressions are experimental
         // see https://github.com/rust-lang/rust/issues/15701
         #[allow(
@@ -140,7 +141,7 @@ pub trait RngSampler<M: MathsCore>: RngCore<M> {
         let index =
             M::floor(self.sample_uniform_closed_open().get() * (length.get() as f64)) as u64;
         // Safety in case of f64 rounding errors
-        index.min(length.get() - 1)
+        index.min(length.sub_one())
     }
 
     #[must_use]

necsim/core/src/landscape/extent.rs

@@ -1,8 +1,8 @@
-use necsim_core_bond::OffByOneU32;
+use necsim_core_bond::{OffByOneU32, OffByOneU64};
 
 use super::Location;
 
-#[allow(clippy::module_name_repetitions)]
+#[allow(clippy::module_name_repetitions, clippy::unsafe_derive_deserialize)]
 #[derive(PartialEq, Eq, Clone, Debug, serde::Deserialize, serde::Serialize, TypeLayout)]
 #[cfg_attr(feature = "cuda", derive(rust_cuda::lend::LendRustToCuda))]
 #[cfg_attr(feature = "cuda", cuda(ignore))]
@@ -40,6 +40,14 @@ impl LandscapeExtent {
         self.height
     }
 
+    #[must_use]
+    pub const fn area(&self) -> OffByOneU64 {
+        // Safety: {1..=2^32} * {1..=2^32} = {1..=2^64}
+        unsafe {
+            OffByOneU64::new_unchecked((self.width.get() as u128) * (self.height.get() as u128))
+        }
+    }
+
     #[must_use]
     pub const fn contains(&self, location: &Location) -> bool {
         location.x() >= self.origin.x()

necsim/impls/no-std/src/cogs/dispersal_sampler/almost_infinite_clark2dt.rs → necsim/impls/no-std/src/cogs/dispersal_sampler/almost_infinite/clark2dt.rs

@@ -132,7 +132,7 @@ impl<M: MathsCore, G: RngCore<M>> SeparableDispersalSampler<M, AlmostInfiniteHab
     }
 }
 
-#[allow(clippy::doc_markdown)]
+#[allow(clippy::doc_markdown, clippy::module_inception)]
 /// Clark2Dt dispersal:
 ///
 /// Clark, J.S., Silman, M., Kern, R., Macklin, E. and HilleRisLambers, J.