diff --git a/necsim/core/src/cogs/dispersal_sampler.rs b/necsim/core/src/cogs/dispersal_sampler.rs index 60e6b7b88..00f7e01ed 100644 --- a/necsim/core/src/cogs/dispersal_sampler.rs +++ b/necsim/core/src/cogs/dispersal_sampler.rs @@ -15,8 +15,8 @@ pub trait DispersalSampler, G: RngCore>: crate::cogs::Backup + core::fmt::Debug { #[must_use] - #[debug_requires(habitat.contains(location), "location is inside habitat")] - #[debug_ensures(old(habitat).contains(&ret), "target is inside habitat")] + #[debug_requires(habitat.is_location_habitable(location), "location is habitable")] + #[debug_ensures(old(habitat).is_location_habitable(&ret), "target is habitable")] fn sample_dispersal_from_location( &self, location: &Location, @@ -33,8 +33,8 @@ pub trait SeparableDispersalSampler, G: RngCore>: DispersalSampler { #[must_use] - #[debug_requires(habitat.contains(location), "location is inside habitat")] - #[debug_ensures(old(habitat).contains(&ret), "target is inside habitat")] + #[debug_requires(habitat.is_location_habitable(location), "location is habitable")] + #[debug_ensures(old(habitat).is_location_habitable(&ret), "target is habitable")] #[debug_ensures(&ret != location, "disperses to a different location")] fn sample_non_self_dispersal_from_location( &self, @@ -44,7 +44,7 @@ pub trait SeparableDispersalSampler, G: RngCore>: ) -> Location; #[must_use] - #[debug_requires(habitat.contains(location), "location is inside habitat")] + #[debug_requires(habitat.is_location_habitable(location), "location is habitable")] fn get_self_dispersal_probability_at_location( &self, location: &Location, diff --git a/necsim/core/src/cogs/habitat.rs b/necsim/core/src/cogs/habitat.rs index b5f27721e..105ed41b6 100644 --- a/necsim/core/src/cogs/habitat.rs +++ b/necsim/core/src/cogs/habitat.rs @@ -22,8 +22,15 @@ pub trait Habitat: crate::cogs::Backup + core::fmt::Debug + Sized fn get_extent(&self) -> &LandscapeExtent; #[must_use] - fn contains(&self, location: &Location) -> bool { - self.get_extent().contains(location) + fn is_location_habitable(&self, location: &Location) -> bool { + self.get_extent().contains(location) && self.get_habitat_at_location(location) > 0_u32 + } + + #[must_use] + fn is_indexed_location_habitable(&self, indexed_location: &IndexedLocation) -> bool { + self.get_extent().contains(indexed_location.location()) + && (indexed_location.index() + < self.get_habitat_at_location(indexed_location.location())) } #[must_use] diff --git a/necsim/core/src/cogs/lineage_store.rs b/necsim/core/src/cogs/lineage_store.rs index 8d589cf69..afc0d319d 100644 --- a/necsim/core/src/cogs/lineage_store.rs +++ b/necsim/core/src/cogs/lineage_store.rs @@ -32,8 +32,8 @@ pub trait LocallyCoherentLineageStore>: { #[must_use] #[debug_requires( - habitat.contains(indexed_location.location()), - "indexed location is inside habitat" + habitat.is_indexed_location_habitable(indexed_location), + "indexed location is habitable" )] fn get_global_lineage_reference_at_indexed_location( &self, @@ -42,8 +42,8 @@ pub trait LocallyCoherentLineageStore>: ) -> Option<&GlobalLineageReference>; #[debug_requires( - habitat.contains(lineage.indexed_location.location()), - "indexed location is inside habitat" + habitat.is_indexed_location_habitable(&lineage.indexed_location), + "indexed location is habitable" )] #[debug_ensures(self.get_lineage_for_local_reference( &ret @@ -68,9 +68,10 @@ pub trait LocallyCoherentLineageStore>: #[debug_requires(self.get_lineage_for_local_reference( &reference ).is_some(), "lineage is active")] - #[debug_ensures(old(habitat).contains( - ret.indexed_location.location() - ), "prior location is inside habitat")] + #[debug_ensures( + old(habitat).is_indexed_location_habitable(&ret.indexed_location), + "prior indexed location is habitable" + )] #[debug_ensures(self.get_lineage_for_local_reference( &old(unsafe { crate::cogs::Backup::backup_unchecked(&reference) }) ).is_none(), "lineage was deactivated")] @@ -104,7 +105,10 @@ pub trait GloballyCoherentLineageStore>: fn iter_active_locations(&self, habitat: &H) -> Self::LocationIterator<'_>; #[must_use] - #[debug_requires(habitat.contains(location), "location is inside habitat")] + #[debug_requires( + habitat.is_location_habitable(location), + "location is habitable" + )] fn get_local_lineage_references_at_location_unordered( &self, location: &Location, diff --git a/necsim/core/src/cogs/speciation_probability.rs b/necsim/core/src/cogs/speciation_probability.rs index 75ef6c6c4..8584228d6 100644 --- a/necsim/core/src/cogs/speciation_probability.rs +++ b/necsim/core/src/cogs/speciation_probability.rs @@ -11,7 +11,10 @@ pub trait SpeciationProbability>: crate::cogs::Backup + core::fmt::Debug { #[must_use] - #[debug_requires(habitat.contains(location), "location is inside habitat")] + #[debug_requires( + habitat.is_location_habitable(location), + "location is habitable" + )] fn get_speciation_probability_at_location( &self, location: &Location, diff --git a/necsim/core/src/cogs/turnover_rate.rs b/necsim/core/src/cogs/turnover_rate.rs index 1763876f0..fce772b01 100644 --- a/necsim/core/src/cogs/turnover_rate.rs +++ b/necsim/core/src/cogs/turnover_rate.rs @@ -11,7 +11,10 @@ pub trait TurnoverRate>: crate::cogs::Backup + core::fmt::Debug { #[must_use] - #[debug_requires(habitat.contains(location), "location is inside habitat")] + #[debug_requires( + habitat.is_location_habitable(location), + "location is habitable" + )] #[debug_ensures( ret != 0.0_f64 || habitat.get_habitat_at_location(location) == 0_u32, "only returns zero if the location is inhabitable" diff --git a/necsim/core/src/event.rs b/necsim/core/src/event.rs index 637bfdace..40108ae85 100644 --- a/necsim/core/src/event.rs +++ b/necsim/core/src/event.rs @@ -12,43 +12,32 @@ use crate::{ }; #[allow(clippy::module_name_repetitions)] -#[derive(Debug, Clone, Serialize, Deserialize)] +#[allow(clippy::unsafe_derive_deserialize)] +#[derive(Debug, Clone, Serialize, Deserialize, TypeLayout)] +#[repr(C)] pub struct PackedEvent { global_lineage_reference: GlobalLineageReference, - prior_time: NonNegativeF64, // time of the previous event - event_time: PositiveF64, // time of this event + // p\e (-) (+) + // (-) N S + // (+) M C + prior_time: f64, // time of the previous event, NonNegativeF64 + event_time: f64, // time of this event, PositiveF64 origin: IndexedLocation, - target: Option, - interaction: LineageInteraction, + target: IndexedLocation, + coalescence: GlobalLineageReference, } impl PackedEvent { #[must_use] #[inline] pub fn event_time(&self) -> PositiveF64 { - self.event_time + unsafe { PositiveF64::new_unchecked(self.event_time.make_positive()) } } } -#[allow(dead_code)] -const EXCESSIVE_OPTION_PACKED_EVENT_ERROR: [(); 1 - { - const ASSERT: bool = - core::mem::size_of::>() == core::mem::size_of::(); - ASSERT -} as usize] = []; - #[allow(dead_code)] const EXCESSIVE_PACKED_EVENT_ERROR: [(); 1 - { - const ASSERT: bool = { - const SPECIATION_SIZE: usize = core::mem::size_of::(); - const DISPERSAL_SIZE: usize = core::mem::size_of::(); - - if SPECIATION_SIZE > DISPERSAL_SIZE { - SPECIATION_SIZE - } else { - DISPERSAL_SIZE - } - } == core::mem::size_of::(); + const ASSERT: bool = core::mem::size_of::() == 56; ASSERT } as usize] = []; @@ -77,16 +66,13 @@ pub struct SpeciationEvent { } #[allow(dead_code)] -const EXCESSIVE_OPTION_SPECIATION_EVENT_ERROR: [(); 1 - { - const ASSERT: bool = - core::mem::size_of::>() == core::mem::size_of::(); +const EXCESSIVE_SPECIATION_EVENT_ERROR: [(); 1 - { + const ASSERT: bool = core::mem::size_of::() == 40; ASSERT } as usize] = []; #[allow(clippy::module_name_repetitions)] -#[allow(clippy::unsafe_derive_deserialize)] -#[derive(Debug, Clone, Serialize, Deserialize, TypeLayout)] -#[repr(C)] +#[derive(Debug, Clone, Serialize, Deserialize)] pub struct DispersalEvent { pub global_lineage_reference: GlobalLineageReference, pub prior_time: NonNegativeF64, // time of the previous event @@ -96,6 +82,12 @@ pub struct DispersalEvent { pub interaction: LineageInteraction, } +#[allow(dead_code)] +const EXCESSIVE_DISPERSAL_EVENT_ERROR: [(); 1 - { + const ASSERT: bool = core::mem::size_of::() == 64; + ASSERT +} as usize] = []; + #[allow(dead_code)] const EXCESSIVE_OPTION_DISPERSAL_EVENT_ERROR: [(); 1 - { const ASSERT: bool = @@ -111,26 +103,62 @@ pub enum TypedEvent { impl From for PackedEvent { fn from(event: SpeciationEvent) -> Self { + // speciation is encoded as p(-), e(+) Self { - global_lineage_reference: event.global_lineage_reference, - prior_time: event.prior_time, - event_time: event.event_time, - origin: event.origin, - target: None, - interaction: LineageInteraction::None, + global_lineage_reference: event.global_lineage_reference.clone(), + prior_time: event.prior_time.get().make_negative(), + event_time: event.event_time.get(), + origin: event.origin.clone(), + target: event.origin, + coalescence: event.global_lineage_reference, } } } impl From for PackedEvent { - fn from(event: DispersalEvent) -> Self { + fn from( + DispersalEvent { + global_lineage_reference, + prior_time, + event_time, + origin, + target, + interaction, + }: DispersalEvent, + ) -> Self { + let prior_time = prior_time.get(); + let event_time = event_time.get(); + + let (coalescence, prior_time, event_time) = match interaction { + LineageInteraction::None => { + // dispersal, no coalescence is encoded as p(-), e(-) + ( + global_lineage_reference.clone(), + prior_time.make_negative(), + event_time.make_negative(), + ) + }, + LineageInteraction::Maybe => { + // dispersal, maybe coalescence is encoded as p(+), e(-) + ( + global_lineage_reference.clone(), + prior_time, + event_time.make_negative(), + ) + }, + LineageInteraction::Coalescence(coalescence) => { + // dispersal, with coalescence is encoded as p(+), e(+) + (coalescence, prior_time, event_time) + }, + }; + Self { - global_lineage_reference: event.global_lineage_reference, - prior_time: event.prior_time, - event_time: event.event_time, - origin: event.origin, - target: Some(event.target), - interaction: event.interaction, + global_lineage_reference, + prior_time, + event_time, + origin, + target, + coalescence, } } } @@ -144,25 +172,70 @@ impl From for PackedEvent { } } +impl From for TypedEvent { + fn from(event: SpeciationEvent) -> Self { + Self::Speciation(event) + } +} + +impl From for TypedEvent { + fn from(event: DispersalEvent) -> Self { + Self::Dispersal(event) + } +} + impl From for TypedEvent { - fn from(event: PackedEvent) -> Self { - #[allow(clippy::option_if_let_else)] - if let Some(target) = event.target { - Self::Dispersal(DispersalEvent { - global_lineage_reference: event.global_lineage_reference, - prior_time: event.prior_time, - event_time: event.event_time, - origin: event.origin, + fn from( + PackedEvent { + global_lineage_reference, + prior_time, + event_time, + origin, + target, + coalescence, + }: PackedEvent, + ) -> Self { + let prior_pos = prior_time.is_sign_positive(); + let event_pos = event_time.is_sign_positive(); + + let prior_time = unsafe { NonNegativeF64::new_unchecked(prior_time.make_positive()) }; + let event_time = unsafe { PositiveF64::new_unchecked(event_time.make_positive()) }; + + match (prior_pos, event_pos) { + // dispersal, no coalescence is encoded as p(-), e(-) + (false, false) => Self::Dispersal(DispersalEvent { + global_lineage_reference, + prior_time, + event_time, + origin, + target, + interaction: LineageInteraction::None, + }), + // speciation is encoded as p(-), e(+) + (false, true) => Self::Speciation(SpeciationEvent { + global_lineage_reference, + prior_time, + event_time, + origin, + }), + // dispersal, maybe coalescence is encoded as p(+), e(-) + (true, false) => Self::Dispersal(DispersalEvent { + global_lineage_reference, + prior_time, + event_time, + origin, + target, + interaction: LineageInteraction::Maybe, + }), + // dispersal, with coalescence is encoded as p(+), e(+) + (true, true) => Self::Dispersal(DispersalEvent { + global_lineage_reference, + prior_time, + event_time, + origin, target, - interaction: event.interaction, - }) - } else { - Self::Speciation(SpeciationEvent { - global_lineage_reference: event.global_lineage_reference, - prior_time: event.prior_time, - event_time: event.event_time, - origin: event.origin, - }) + interaction: LineageInteraction::Coalescence(coalescence), + }), } } } @@ -171,13 +244,13 @@ impl Eq for PackedEvent {} impl PartialEq for PackedEvent { // `Event`s are equal when they have the same `origin`, `event_time`, - // `target` and `interaction` - // (`global_lineage_reference` and `prior_time` are ignored) + // and `target` + // (`global_lineage_reference`, `prior_time`, and `coalescence` are ignored) fn eq(&self, other: &Self) -> bool { self.origin == other.origin - && self.event_time == other.event_time + && unsafe { PositiveF64::new_unchecked(self.event_time) } + == unsafe { PositiveF64::new_unchecked(other.event_time) } && self.target == other.target - && self.interaction == other.interaction } } @@ -186,24 +259,22 @@ impl Ord for PackedEvent { // Order `Event`s in lexicographical order: // (1) event_time /=\ // (2) origin different | events - // (3) target and interaction \=/ + // (3) target \=/ // (4) prior_time parent + offspring // (5) global_lineage_reference - + // (coalescence is ignored) ( - &self.event_time, + &unsafe { PositiveF64::new_unchecked(self.event_time.make_positive()) }, &self.origin, &self.target, - &self.interaction, - &self.prior_time, + &unsafe { NonNegativeF64::new_unchecked(self.prior_time.make_positive()) }, &self.global_lineage_reference, ) .cmp(&( - &other.event_time, + &unsafe { PositiveF64::new_unchecked(other.event_time.make_positive()) }, &other.origin, &other.target, - &other.interaction, - &other.prior_time, + &unsafe { NonNegativeF64::new_unchecked(other.prior_time.make_positive()) }, &other.global_lineage_reference, )) } @@ -217,13 +288,12 @@ impl PartialOrd for PackedEvent { impl Hash for PackedEvent { // `Event`s are equal when they have the same `origin`, `event_time`, - // `target` and `interaction` - // (`global_lineage_reference` and `prior_time` are ignored) + // and `target`. + // (`global_lineage_reference`, `prior_time`, and `coalescence` are ignored) fn hash(&self, state: &mut S) { self.origin.hash(state); - self.event_time.hash(state); + unsafe { PositiveF64::new_unchecked(self.event_time.make_positive()) }.hash(state); self.target.hash(state); - self.interaction.hash(state); } } @@ -241,12 +311,45 @@ impl Eq for DispersalEvent {} impl PartialEq for DispersalEvent { // `SpeciationEvent`s are equal when they have the same `origin`, - // `event_time`, `target` and `interaction` - // (`global_lineage_reference` and `prior_time` are ignored) + // `event_time`, and `target` + // (`global_lineage_reference`, `prior_time`, and `interaction` are ignored) fn eq(&self, other: &Self) -> bool { self.origin == other.origin && self.event_time == other.event_time && self.target == other.target - && self.interaction == other.interaction + } +} + +trait F64 { + fn make_positive(self) -> Self; + fn make_negative(self) -> Self; +} + +impl F64 for f64 { + fn make_positive(self) -> Self { + f64::from_bits(self.to_bits() & 0x7fff_ffff_ffff_ffff_u64) + } + + fn make_negative(self) -> Self { + f64::from_bits(self.to_bits() | 0x8000_0000_0000_0000_u64) + } +} + +#[cfg(test)] +mod tests { + use super::F64; + + #[test] + fn test_make_positive() { + assert_eq!((-0.0_f64).make_positive().to_bits(), (0.0_f64).to_bits()); + assert_eq!((-42.0_f64).make_positive().to_bits(), (42.0_f64).to_bits()); + assert_eq!((24.0_f64).make_positive().to_bits(), (24.0_f64).to_bits()); + } + + #[test] + fn test_make_negative() { + assert_eq!((0.0_f64).make_negative().to_bits(), (-0.0_f64).to_bits()); + assert_eq!((42.0_f64).make_negative().to_bits(), (-42.0_f64).to_bits()); + assert_eq!((-24.0_f64).make_negative().to_bits(), (-24.0_f64).to_bits()); } } diff --git a/necsim/core/src/landscape/extent.rs b/necsim/core/src/landscape/extent.rs index f981184ab..09ac62d76 100644 --- a/necsim/core/src/landscape/extent.rs +++ b/necsim/core/src/landscape/extent.rs @@ -16,11 +16,7 @@ pub struct LandscapeExtent { impl LandscapeExtent { #[must_use] - #[debug_ensures(ret.x() == x, "stores x")] - #[debug_ensures(ret.y() == y, "stores y")] - #[debug_ensures(ret.width() == width, "stores width")] - #[debug_ensures(ret.height() == height, "stores height")] - pub fn new(x: u32, y: u32, width: OffByOneU32, height: OffByOneU32) -> Self { + pub const fn new(x: u32, y: u32, width: OffByOneU32, height: OffByOneU32) -> Self { Self { x, y, @@ -30,27 +26,27 @@ impl LandscapeExtent { } #[must_use] - pub fn x(&self) -> u32 { + pub const fn x(&self) -> u32 { self.x } #[must_use] - pub fn y(&self) -> u32 { + pub const fn y(&self) -> u32 { self.y } #[must_use] - pub fn width(&self) -> OffByOneU32 { + pub const fn width(&self) -> OffByOneU32 { self.width } #[must_use] - pub fn height(&self) -> OffByOneU32 { + pub const fn height(&self) -> OffByOneU32 { self.height } #[must_use] - pub fn contains(&self, location: &Location) -> bool { + pub const fn contains(&self, location: &Location) -> bool { location.x() >= self.x && location.x() <= self.width.add_incl(self.x) && location.y() >= self.y diff --git a/necsim/core/src/landscape/location.rs b/necsim/core/src/landscape/location.rs index 0ba841aeb..c3686e5c6 100644 --- a/necsim/core/src/landscape/location.rs +++ b/necsim/core/src/landscape/location.rs @@ -1,6 +1,4 @@ -use core::num::NonZeroU32; - -use serde::{Deserialize, Deserializer, Serialize, Serializer}; +use serde::{Deserialize, Serialize}; use crate::cogs::Backup; @@ -24,19 +22,17 @@ impl Backup for Location { impl Location { #[must_use] - #[debug_ensures(ret.x() == x, "stores x")] - #[debug_ensures(ret.y() == y, "stores y")] - pub fn new(x: u32, y: u32) -> Self { + pub const fn new(x: u32, y: u32) -> Self { Self { x, y } } #[must_use] - pub fn x(&self) -> u32 { + pub const fn x(&self) -> u32 { self.x } #[must_use] - pub fn y(&self) -> u32 { + pub const fn y(&self) -> u32 { self.y } } @@ -47,8 +43,6 @@ impl From for Location { } } -// IndexedLocation uses a NonZeroU32 index internally to enable same-size -// Option optimisation #[derive( Eq, PartialEq, PartialOrd, Ord, Clone, Hash, Debug, Serialize, Deserialize, TypeLayout, )] @@ -57,59 +51,23 @@ impl From for Location { #[repr(C)] pub struct IndexedLocation { location: Location, - index: LocationIndex, + index: u32, } impl IndexedLocation { #[must_use] - #[debug_ensures( - ret.location.x() == old(location.x()) && ret.location.y() == old(location.y()), - "stores location" - )] - #[debug_ensures(ret.index() == index, "stores index")] - pub fn new(location: Location, index: u32) -> Self { - Self { - location, - index: LocationIndex::new(index), - } + pub const fn new(location: Location, index: u32) -> Self { + Self { location, index } } #[must_use] - pub fn location(&self) -> &Location { + pub const fn location(&self) -> &Location { &self.location } #[must_use] - pub fn index(&self) -> u32 { - self.index.get() - } -} - -#[derive(Eq, PartialEq, PartialOrd, Ord, Clone, Copy, Hash, Debug, TypeLayout)] -#[repr(transparent)] -struct LocationIndex(NonZeroU32); - -impl LocationIndex { - fn new(index: u32) -> Self { - Self(unsafe { NonZeroU32::new_unchecked(index + 1) }) - } - - fn get(self) -> u32 { - self.0.get() - 1 - } -} - -impl Serialize for LocationIndex { - fn serialize(&self, serializer: S) -> Result { - (self.0.get() - 1).serialize(serializer) - } -} - -impl<'de> Deserialize<'de> for LocationIndex { - fn deserialize>(deserializer: D) -> Result { - let inner = u32::deserialize(deserializer)?; - - Ok(Self(unsafe { NonZeroU32::new_unchecked(inner + 1) })) + pub const fn index(&self) -> u32 { + self.index } } @@ -121,7 +79,7 @@ struct IndexedLocationRaw { x: u32, y: u32, #[serde(alias = "i")] - index: LocationIndex, + index: u32, } impl From for IndexedLocationRaw { diff --git a/necsim/core/src/lineage.rs b/necsim/core/src/lineage.rs index e7962e8ab..14320e63e 100644 --- a/necsim/core/src/lineage.rs +++ b/necsim/core/src/lineage.rs @@ -5,7 +5,7 @@ use core::{ use serde::{Deserialize, Deserializer, Serialize, Serializer}; -use necsim_core_bond::{NonNegativeF64, NonZeroOneU64, PositiveF64}; +use necsim_core_bond::{NonNegativeF64, PositiveF64}; use crate::{ cogs::{ @@ -17,31 +17,31 @@ use crate::{ #[derive(Clone, Debug, PartialEq, Eq, PartialOrd, Ord, Hash, TypeLayout)] #[repr(transparent)] -pub struct GlobalLineageReference(NonZeroOneU64); +pub struct GlobalLineageReference(u64); impl GlobalLineageReference { #[doc(hidden)] #[must_use] - pub unsafe fn into_inner(self) -> NonZeroOneU64 { + pub unsafe fn into_inner(self) -> u64 { self.0 } #[doc(hidden)] #[must_use] - pub unsafe fn from_inner(inner: NonZeroOneU64) -> Self { + pub unsafe fn from_inner(inner: u64) -> Self { Self(inner) } } impl fmt::Display for GlobalLineageReference { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - write!(f, "{}", self.0.get() - 2) + write!(f, "{}", self.0) } } impl Serialize for GlobalLineageReference { fn serialize(&self, serializer: S) -> Result { - (self.0.get() - 2).serialize(serializer) + self.0.serialize(serializer) } } @@ -49,7 +49,7 @@ impl<'de> Deserialize<'de> for GlobalLineageReference { fn deserialize>(deserializer: D) -> Result { let inner = u64::deserialize(deserializer)?; - Ok(Self(unsafe { NonZeroOneU64::new_unchecked(inner + 2) })) + Ok(Self(inner)) } } @@ -62,33 +62,25 @@ impl Backup for GlobalLineageReference { impl> LineageReference for GlobalLineageReference {} -#[allow(clippy::module_name_repetitions, clippy::unsafe_derive_deserialize)] -#[derive(Debug, Clone, Serialize, Deserialize, PartialOrd, Ord, TypeLayout)] -#[repr(transparent)] -pub struct LineageInteraction(u64); +#[allow(clippy::module_name_repetitions)] +#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord, Hash, Serialize, Deserialize)] +pub enum LineageInteraction { + None, + Maybe, + Coalescence(GlobalLineageReference), +} impl LineageInteraction { - #[allow(non_upper_case_globals)] - pub const Maybe: Self = Self(1_u64); - #[allow(non_upper_case_globals)] - pub const None: Self = Self(0_u64); - - #[allow(non_snake_case, clippy::needless_pass_by_value)] - #[must_use] - pub const fn Coalescence(parent: GlobalLineageReference) -> Self { - Self(parent.0.get()) - } - #[must_use] pub const fn is_coalescence(&self) -> bool { - self.0 > Self::Maybe.0 + matches!(self, Self::Coalescence(_)) } #[must_use] pub const fn parent(&self) -> Option { - match NonZeroOneU64::new(self.0) { - Ok(parent) => Some(GlobalLineageReference(parent)), - Err(_) => None, + match self { + Self::Coalescence(parent) => Some(GlobalLineageReference(parent.0)), + _ => None, } } } @@ -102,21 +94,6 @@ impl From> for LineageInteraction { } } -// Note: manually implementing PartialEq and Eq disables pattern matching -impl PartialEq for LineageInteraction { - fn eq(&self, other: &Self) -> bool { - self.0 == other.0 - } -} - -impl Eq for LineageInteraction {} - -impl core::hash::Hash for LineageInteraction { - fn hash(&self, state: &mut H) { - self.0.hash(state); - } -} - #[allow(clippy::unsafe_derive_deserialize)] #[derive(Debug, Clone, PartialEq, Eq, Hash, Serialize, Deserialize, TypeLayout)] #[serde(deny_unknown_fields)] @@ -130,12 +107,6 @@ pub struct Lineage { pub indexed_location: IndexedLocation, } -#[allow(dead_code)] -const EXCESSIVE_OPTION_LINEAGE_ERROR: [(); 1 - { - const ASSERT: bool = core::mem::size_of::>() == core::mem::size_of::(); - ASSERT -} as usize] = []; - impl Lineage { #[must_use] #[debug_ensures( @@ -148,11 +119,9 @@ impl Lineage { habitat: &H, ) -> Self { Self { - global_reference: GlobalLineageReference(unsafe { - NonZeroOneU64::new_unchecked( - habitat.map_indexed_location_to_u64_injective(&indexed_location) + 2, - ) - }), + global_reference: GlobalLineageReference( + habitat.map_indexed_location_to_u64_injective(&indexed_location), + ), last_event_time: NonNegativeF64::zero(), indexed_location, } diff --git a/necsim/core/src/reporter/impl.rs b/necsim/core/src/reporter/impl.rs index b95062126..315086200 100644 --- a/necsim/core/src/reporter/impl.rs +++ b/necsim/core/src/reporter/impl.rs @@ -2,67 +2,67 @@ #[allow(clippy::module_name_repetitions)] macro_rules! impl_report { // Special case: Ignored = MaybeUsed - ($(#[$metas:meta])* $target:ident(&mut $this:ident, $value:ident: Ignored) {}) => + ($(#[$metas:meta])* $([$default:tt])? $target:ident(&mut $this:ident, $value:ident: Ignored) {}) => { impl_report!{ $(#[$metas])* - $target(&mut $this, $value: MaybeUsed<$crate::reporter::boolean::False>) {} + $([$default])? $target(&mut $this, $value: MaybeUsed<$crate::reporter::boolean::False>) {} } }; // Special case: Used = MaybeUsed - ($(#[$metas:meta])* $target:ident(&mut $this:ident, $value:ident: Used) $code:block) => + ($(#[$metas:meta])* $([$default:tt])? $target:ident(&mut $this:ident, $value:ident: Used) $code:block) => { impl_report!{ $(#[$metas])* - $target(&mut $this, $value: MaybeUsed<$crate::reporter::boolean::True>) + $([$default])? $target(&mut $this, $value: MaybeUsed<$crate::reporter::boolean::True>) $code } }; // Dispatch case: MaybeUsed + speciation - ($(#[$metas:meta])* speciation(&mut $this:ident, $value:ident: MaybeUsed<$Usage:ty>) + ($(#[$metas:meta])* $([$default:tt])? speciation(&mut $this:ident, $value:ident: MaybeUsed<$Usage:ty>) $code:block) => { impl_report!{ $(#[$metas])* - fn report_speciation(&mut $this, $value: MaybeUsed< + $([$default])? fn report_speciation(&mut $this, $value: MaybeUsed< $crate::event::SpeciationEvent, ReportSpeciation = $Usage >) $code } }; // Dispatch case: MaybeUsed + dispersal - ($(#[$metas:meta])* dispersal(&mut $this:ident, $value:ident: MaybeUsed<$Usage:ty>) + ($(#[$metas:meta])* $([$default:tt])? dispersal(&mut $this:ident, $value:ident: MaybeUsed<$Usage:ty>) $code:block) => { impl_report!{ $(#[$metas])* - fn report_dispersal(&mut $this, $value: MaybeUsed< + $([$default])? fn report_dispersal(&mut $this, $value: MaybeUsed< $crate::event::DispersalEvent, ReportDispersal = $Usage >) $code } }; // Dispatch case: MaybeUsed + progress - ($(#[$metas:meta])* progress(&mut $this:ident, $value:ident: MaybeUsed<$Usage:ty>) + ($(#[$metas:meta])* $([$default:tt])? progress(&mut $this:ident, $value:ident: MaybeUsed<$Usage:ty>) $code:block) => { impl_report!{ $(#[$metas])* - fn report_progress(&mut $this, $value: MaybeUsed) $code + $([$default])? fn report_progress(&mut $this, $value: MaybeUsed) $code } }; // Impl case: MaybeUsed - ($(#[$metas:meta])* fn $target:ident(&mut $this:ident, $value:ident: MaybeUsed< + ($(#[$metas:meta])* $([$default:tt])? fn $target:ident(&mut $this:ident, $value:ident: MaybeUsed< $EventTy:ty, $UsageIdent:ident = $UsageTy:ty >) $code:block) => { $(#[$metas])* - fn $target( + $($default)? fn $target( &mut $this, $value: &$crate::reporter::used::MaybeUsed<$EventTy, Self::$UsageIdent>, ) { $value.maybe_use_in(|$value| $code) } - type $UsageIdent = $UsageTy; + $($default)? type $UsageIdent = $UsageTy; }; } diff --git a/necsim/impls/cuda/src/event_buffer.rs b/necsim/impls/cuda/src/event_buffer.rs index 394bde17f..b90c9d260 100644 --- a/necsim/impls/cuda/src/event_buffer.rs +++ b/necsim/impls/cuda/src/event_buffer.rs @@ -1,4 +1,4 @@ -use core::{fmt, marker::PhantomData}; +use core::fmt; #[cfg(not(target_os = "cuda"))] use rust_cuda::rustacuda::{ @@ -7,13 +7,15 @@ use rust_cuda::rustacuda::{ }; use rust_cuda::utils::{ - aliasing::{SplitSliceOverCudaThreadsConstStride, SplitSliceOverCudaThreadsDynamicStride}, - exchange::buffer::CudaExchangeBuffer, + aliasing::SplitSliceOverCudaThreadsDynamicStride, exchange::buffer::CudaExchangeBuffer, }; use necsim_core::{ - event::{DispersalEvent, SpeciationEvent}, - reporter::{boolean::Boolean, Reporter}, + event::{PackedEvent, SpeciationEvent, TypedEvent}, + reporter::{ + boolean::{Boolean, False, True}, + Reporter, + }, }; #[cfg(target_os = "cuda")] @@ -21,34 +23,60 @@ use necsim_core::impl_report; use super::utils::MaybeSome; -#[allow(clippy::module_name_repetitions)] +#[allow(clippy::module_name_repetitions, clippy::type_complexity)] #[derive(rust_cuda::common::LendRustToCuda)] #[cuda(free = "ReportSpeciation", free = "ReportDispersal")] pub struct EventBuffer { #[cuda(embed)] - speciation_mask: - SplitSliceOverCudaThreadsConstStride, 1_usize>, + event_mask: SplitSliceOverCudaThreadsDynamicStride>, #[cuda(embed)] - speciation_buffer: SplitSliceOverCudaThreadsConstStride< - CudaExchangeBuffer, false, true>, - 1_usize, - >, - #[cuda(embed)] - dispersal_mask: SplitSliceOverCudaThreadsDynamicStride>, - #[cuda(embed)] - dispersal_buffer: SplitSliceOverCudaThreadsDynamicStride< - CudaExchangeBuffer, false, true>, + event_buffer: SplitSliceOverCudaThreadsDynamicStride< + CudaExchangeBuffer< + MaybeSome< as EventType>::Event>, + false, + true, + >, >, max_events: usize, event_counter: usize, - marker: PhantomData<(ReportSpeciation, ReportDispersal)>, +} + +pub trait EventType { + type Event: 'static + + ~const rust_cuda::const_type_layout::TypeGraphLayout + + rust_cuda::safety::StackOnly + + Into + + Into + + Clone; +} + +impl EventType + for EventBuffer +{ + default type Event = PackedEvent; +} + +impl EventType for EventBuffer { + type Event = PackedEvent; +} + +impl EventType for EventBuffer { + type Event = PackedEvent; +} + +impl EventType for EventBuffer { + type Event = SpeciationEvent; +} + +impl EventType for EventBuffer { + type Event = PackedEvent; } impl fmt::Debug for EventBuffer { fn fmt(&self, fmt: &mut fmt::Formatter) -> fmt::Result { - fmt.debug_struct(stringify!(EventBuffer)) + fmt.debug_struct("EventBuffer") .field("max_events", &self.max_events) .field("event_counter", &self.event_counter) .finish() @@ -69,72 +97,75 @@ impl let block_size = (block_size.x * block_size.y * block_size.z) as usize; let grid_size = (grid_size.x * grid_size.y * grid_size.z) as usize; - let speciation_capacity = if ReportSpeciation::VALUE { - block_size * grid_size - } else { - 0_usize - }; - let dispersal_capacity = if ReportDispersal::VALUE { - max_events * block_size * grid_size + #[allow(clippy::bool_to_int_with_if)] + let max_events = if ReportDispersal::VALUE { + max_events + } else if ReportSpeciation::VALUE { + 1_usize } else { 0_usize }; - let mut speciation_buffer = alloc::vec::Vec::with_capacity(speciation_capacity); - speciation_buffer.resize_with(speciation_capacity, || MaybeSome::None); + let event_capacity = max_events * block_size * grid_size; - let mut dispersal_buffer = alloc::vec::Vec::with_capacity(dispersal_capacity); - dispersal_buffer.resize_with(dispersal_capacity, || MaybeSome::None); + let mut event_buffer = alloc::vec::Vec::with_capacity(event_capacity); + event_buffer.resize_with(event_capacity, || MaybeSome::None); Ok(Self { - speciation_mask: SplitSliceOverCudaThreadsConstStride::new(CudaExchangeBuffer::new( - &false, - speciation_capacity, - )?), - speciation_buffer: SplitSliceOverCudaThreadsConstStride::new( - CudaExchangeBuffer::from_vec(speciation_buffer)?, - ), - dispersal_mask: SplitSliceOverCudaThreadsDynamicStride::new( - CudaExchangeBuffer::new(&false, dispersal_capacity)?, + event_mask: SplitSliceOverCudaThreadsDynamicStride::new( + CudaExchangeBuffer::new(&false, event_capacity)?, max_events, ), - dispersal_buffer: SplitSliceOverCudaThreadsDynamicStride::new( - CudaExchangeBuffer::from_vec(dispersal_buffer)?, + event_buffer: SplitSliceOverCudaThreadsDynamicStride::new( + CudaExchangeBuffer::from_vec(event_buffer)?, max_events, ), max_events, event_counter: 0_usize, - marker: PhantomData::<(ReportSpeciation, ReportDispersal)>, }) } - pub fn report_events

(&mut self, reporter: &mut P) + #[allow(clippy::missing_panics_doc)] // TODO: remove + pub fn report_events_unordered

(&mut self, reporter: &mut P) where P: Reporter, { - for (mask, dispersal) in self - .dispersal_mask - .iter_mut() - .zip(self.dispersal_buffer.iter()) - { - if ReportDispersal::VALUE && *mask.read() { - reporter.report_dispersal(unsafe { dispersal.read().assume_some_ref() }.into()); - } + // let mut last_time = 0.0_f64; - mask.write(false); - } + // let mut times = alloc::vec::Vec::new(); - for (mask, speciation) in self - .speciation_mask - .iter_mut() - .zip(self.speciation_buffer.iter()) - { - if ReportSpeciation::VALUE && *mask.read() { - reporter.report_speciation(unsafe { speciation.read().assume_some_ref() }.into()); - } + for (mask, event) in self.event_mask.iter_mut().zip(self.event_buffer.iter()) { + if *mask.read() { + let event: TypedEvent = unsafe { event.read().assume_some_read() }.into(); + // let new_time: f64 = match &event { + // TypedEvent::Speciation(speciation) => speciation.event_time, + // TypedEvent::Dispersal(dispersal) => dispersal.event_time, + // } + // .get(); + // times.push(Some(new_time)); + // assert!(new_time >= last_time, "{new_time} {last_time}"); + // last_time = new_time; + + match event { + TypedEvent::Speciation(ref speciation) => { + reporter.report_speciation(speciation.into()); + }, + TypedEvent::Dispersal(ref dispersal) => { + reporter.report_dispersal(dispersal.into()); + }, + } + } /*else { + times.push(None); + }*/ mask.write(false); } + + // panic!("{:?}", times); + } + + pub fn max_events_per_individual(&self) -> usize { + self.max_events } } @@ -142,23 +173,72 @@ impl impl Reporter for EventBuffer { + impl_report!([default] speciation(&mut self, _event: Ignored) {}); + + impl_report!([default] dispersal(&mut self, _event: Ignored) {}); + + impl_report!([default] progress(&mut self, _progress: Ignored) {}); +} + +#[cfg(target_os = "cuda")] +impl Reporter for EventBuffer { impl_report!( #[debug_requires( - !self.speciation_mask.get(0).map( - rust_cuda::utils::exchange::buffer::CudaExchangeItem::read - ).copied().unwrap_or(true), - "does not report extraneous speciation event" + self.event_counter < self.max_events, + "does not report extraneous dispersal events" + )] + dispersal(&mut self, event: Used) { + if let Some(mask) = self.event_mask.get_mut(self.event_counter) { + mask.write(true); + + unsafe { + self.event_buffer.get_unchecked_mut(self.event_counter) + }.write(MaybeSome::Some(event.clone().into())); + } + + self.event_counter += 1; + } + ); +} + +#[cfg(target_os = "cuda")] +impl Reporter for EventBuffer { + impl_report!( + #[debug_requires( + self.event_counter == 0, + "does not report extraneous speciation events" )] speciation(&mut self, event: Used) { - if ReportSpeciation::VALUE { - if let Some(mask) = self.speciation_mask.get_mut(0) { - mask.write(true); + if let Some(mask) = self.event_mask.get_mut(0) { + mask.write(true); - unsafe { - self.speciation_buffer.get_unchecked_mut(0) - }.write(MaybeSome::Some(event.clone())); - } + unsafe { + self.event_buffer.get_unchecked_mut(0) + }.write(MaybeSome::Some(event.clone())); + } + + self.event_counter = self.max_events; + } + ); +} + +#[cfg(target_os = "cuda")] +impl Reporter for EventBuffer { + impl_report!( + #[debug_requires( + self.event_counter < self.max_events, + "does not report extraneous speciation events" + )] + speciation(&mut self, event: Used) { + if let Some(mask) = self.event_mask.get_mut(self.event_counter) { + mask.write(true); + + unsafe { + self.event_buffer.get_unchecked_mut(self.event_counter) + }.write(MaybeSome::Some(event.clone().into())); } + + self.event_counter = self.max_events; } ); @@ -168,19 +248,15 @@ impl Reporter "does not report extraneous dispersal events" )] dispersal(&mut self, event: Used) { - if ReportDispersal::VALUE { - if let Some(mask) = self.dispersal_mask.get_mut(self.event_counter) { - mask.write(true); + if let Some(mask) = self.event_mask.get_mut(self.event_counter) { + mask.write(true); - unsafe { - self.dispersal_buffer.get_unchecked_mut(self.event_counter) - }.write(MaybeSome::Some(event.clone())); - } - - self.event_counter += 1; + unsafe { + self.event_buffer.get_unchecked_mut(self.event_counter) + }.write(MaybeSome::Some(event.clone().into())); } + + self.event_counter += 1; } ); - - impl_report!(progress(&mut self, _progress: Ignored) {}); } diff --git a/necsim/impls/cuda/src/lib.rs b/necsim/impls/cuda/src/lib.rs index 8e9019c11..44c4c984d 100644 --- a/necsim/impls/cuda/src/lib.rs +++ b/necsim/impls/cuda/src/lib.rs @@ -8,6 +8,8 @@ #![cfg_attr(target_os = "cuda", feature(asm_experimental_arch))] #![cfg_attr(target_os = "cuda", feature(asm_const))] #![cfg_attr(target_os = "cuda", feature(const_float_bits_conv))] +#![allow(incomplete_features)] +#![feature(specialization)] extern crate alloc; diff --git a/necsim/impls/no-std/src/cogs/active_lineage_sampler/alias/individual/mod.rs b/necsim/impls/no-std/src/cogs/active_lineage_sampler/alias/individual/mod.rs index ec7c90257..94049419d 100644 --- a/necsim/impls/no-std/src/cogs/active_lineage_sampler/alias/individual/mod.rs +++ b/necsim/impls/no-std/src/cogs/active_lineage_sampler/alias/individual/mod.rs @@ -90,7 +90,7 @@ impl< while let Some(lineage) = origin_sampler.next() { if !origin_sampler .habitat() - .contains(lineage.indexed_location.location()) + .is_location_habitable(lineage.indexed_location.location()) { exceptional_lineages.push(ExceptionalLineage::OutOfHabitat(lineage)); continue; diff --git a/necsim/impls/no-std/src/cogs/active_lineage_sampler/alias/location/mod.rs b/necsim/impls/no-std/src/cogs/active_lineage_sampler/alias/location/mod.rs index 18baad763..b8c4e29fc 100644 --- a/necsim/impls/no-std/src/cogs/active_lineage_sampler/alias/location/mod.rs +++ b/necsim/impls/no-std/src/cogs/active_lineage_sampler/alias/location/mod.rs @@ -109,7 +109,7 @@ impl< while let Some(lineage) = origin_sampler.next() { if !origin_sampler .habitat() - .contains(lineage.indexed_location.location()) + .is_location_habitable(lineage.indexed_location.location()) { exceptional_lineages.push(ExceptionalLineage::OutOfHabitat(lineage)); continue; diff --git a/necsim/impls/no-std/src/cogs/active_lineage_sampler/classical/mod.rs b/necsim/impls/no-std/src/cogs/active_lineage_sampler/classical/mod.rs index bb9cda39c..2bd683d45 100644 --- a/necsim/impls/no-std/src/cogs/active_lineage_sampler/classical/mod.rs +++ b/necsim/impls/no-std/src/cogs/active_lineage_sampler/classical/mod.rs @@ -77,7 +77,7 @@ impl< while let Some(lineage) = origin_sampler.next() { if !origin_sampler .habitat() - .contains(lineage.indexed_location.location()) + .is_location_habitable(lineage.indexed_location.location()) { exceptional_lineages.push(ExceptionalLineage::OutOfHabitat(lineage)); continue; diff --git a/necsim/impls/no-std/src/cogs/active_lineage_sampler/independent/mod.rs b/necsim/impls/no-std/src/cogs/active_lineage_sampler/independent/mod.rs index 9bd4927cf..1aafbee33 100644 --- a/necsim/impls/no-std/src/cogs/active_lineage_sampler/independent/mod.rs +++ b/necsim/impls/no-std/src/cogs/active_lineage_sampler/independent/mod.rs @@ -101,7 +101,7 @@ impl< while let Some(lineage) = origin_sampler.next() { if !origin_sampler .habitat() - .contains(lineage.indexed_location.location()) + .is_location_habitable(lineage.indexed_location.location()) { exceptional_lineages.push(ExceptionalLineage::OutOfHabitat(lineage)); continue; diff --git a/necsim/impls/no-std/src/cogs/coalescence_sampler/independent.rs b/necsim/impls/no-std/src/cogs/coalescence_sampler/independent.rs index e4c27f72f..0e9a16f6a 100644 --- a/necsim/impls/no-std/src/cogs/coalescence_sampler/independent.rs +++ b/necsim/impls/no-std/src/cogs/coalescence_sampler/independent.rs @@ -8,7 +8,10 @@ use necsim_core::{ lineage::LineageInteraction, }; -use crate::cogs::lineage_store::independent::IndependentLineageStore; +use crate::cogs::lineage_store::{ + coherent::globally::singleton_demes::SingletonDemesHabitat, + independent::IndependentLineageStore, +}; #[allow(clippy::module_name_repetitions)] #[derive(Debug)] @@ -34,19 +37,44 @@ impl> CoalescenceSampler { #[must_use] - #[debug_ensures(ret.1 == LineageInteraction::Maybe, "always reports maybe")] - fn sample_interaction_at_location( + // #[debug_ensures(ret.1 == LineageInteraction::Maybe, "always reports maybe")] + default fn sample_interaction_at_location( &self, location: Location, habitat: &H, _lineage_store: &IndependentLineageStore, coalescence_rng_sample: CoalescenceRngSample, ) -> (IndexedLocation, LineageInteraction) { - let chosen_coalescence_index = coalescence_rng_sample - .sample_coalescence_index::(habitat.get_habitat_at_location(&location)); + let population = habitat.get_habitat_at_location(&location); + + let chosen_coalescence_index = + coalescence_rng_sample.sample_coalescence_index::(population); let indexed_location = IndexedLocation::new(location, chosen_coalescence_index); (indexed_location, LineageInteraction::Maybe) } } + +// Specialise for SingletonDemesHabitat as the compiler cannot yet optimise out +// the call to `habitat.get_habitat_at_location(&location)`. +#[contract_trait] +impl> + CoalescenceSampler> + for IndependentCoalescenceSampler +{ + #[must_use] + #[debug_ensures(ret.1 == LineageInteraction::Maybe, "always reports maybe")] + fn sample_interaction_at_location( + &self, + location: Location, + _habitat: &H, + _lineage_store: &IndependentLineageStore, + _coalescence_rng_sample: CoalescenceRngSample, + ) -> (IndexedLocation, LineageInteraction) { + ( + IndexedLocation::new(location, 0_u32), + LineageInteraction::Maybe, + ) + } +} diff --git a/necsim/impls/no-std/src/cogs/dispersal_sampler/spatially_implicit.rs b/necsim/impls/no-std/src/cogs/dispersal_sampler/spatially_implicit.rs index 2b0cc3ef4..9664e50bb 100644 --- a/necsim/impls/no-std/src/cogs/dispersal_sampler/spatially_implicit.rs +++ b/necsim/impls/no-std/src/cogs/dispersal_sampler/spatially_implicit.rs @@ -49,13 +49,11 @@ impl> DispersalSampler { #[must_use] - #[debug_requires( - habitat.local().contains(location) || habitat.meta().contains(location), - "location is inside either the local or meta habitat extent" - )] - #[debug_ensures(habitat.meta().contains(&ret) || if old(habitat.local().contains(location)) { - habitat.local().contains(&ret) - } else { false }, "target is inside the meta habitat extent, \ + #[debug_ensures(habitat.meta().get_extent().contains(&ret) || ( + if old(habitat.local().get_extent().contains(location)) { + habitat.local().get_extent().contains(&ret) + } else { false } + ), "target is inside the meta habitat extent, \ or -- iff the location was local -- in the local habitat extent" )] fn sample_dispersal_from_location( @@ -66,7 +64,9 @@ impl> DispersalSampler Location { use necsim_core::cogs::RngSampler; - if habitat.local().contains(location) { + // By PRE, location must be habitable, i.e. either in the local or the meta + // habitat + if habitat.local().get_extent().contains(location) { if rng.sample_event(self.local_migration_probability_per_generation.into()) { // Provide a dummpy Location in the meta community to disperse from self.meta.sample_dispersal_from_location( @@ -93,13 +93,11 @@ impl> SeparableDispersalSampler { #[must_use] - #[debug_requires( - habitat.local().contains(location) || habitat.meta().contains(location), - "location is inside either the local or meta habitat extent" - )] - #[debug_ensures(habitat.meta().contains(&ret) || if old(habitat.local().contains(location)) { - habitat.local().contains(&ret) - } else { false }, "target is inside the meta habitat extent, \ + #[debug_ensures(habitat.meta().get_extent().contains(&ret) || ( + if old(habitat.local().get_extent().contains(location)) { + habitat.local().get_extent().contains(&ret) + } else { false } + ), "target is inside the meta habitat extent, \ or -- iff the location was local -- in the local habitat extent" )] fn sample_non_self_dispersal_from_location( @@ -110,7 +108,9 @@ impl> SeparableDispersalSampler Location { use necsim_core::cogs::RngSampler; - if habitat.local().contains(location) { + // By PRE, location must be habitable, i.e. either in the local or the meta + // habitat + if habitat.local().get_extent().contains(location) { if rng.sample_event(self.local_migration_probability_per_generation.into()) { // Provide a dummpy Location in the meta community to disperse from // As the individual is dispersing to a different community, @@ -134,16 +134,14 @@ impl> SeparableDispersalSampler, ) -> ClosedUnitF64 { - if habitat.local().contains(location) { + // By PRE, location must be habitable, i.e. either in the local or the meta + // habitat + if habitat.local().get_extent().contains(location) { self.local .get_self_dispersal_probability_at_location(location, habitat.local()) * ClosedUnitF64::from(self.local_migration_probability_per_generation).one_minus() diff --git a/necsim/impls/no-std/src/cogs/dispersal_sampler/trespassing/mod.rs b/necsim/impls/no-std/src/cogs/dispersal_sampler/trespassing/mod.rs index b98b75fd7..996dc2684 100644 --- a/necsim/impls/no-std/src/cogs/dispersal_sampler/trespassing/mod.rs +++ b/necsim/impls/no-std/src/cogs/dispersal_sampler/trespassing/mod.rs @@ -16,8 +16,8 @@ pub trait AntiTrespassingDispersalSampler, G: RngCor Backup + core::fmt::Debug { #[must_use] - #[debug_requires(!habitat.contains(location), "location is outside habitat")] - #[debug_ensures(old(habitat).contains(&ret), "target is inside habitat")] + #[debug_requires(!habitat.is_location_habitable(location), "location is inhabitable")] + #[debug_ensures(old(habitat).is_location_habitable(&ret), "target is habitable")] fn sample_anti_trespassing_dispersal_from_location( &self, location: &Location, @@ -91,14 +91,15 @@ impl< #[must_use] #[inline] #[allow(clippy::no_effect_underscore_binding)] - #[debug_ensures(old(habitat).contains(&ret), "target is inside habitat")] + #[debug_ensures(old(habitat).is_location_habitable(&ret), "target is habitable")] fn sample_dispersal_from_location( &self, location: &Location, habitat: &H, rng: &mut G, ) -> Location { - // Explicitly circumvent the precondition that `habitat.contains(location)` + // Explicitly circumvent the precondition that + // `habitat.is_location_habitable(location)` self.__contracts_impl_sample_dispersal_from_location(location, habitat, rng) } @@ -110,12 +111,12 @@ impl< habitat: &H, rng: &mut G, ) -> Location { - if habitat.contains(location) { - // Re-establish the precondition that `habitat.contains(location)` + if habitat.is_location_habitable(location) { + // Re-establish the precondition that `habitat.is_location_habitable(location)` self.legal_dispersal_sampler .sample_dispersal_from_location(location, habitat, rng) } else { - // Establish the precondition that `!habitat.contains(location)` + // Establish the precondition that `!habitat.is_location_habitable(location)` self.trespassing_dispersal_sampler .sample_anti_trespassing_dispersal_from_location(location, habitat, rng) } @@ -133,7 +134,7 @@ impl< #[must_use] #[inline] #[allow(clippy::no_effect_underscore_binding)] - #[debug_ensures(old(habitat).contains(&ret), "target is inside habitat")] + #[debug_ensures(old(habitat).is_location_habitable(&ret), "target is habitable")] #[debug_ensures(&ret != location, "disperses to a different location")] fn sample_non_self_dispersal_from_location( &self, @@ -141,7 +142,8 @@ impl< habitat: &H, rng: &mut G, ) -> Location { - // Explicitly circumvent the precondition that `habitat.contains(location)` + // Explicitly circumvent the precondition that + // `habitat.is_location_habitable(location)` self.__contracts_impl_sample_non_self_dispersal_from_location(location, habitat, rng) } @@ -153,12 +155,15 @@ impl< habitat: &H, rng: &mut G, ) -> Location { - if habitat.contains(location) { - // Re-establish the precondition that `habitat.contains(location)` + if habitat.is_location_habitable(location) { + // Re-establish the precondition that `habitat.is_location_habitable(location)` self.legal_dispersal_sampler .sample_non_self_dispersal_from_location(location, habitat, rng) } else { - // Establish the precondition that `!habitat.contains(location)` + // Establish the precondition that `!habitat.is_location_habitable(location)` + // Since the origin location is inhabitable but the target location is + // habitable, the target location must be different from the origin + // location self.trespassing_dispersal_sampler .sample_anti_trespassing_dispersal_from_location(location, habitat, rng) } @@ -171,7 +176,8 @@ impl< location: &Location, habitat: &H, ) -> ClosedUnitF64 { - // Explicitly circumvent the precondition that `habitat.contains(location)` + // Explicitly circumvent the precondition that + // `habitat.is_location_habitable(location)` self.__contracts_impl_get_self_dispersal_probability_at_location(location, habitat) } @@ -182,14 +188,14 @@ impl< location: &Location, habitat: &H, ) -> ClosedUnitF64 { - if habitat.contains(location) { - // Re-establish the precondition that `habitat.contains(location)` + if habitat.is_location_habitable(location) { + // Re-establish the precondition that `habitat.is_location_habitable(location)` self.legal_dispersal_sampler .get_self_dispersal_probability_at_location(location, habitat) } else { - // The `AntiTrespassingDispersalSampler` always jumps from outside - // the habitat to inside the habitat, i.e. there is never any - // self-dispersal outside the habitat + // The `AntiTrespassingDispersalSampler` always jumps from an inhabitable + // location to a habitable location, i.e. there is never any self-dispersal + // when starting at an inhabitable location ClosedUnitF64::zero() } } diff --git a/necsim/impls/no-std/src/cogs/dispersal_sampler/wrapping_noise.rs b/necsim/impls/no-std/src/cogs/dispersal_sampler/wrapping_noise.rs index d85ad4114..5f38306db 100644 --- a/necsim/impls/no-std/src/cogs/dispersal_sampler/wrapping_noise.rs +++ b/necsim/impls/no-std/src/cogs/dispersal_sampler/wrapping_noise.rs @@ -95,22 +95,20 @@ impl> SeparableDispersalSampler, ) -> ClosedUnitF64 { - if habitat.get_habitat_at_location(location) == 0 { - ClosedUnitF64::zero() - } else { - let p_self_dispersal = self - .inner - .get_self_dispersal_probability_at_location(location, habitat.get_inner()); - let p_out_dispersal = p_self_dispersal.one_minus() * habitat.coverage(); + // By PRE, the location is habitable, i.e. self-dispersal is possible + + let p_self_dispersal = self + .inner + .get_self_dispersal_probability_at_location(location, habitat.get_inner()); + let p_out_dispersal = p_self_dispersal.one_minus() * habitat.coverage(); - // Safety: - // - p_self_dispersal and p_out_dispersal are both in [0, 1] - // - a / (a + [0, 1]) = [0, 1] - unsafe { - ClosedUnitF64::new_unchecked( - p_self_dispersal.get() / (p_self_dispersal.get() + p_out_dispersal.get()), - ) - } + // Safety: + // - p_self_dispersal and p_out_dispersal are both in [0, 1] + // - a / (a + [0, 1]) = [0, 1] + unsafe { + ClosedUnitF64::new_unchecked( + p_self_dispersal.get() / (p_self_dispersal.get() + p_out_dispersal.get()), + ) } } } diff --git a/necsim/impls/no-std/src/cogs/event_sampler/tracking.rs b/necsim/impls/no-std/src/cogs/event_sampler/tracking.rs index 456a427ed..8b5c1cccd 100644 --- a/necsim/impls/no-std/src/cogs/event_sampler/tracking.rs +++ b/necsim/impls/no-std/src/cogs/event_sampler/tracking.rs @@ -34,13 +34,6 @@ pub struct SpeciationSample { sample_location: IndexedLocation, } -#[allow(dead_code)] -const EXCESSIVE_OPTION_SPECIATION_SAMPLE_ERROR: [(); 1 - { - const ASSERT: bool = core::mem::size_of::>() - == core::mem::size_of::(); - ASSERT -} as usize] = []; - impl SpeciationSample { pub fn update_min( min_spec_sample: &mut Option, diff --git a/necsim/impls/no-std/src/cogs/habitat/almost_infinite.rs b/necsim/impls/no-std/src/cogs/habitat/almost_infinite.rs index 91cfc423c..914672dbc 100644 --- a/necsim/impls/no-std/src/cogs/habitat/almost_infinite.rs +++ b/necsim/impls/no-std/src/cogs/habitat/almost_infinite.rs @@ -1,4 +1,4 @@ -use core::marker::PhantomData; +use core::{fmt, marker::PhantomData}; use necsim_core::{ cogs::{Backup, Habitat, MathsCore, RngCore, UniformlySampleableHabitat}, @@ -8,19 +8,25 @@ use necsim_core_bond::{OffByOneU32, OffByOneU64}; use crate::cogs::lineage_store::coherent::globally::singleton_demes::SingletonDemesHabitat; +const ALMOST_INFINITE_EXTENT: LandscapeExtent = + LandscapeExtent::new(0, 0, OffByOneU32::max(), OffByOneU32::max()); + #[allow(clippy::module_name_repetitions)] -#[derive(Debug)] #[cfg_attr(feature = "cuda", derive(rust_cuda::common::LendRustToCuda))] #[cfg_attr(feature = "cuda", cuda(free = "M"))] pub struct AlmostInfiniteHabitat { - extent: LandscapeExtent, marker: PhantomData, } +impl fmt::Debug for AlmostInfiniteHabitat { + fn fmt(&self, fmt: &mut fmt::Formatter) -> fmt::Result { + fmt.debug_struct(stringify!(AlmostInfiniteHabitat)).finish() + } +} + impl Default for AlmostInfiniteHabitat { fn default() -> Self { Self { - extent: LandscapeExtent::new(0_u32, 0_u32, OffByOneU32::max(), OffByOneU32::max()), marker: PhantomData::, } } @@ -30,7 +36,6 @@ impl Default for AlmostInfiniteHabitat { impl Backup for AlmostInfiniteHabitat { unsafe fn backup_unchecked(&self) -> Self { Self { - extent: self.extent.clone(), marker: PhantomData::, } } @@ -47,7 +52,7 @@ impl Habitat for AlmostInfiniteHabitat { #[must_use] fn get_extent(&self) -> &LandscapeExtent { - &self.extent + &ALMOST_INFINITE_EXTENT } #[must_use] diff --git a/necsim/impls/no-std/src/cogs/habitat/spatially_implicit.rs b/necsim/impls/no-std/src/cogs/habitat/spatially_implicit.rs index 3736927b5..5f78012e9 100644 --- a/necsim/impls/no-std/src/cogs/habitat/spatially_implicit.rs +++ b/necsim/impls/no-std/src/cogs/habitat/spatially_implicit.rs @@ -8,6 +8,9 @@ use necsim_core_bond::{OffByOneU32, OffByOneU64}; use crate::cogs::habitat::non_spatial::NonSpatialHabitat; +const SPATIALLY_IMPLICIT_EXTENT: LandscapeExtent = + LandscapeExtent::new(0, 0, OffByOneU32::max(), OffByOneU32::max()); + #[allow(clippy::module_name_repetitions)] #[derive(Debug)] #[cfg_attr(feature = "cuda", derive(rust_cuda::common::LendRustToCuda))] @@ -17,7 +20,6 @@ pub struct SpatiallyImplicitHabitat { local: NonSpatialHabitat, #[cfg_attr(feature = "cuda", cuda(embed))] meta: NonSpatialHabitat, - extent: LandscapeExtent, } impl SpatiallyImplicitHabitat { @@ -47,11 +49,7 @@ impl SpatiallyImplicitHabitat { meta_deme, ); - Self { - extent: LandscapeExtent::new(0, 0, OffByOneU32::max(), OffByOneU32::max()), - local, - meta, - } + Self { local, meta } } #[must_use] @@ -71,7 +69,6 @@ impl Backup for SpatiallyImplicitHabitat { Self { local: self.local.backup_unchecked(), meta: self.meta.backup_unchecked(), - extent: self.extent.clone(), } } } @@ -87,7 +84,7 @@ impl Habitat for SpatiallyImplicitHabitat { #[must_use] fn get_extent(&self) -> &LandscapeExtent { - &self.extent + &SPATIALLY_IMPLICIT_EXTENT } #[must_use] diff --git a/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/mod.rs b/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/mod.rs index 86627bfdb..e86345216 100644 --- a/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/mod.rs +++ b/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/mod.rs @@ -1,5 +1,5 @@ use alloc::boxed::Box; -use core::{f64::consts::PI, fmt, num::NonZeroUsize}; +use core::{fmt, num::NonZeroUsize}; use necsim_core_bond::{ClosedUnitF64, OffByOneU64, OpenClosedUnitF64 as PositiveUnitF64}; use r#final::Final; @@ -14,7 +14,7 @@ use necsim_core::{ use crate::cogs::{ habitat::almost_infinite::AlmostInfiniteHabitat, - lineage_store::coherent::globally::singleton_demes::SingletonDemesHabitat, + lineage_store::coherent::globally::singleton_demes::SingletonDemesHabitat, rng::wyhash::WyHash, }; #[allow(clippy::module_name_repetitions)] @@ -35,7 +35,6 @@ pub struct WrappingNoiseHabitat { impl fmt::Debug for WrappingNoiseHabitat { fn fmt(&self, fmt: &mut fmt::Formatter) -> fmt::Result { fmt.debug_struct(stringify!(WrappingNoiseHabitat)) - .field("inner", &self.inner) .field("coverage", &self.coverage) .field("scale", &self.scale) .field("persistence", &self.persistence) @@ -53,27 +52,29 @@ impl WrappingNoiseHabitat { persistence: PositiveUnitF64, octaves: NonZeroUsize, ) -> Self { - const THRESHOLD_QUALITY_LOG2: u32 = 8; - let noise = Box::new(OpenSimplexNoise::new(Some(seed))); // Emperically determine a threshold to uniformly sample habitat // from the generated Simplex Noise let mut samples = alloc::vec::Vec::new(); - for y in 0..(1 << THRESHOLD_QUALITY_LOG2) { - for x in 0..(1 << THRESHOLD_QUALITY_LOG2) { - samples.push(sum_noise_octaves::( - &noise, - &Location::new( - x << (32 - THRESHOLD_QUALITY_LOG2), - y << (32 - THRESHOLD_QUALITY_LOG2), - ), - persistence, - scale, - octaves, - )); - } + // Utilise a PRNG to avoid sampling degeneracies for finding the + // threshold which would poison the entire sampler + let mut rng: WyHash = WyHash::from_seed(seed.to_le_bytes()); + + for _ in 0..(1_usize << 16) { + let location = rng.sample_u64(); + + samples.push(sum_noise_octaves::( + &noise, + &Location::new( + (location & 0x0000_0000_FFFF_FFFF) as u32, + ((location >> 32) & 0x0000_0000_FFFF_FFFF) as u32, + ), + persistence, + scale, + octaves, + )); } samples.sort_unstable_by(f64::total_cmp); @@ -211,8 +212,6 @@ impl> UniformlySampleableHabitat for WrappingN impl SingletonDemesHabitat for WrappingNoiseHabitat {} -const U32_MAX_AS_F64: f64 = (u32::MAX as f64) + 1.0_f64; - // Adapted from Christian Maher's article "Working with Simplex Noise" // Licensed under CC BY 3.0 // Published at https://cmaher.github.io/posts/working-with-simplex-noise/ @@ -223,15 +222,26 @@ fn sum_noise_octaves( scale: PositiveUnitF64, octaves: NonZeroUsize, ) -> f64 { + const F64_2_32: f64 = (u32::MAX as f64) + 1.0_f64; + let mut max_amplitude = 0.0_f64; let mut amplitude = 1.0_f64; - let mut frequency = scale.get(); + // Pre-scale the frequency to avoid pow2 degeneracies + let mut frequency = scale.get() * core::f64::consts::FRAC_1_PI * 3.0; let mut result = 0.0_f64; for _ in 0..octaves.get() { - let (x, y, z, w) = location_to_wrapping_4d::(location, frequency); - result += noise.eval_4d::(x, y, z, w) * amplitude; + // Ensure wrapping occurs at an integer boundary + let wrap = M::round(F64_2_32 * frequency); + let fixed_frequency = wrap / F64_2_32; + + let (x, y) = ( + f64::from(location.x()) * fixed_frequency, + f64::from(location.y()) * fixed_frequency, + ); + + result += noise.eval_2d::(x, y, wrap) * amplitude; max_amplitude += amplitude; amplitude *= persistence.get(); frequency *= 2.0_f64; @@ -239,19 +249,3 @@ fn sum_noise_octaves( result / max_amplitude } - -// Adapted from JTippetts' Seamless Noise article on gamedev.net: -// https://www.gamedev.net/blog/33/entry-2138456-seamless-noise/ -fn location_to_wrapping_4d(location: &Location, scale: f64) -> (f64, f64, f64, f64) { - let s = f64::from(location.x()) / U32_MAX_AS_F64; - let t = f64::from(location.y()) / U32_MAX_AS_F64; - - let scale = scale * U32_MAX_AS_F64; - - let nx = M::cos(s * 2.0_f64 * PI) * scale; - let ny = M::cos(t * 2.0_f64 * PI) * scale; - let nz = M::sin(s * 2.0_f64 * PI) * scale; - let nw = M::sin(t * 2.0_f64 * PI) * scale; - - (nx, ny, nz, nw) -} diff --git a/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/opensimplex_noise/mod.rs b/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/opensimplex_noise/mod.rs index 811973254..d19a1a75d 100644 --- a/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/opensimplex_noise/mod.rs +++ b/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/opensimplex_noise/mod.rs @@ -7,6 +7,8 @@ // OpenSimplex noise implementation vendored from // https://github.com/Mapet13/opensimplex_noise_rust // by Jakub Sordyl, licensed under the Unlicense +// +// Adapted to use MathsCore and support wrapping use necsim_core_maths::MathsCore; @@ -47,19 +49,20 @@ impl OpenSimplexNoise { #[allow(dead_code)] #[must_use] - pub fn eval_2d(&self, x: f64, y: f64) -> f64 { - OpenSimplexNoise2D::eval::(Vec2::new(x, y), &self.perm) + pub fn eval_2d(&self, x: f64, y: f64, wrap: f64) -> f64 { + OpenSimplexNoise2D::eval::(Vec2::new(x, y), &self.perm, wrap) } #[allow(dead_code)] #[must_use] - pub fn eval_3d(&self, x: f64, y: f64, z: f64) -> f64 { - OpenSimplexNoise3D::eval::(Vec3::new(x, y, z), &self.perm) + pub fn eval_3d(&self, x: f64, y: f64, z: f64, wrap: f64) -> f64 { + OpenSimplexNoise3D::eval::(Vec3::new(x, y, z), &self.perm, wrap) } + #[allow(dead_code)] #[must_use] - pub fn eval_4d(&self, x: f64, y: f64, z: f64, w: f64) -> f64 { - OpenSimplexNoise4D::eval::(Vec4::new(x, y, z, w), &self.perm) + pub fn eval_4d(&self, x: f64, y: f64, z: f64, w: f64, wrap: f64) -> f64 { + OpenSimplexNoise4D::eval::(Vec4::new(x, y, z, w), &self.perm, wrap) } } @@ -67,8 +70,8 @@ pub trait NoiseEvaluator> { const STRETCH_POINT: T; const SQUISH_POINT: T; - fn eval(point: T, perm: &PermTable) -> f64; - fn extrapolate(grid: T, delta: T, perm: &PermTable) -> f64; + fn eval(point: T, perm: &PermTable, wrap: f64) -> f64; + fn extrapolate(grid: T, delta: T, perm: &PermTable, wrap: f64) -> f64; } fn generate_perm_array(seed: i64) -> PermTable { diff --git a/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/opensimplex_noise/open_simplex_noise_2d.rs b/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/opensimplex_noise/open_simplex_noise_2d.rs index 0841b66e2..b2b0b9363 100644 --- a/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/opensimplex_noise/open_simplex_noise_2d.rs +++ b/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/opensimplex_noise/open_simplex_noise_2d.rs @@ -11,29 +11,54 @@ const SQUISH: f64 = 0.366_025_403_784_439; // (sqrt(2 + 1) - 1) / 2 const NORMALIZING_SCALAR: f64 = 47.0; -const GRAD_TABLE: [Vec2; 8] = [ - Vec2::new(5.0, 2.0), - Vec2::new(2.0, 5.0), - Vec2::new(-5.0, 2.0), - Vec2::new(-2.0, 5.0), - Vec2::new(5.0, -2.0), - Vec2::new(2.0, -5.0), - Vec2::new(-5.0, -2.0), - Vec2::new(-2.0, -5.0), -]; - pub enum OpenSimplexNoise2D {} impl NoiseEvaluator> for OpenSimplexNoise2D { const SQUISH_POINT: Vec2 = Vec2::new(SQUISH, SQUISH); const STRETCH_POINT: Vec2 = Vec2::new(STRETCH, STRETCH); - fn extrapolate(grid: Vec2, delta: Vec2, perm: &PermTable) -> f64 { - let point = GRAD_TABLE[Self::get_grad_table_index(grid, perm)]; + fn extrapolate( + grid: Vec2, + delta: Vec2, + perm: &PermTable, + wrap: f64, + ) -> f64 { + // Wrap the grid to put in the range [0; wrap), then snap to grid points + let grid = grid.map(|i| i - M::floor(i / wrap) * wrap).map(M::floor); + + let idx = Self::get_grad_table_index(grid, perm); + + // Calculate the gradient directly to avoid a GRAD_TABLE lookup + // const GRAD_TABLE: [Vec2; 8] = [ + // Vec2::new(5.0, 2.0), + // Vec2::new(2.0, 5.0), + // Vec2::new(-5.0, 2.0), + // Vec2::new(-2.0, 5.0), + // Vec2::new(5.0, -2.0), + // Vec2::new(2.0, -5.0), + // Vec2::new(-5.0, -2.0), + // Vec2::new(-2.0, -5.0), + // ]; + let mut point = if (idx & 1) == 0 { + Vec2::new(5.0, 2.0) + } else { + Vec2::new(2.0, 5.0) + }; + if (idx & 2) != 0 { + point.x *= -1.0; + } + if (idx & 4) != 0 { + point.y *= -1.0; + } + + // let point = GRAD_TABLE[Self::get_grad_table_index(grid, perm)]; point.x * delta.x + point.y * delta.y } - fn eval(input: Vec2, perm: &PermTable) -> f64 { + fn eval(input: Vec2, perm: &PermTable, wrap: f64) -> f64 { + // Pre-squish the input to allow wrapping in extrapolate + let input = input + (Self::SQUISH_POINT * input.sum()); + let stretch: Vec2 = input + (Self::STRETCH_POINT * input.sum()); let grid = stretch.map(M::floor); @@ -41,14 +66,27 @@ impl NoiseEvaluator> for OpenSimplexNoise2D { let ins = stretch - grid; let origin = input - squashed; - OpenSimplexNoise2D::get_value(grid, origin, ins, perm) + OpenSimplexNoise2D::get_value::(grid, origin, ins, perm, wrap) } } impl OpenSimplexNoise2D { - fn get_value(grid: Vec2, origin: Vec2, ins: Vec2, perm: &PermTable) -> f64 { + #[cfg_attr(target_os = "cuda", inline)] + fn get_value( + grid: Vec2, + origin: Vec2, + ins: Vec2, + perm: &PermTable, + wrap: f64, + ) -> f64 { let contribute = |x, y| -> f64 { - utils::contribute::>(Vec2::new(x, y), origin, grid, perm) + utils::contribute::, M>( + Vec2::new(x, y), + origin, + grid, + perm, + wrap, + ) }; let value = contribute(1.0, 0.0) diff --git a/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/opensimplex_noise/open_simplex_noise_3d.rs b/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/opensimplex_noise/open_simplex_noise_3d.rs index bb860be89..11f9b4835 100644 --- a/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/opensimplex_noise/open_simplex_noise_3d.rs +++ b/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/opensimplex_noise/open_simplex_noise_3d.rs @@ -44,13 +44,24 @@ impl NoiseEvaluator> for OpenSimplexNoise3D { const SQUISH_POINT: Vec3 = Vec3::new(SQUISH, SQUISH, SQUISH); const STRETCH_POINT: Vec3 = Vec3::new(STRETCH, STRETCH, STRETCH); - fn extrapolate(grid: Vec3, delta: Vec3, perm: &PermTable) -> f64 { + fn extrapolate( + grid: Vec3, + delta: Vec3, + perm: &PermTable, + wrap: f64, + ) -> f64 { + // Wrap the grid to put in the range [0; wrap), then snap to grid points + let grid = grid.map(|i| i - M::floor(i / wrap) * wrap).map(M::floor); + let point = GRAD_TABLE[Self::get_grad_table_index(grid, perm)]; point.x * delta.x + point.y * delta.y + point.z * delta.z } - fn eval(input: Vec3, perm: &PermTable) -> f64 { + fn eval(input: Vec3, perm: &PermTable, wrap: f64) -> f64 { + // Pre-squish the input to allow wrapping in extrapolate + let input = input + (Self::SQUISH_POINT * input.sum()); + let stretch: Vec3 = input + (Self::STRETCH_POINT * input.sum()); let grid = stretch.map(M::floor); @@ -58,18 +69,25 @@ impl NoiseEvaluator> for OpenSimplexNoise3D { let ins = stretch - grid; let origin = input - squashed; - Self::get_value(grid, origin, ins, perm) + Self::get_value::(grid, origin, ins, perm, wrap) } } impl OpenSimplexNoise3D { - fn get_value(grid: Vec3, origin: Vec3, ins: Vec3, perm: &PermTable) -> f64 { + fn get_value( + grid: Vec3, + origin: Vec3, + ins: Vec3, + perm: &PermTable, + wrap: f64, + ) -> f64 { let contribute = |x: f64, y: f64, z: f64| { - utils::contribute::>( + utils::contribute::, M>( Vec3::new(x, y, z), origin, grid, perm, + wrap, ) }; diff --git a/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/opensimplex_noise/open_simplex_noise_4d.rs b/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/opensimplex_noise/open_simplex_noise_4d.rs index b1513bdee..fb330f108 100644 --- a/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/opensimplex_noise/open_simplex_noise_4d.rs +++ b/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/opensimplex_noise/open_simplex_noise_4d.rs @@ -84,13 +84,24 @@ impl NoiseEvaluator> for OpenSimplexNoise4D { const SQUISH_POINT: Vec4 = Vec4::new(SQUISH, SQUISH, SQUISH, SQUISH); const STRETCH_POINT: Vec4 = Vec4::new(STRETCH, STRETCH, STRETCH, STRETCH); - fn extrapolate(grid: Vec4, delta: Vec4, perm: &PermTable) -> f64 { + fn extrapolate( + grid: Vec4, + delta: Vec4, + perm: &PermTable, + wrap: f64, + ) -> f64 { + // Wrap the grid to put in the range [0; wrap), then snap to grid points + let grid = grid.map(|i| i - M::floor(i / wrap) * wrap).map(M::floor); + let point = GRAD_TABLE[Self::get_grad_table_index(grid, perm)]; point.x * delta.x + point.y * delta.y + point.z * delta.z + point.w * delta.w } - fn eval(input: Vec4, perm: &PermTable) -> f64 { + fn eval(input: Vec4, perm: &PermTable, wrap: f64) -> f64 { + // Pre-squish the input to allow wrapping in extrapolate + let input = input + (Self::SQUISH_POINT * input.sum()); + let stretch: Vec4 = input + (Self::STRETCH_POINT * input.sum()); let grid = stretch.map(M::floor); @@ -98,7 +109,7 @@ impl NoiseEvaluator> for OpenSimplexNoise4D { let ins = stretch - grid; let origin = input - squashed; - Self::get_value(grid, origin, ins, perm) + Self::get_value::(grid, origin, ins, perm, wrap) } } @@ -657,13 +668,20 @@ impl OpenSimplexNoise4D { + contribute(0.0, 0.0, 1.0, 1.0) } - fn get_value(grid: Vec4, origin: Vec4, ins: Vec4, perm: &PermTable) -> f64 { + fn get_value( + grid: Vec4, + origin: Vec4, + ins: Vec4, + perm: &PermTable, + wrap: f64, + ) -> f64 { let contribute = |x: f64, y: f64, z: f64, w: f64| { - utils::contribute::>( + utils::contribute::, M>( Vec4::new(x, y, z, w), origin, grid, perm, + wrap, ) }; diff --git a/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/opensimplex_noise/utils.rs b/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/opensimplex_noise/utils.rs index 5c9fd568f..85e3ca80b 100644 --- a/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/opensimplex_noise/utils.rs +++ b/necsim/impls/no-std/src/cogs/habitat/wrapping_noise/opensimplex_noise/utils.rs @@ -1,10 +1,13 @@ +use necsim_core_maths::MathsCore; + use super::{vector::VecType, NoiseEvaluator, PermTable}; -pub fn contribute, Vec: VecType>( +pub fn contribute, Vec: VecType, M: MathsCore>( delta: Vec, origin: Vec, grid: Vec, perm: &PermTable, + wrap: f64, ) -> f64 { let shifted: Vec = origin - delta - NoiseEvaluatorType::SQUISH_POINT * delta.sum(); let attn: f64 = 2.0 - shifted.get_attenuation_factor(); @@ -13,5 +16,8 @@ pub fn contribute, Vec: VecType>( return 0.0; } - attn * attn * attn * attn * NoiseEvaluatorType::extrapolate(grid + delta, shifted, perm) + let attn2 = attn * attn; + let attn4 = attn2 * attn2; + + attn4 * NoiseEvaluatorType::extrapolate::(grid + delta, shifted, perm, wrap) } diff --git a/necsim/impls/no-std/src/cogs/origin_sampler/decomposition.rs b/necsim/impls/no-std/src/cogs/origin_sampler/decomposition.rs index d03bd99db..518e26a61 100644 --- a/necsim/impls/no-std/src/cogs/origin_sampler/decomposition.rs +++ b/necsim/impls/no-std/src/cogs/origin_sampler/decomposition.rs @@ -84,12 +84,7 @@ impl<'d, M: MathsCore, O: UntrustedOriginSampler<'d, M>, D: Decomposition= self - .origin_sampler - .habitat() - .get_habitat_at_location(lineage.indexed_location.location()) + .is_indexed_location_habitable(&lineage.indexed_location) { if self.decomposition.get_subdomain().is_root() { return Some(lineage); diff --git a/necsim/impls/no-std/src/cogs/speciation_probability/spatially_implicit.rs b/necsim/impls/no-std/src/cogs/speciation_probability/spatially_implicit.rs index 21d4f81c8..d50e77707 100644 --- a/necsim/impls/no-std/src/cogs/speciation_probability/spatially_implicit.rs +++ b/necsim/impls/no-std/src/cogs/speciation_probability/spatially_implicit.rs @@ -36,17 +36,15 @@ impl SpeciationProbability> for SpatiallyImplicitSpeciationProbability { #[must_use] - #[debug_requires( - habitat.local().contains(location) || habitat.meta().contains(location), - "location is inside either the local or meta habitat extent" - )] #[inline] fn get_speciation_probability_at_location( &self, location: &Location, habitat: &SpatiallyImplicitHabitat, ) -> ClosedUnitF64 { - if habitat.local().contains(location) { + // By PRE, location must be habitable, i.e. either in the local or the meta + // habitat + if habitat.local().get_extent().contains(location) { ClosedUnitF64::zero() } else { self.meta_speciation_probability.into() diff --git a/necsim/impls/no-std/src/decomposition/mod.rs b/necsim/impls/no-std/src/decomposition/mod.rs index 8c6a8186d..7791d99c5 100644 --- a/necsim/impls/no-std/src/decomposition/mod.rs +++ b/necsim/impls/no-std/src/decomposition/mod.rs @@ -14,10 +14,7 @@ pub mod radial; pub trait Decomposition>: Backup + Sized + core::fmt::Debug { fn get_subdomain(&self) -> Partition; - #[debug_requires( - habitat.contains(location), - "location is contained inside habitat" - )] + #[debug_requires(habitat.is_location_habitable(location), "location is habitable")] #[debug_ensures( ret < self.get_subdomain().size().get(), "subdomain rank is in range [0, self.get_subdomain().size())" diff --git a/necsim/plugins/species/src/identity.rs b/necsim/plugins/species/src/identity.rs index 01673f5dc..8de1999d9 100644 --- a/necsim/plugins/species/src/identity.rs +++ b/necsim/plugins/species/src/identity.rs @@ -8,7 +8,7 @@ use necsim_core::{ landscape::{IndexedLocation, Location}, lineage::GlobalLineageReference, }; -use necsim_core_bond::{NonZeroOneU64, PositiveF64}; +use necsim_core_bond::PositiveF64; #[allow(clippy::module_name_repetitions)] #[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord, Hash)] @@ -50,9 +50,9 @@ impl SpeciesIdentity { lineage: GlobalLineageReference, anchor: GlobalLineageReference, ) -> SpeciesIdentity { - let lineage = unsafe { lineage.into_inner().get() - 2 }; + let lineage = unsafe { lineage.into_inner() }; let marker = 0x0; - let anchor = unsafe { anchor.into_inner().get() - 2 }; + let anchor = unsafe { anchor.into_inner() }; Self::from_raw(lineage, marker, anchor) } @@ -92,14 +92,8 @@ impl SpeciesIdentity { return Err(self); } - let lineage = match NonZeroOneU64::new(lineage.wrapping_add(2)) { - Ok(lineage) => unsafe { GlobalLineageReference::from_inner(lineage) }, - Err(_) => return Err(self), - }; - let anchor = match NonZeroOneU64::new(anchor.wrapping_add(2)) { - Ok(anchor) => unsafe { GlobalLineageReference::from_inner(anchor) }, - Err(_) => return Err(self), - }; + let lineage = unsafe { GlobalLineageReference::from_inner(lineage) }; + let anchor = unsafe { GlobalLineageReference::from_inner(anchor) }; Ok((lineage, anchor)) } @@ -237,7 +231,7 @@ mod tests { landscape::{IndexedLocation, Location}, lineage::GlobalLineageReference, }; - use necsim_core_bond::{NonZeroOneU64, PositiveF64}; + use necsim_core_bond::PositiveF64; use super::SpeciesIdentity; @@ -274,27 +268,8 @@ mod tests { let mut rng = StdRng::from_entropy(); for _ in 0..1_000_000 { - let lineage = loop { - let l = rng.next_u64(); - - match NonZeroOneU64::new(l) { - Ok(l) => break unsafe { GlobalLineageReference::from_inner(l) }, - Err(_) => continue, - } - }; - - let anchor = loop { - let a = rng.next_u64(); - - if a > (u64::MAX >> 1) { - continue; - } - - match NonZeroOneU64::new(a) { - Ok(a) => break unsafe { GlobalLineageReference::from_inner(a) }, - Err(_) => continue, - } - }; + let lineage = unsafe { GlobalLineageReference::from_inner(rng.next_u64()) }; + let anchor = unsafe { GlobalLineageReference::from_inner(rng.next_u64()) }; let identity = SpeciesIdentity::from_unspeciated(lineage.clone(), anchor.clone()); diff --git a/necsim/plugins/species/src/individual/feather/dataframe.rs b/necsim/plugins/species/src/individual/feather/dataframe.rs index 57d0d8ff9..52c56b364 100644 --- a/necsim/plugins/species/src/individual/feather/dataframe.rs +++ b/necsim/plugins/species/src/individual/feather/dataframe.rs @@ -118,7 +118,7 @@ impl IndividualSpeciesFeatherReporter { let mut parents = Vec::with_capacity(self.origins.len()); for (lineage, origin) in &self.origins { - ids.push(unsafe { lineage.clone().into_inner().get() - 2 }); + ids.push(unsafe { lineage.clone().into_inner() }); xs.push(origin.location().x()); ys.push(origin.location().y()); @@ -130,8 +130,6 @@ impl IndividualSpeciesFeatherReporter { .unwrap_or(lineage) .clone() .into_inner() - .get() - - 2 }); } diff --git a/necsim/plugins/species/src/individual/feather/mod.rs b/necsim/plugins/species/src/individual/feather/mod.rs index 4f2154999..845834115 100644 --- a/necsim/plugins/species/src/individual/feather/mod.rs +++ b/necsim/plugins/species/src/individual/feather/mod.rs @@ -12,7 +12,7 @@ use necsim_core::{ landscape::{IndexedLocation, Location}, lineage::GlobalLineageReference, }; -use necsim_core_bond::{NonNegativeF64, NonZeroOneU64}; +use necsim_core_bond::NonNegativeF64; use crate::{LastEventState, SpeciesIdentity}; @@ -191,19 +191,13 @@ impl<'de> Deserialize<'de> for IndividualSpeciesFeatherReporter { .zip(parents.values_iter()) .zip(species.iter()) { - let id = unsafe { - GlobalLineageReference::from_inner(NonZeroOneU64::new_unchecked(*id + 2)) - }; + let id = unsafe { GlobalLineageReference::from_inner(*id) }; // Populate the individual `origins` lookup self_origins .insert(id.clone(), IndexedLocation::new(Location::new(*x, *y), *i)); - let parent = unsafe { - GlobalLineageReference::from_inner(NonZeroOneU64::new_unchecked( - *parent + 2, - )) - }; + let parent = unsafe { GlobalLineageReference::from_inner(*parent) }; // Populate the individual `parents` lookup // parent == id -> individual does NOT have a parent diff --git a/necsim/plugins/species/src/individual/sqlite/database.rs b/necsim/plugins/species/src/individual/sqlite/database.rs index 11fbf2565..804886383 100644 --- a/necsim/plugins/species/src/individual/sqlite/database.rs +++ b/necsim/plugins/species/src/individual/sqlite/database.rs @@ -2,7 +2,7 @@ use necsim_core::{ landscape::{IndexedLocation, Location}, lineage::GlobalLineageReference, }; -use necsim_core_bond::{NonZeroOneU64, PositiveF64}; +use necsim_core_bond::PositiveF64; use rusqlite::{named_params, types::Value}; @@ -240,9 +240,7 @@ impl IndividualSpeciesSQLiteReporter { let parent: Option = row.get("parent")?; let species: Option = row.get("species")?; - let id = unsafe { - GlobalLineageReference::from_inner(NonZeroOneU64::new_unchecked(from_i64(id) + 2)) - }; + let id = unsafe { GlobalLineageReference::from_inner(from_i64(id)) }; // Populate the individual `origins` lookup self.origins.insert( @@ -251,11 +249,7 @@ impl IndividualSpeciesSQLiteReporter { ); if let Some(parent) = parent { - let parent = unsafe { - GlobalLineageReference::from_inner(NonZeroOneU64::new_unchecked( - from_i64(parent) + 2, - )) - }; + let parent = unsafe { GlobalLineageReference::from_inner(from_i64(parent)) }; // Populate the individual `parents` lookup self.parents.insert(id.clone(), parent); @@ -352,14 +346,14 @@ impl IndividualSpeciesSQLiteReporter { // Positional parameters boost performance insertion.execute(rusqlite::params![ - /* :id */ to_i64(unsafe { lineage.clone().into_inner().get() - 2 }), + /* :id */ to_i64(unsafe { lineage.clone().into_inner() }), /* :x */ to_i32(origin.location().x()), /* :y */ to_i32(origin.location().y()), /* :i */ to_i32(origin.index()), /* :parent */ self.parents .get(&lineage) - .map(|parent| to_i64(unsafe { parent.clone().into_inner().get() - 2 })), + .map(|parent| to_i64(unsafe { parent.clone().into_inner() })), /* :species */ self.species .get(&ancestor) diff --git a/necsim/plugins/tskit/src/tree/metadata.rs b/necsim/plugins/tskit/src/tree/metadata.rs index ffb883659..b3cc23375 100644 --- a/necsim/plugins/tskit/src/tree/metadata.rs +++ b/necsim/plugins/tskit/src/tree/metadata.rs @@ -1,10 +1,5 @@ -use std::{ - array::TryFromSliceError, - convert::{TryFrom, TryInto}, - io, -}; +use std::{array::TryFromSliceError, convert::TryInto, io}; -use necsim_core_bond::NonZeroOneU64; use tskit::metadata::{IndividualMetadata, MetadataError, MetadataRoundtrip, NodeMetadata}; use necsim_core::lineage::GlobalLineageReference; @@ -15,8 +10,8 @@ pub struct GlobalLineageMetadata(GlobalLineageReference); impl MetadataRoundtrip for GlobalLineageMetadata { fn encode(&self) -> Result, MetadataError> { - // Store the internal u64 without the +2 offset - Ok((unsafe { self.0.clone().into_inner() }.get() - 2) + // Store the internal u64 + Ok(unsafe { self.0.clone().into_inner() } .to_le_bytes() .to_vec()) } @@ -32,17 +27,9 @@ impl MetadataRoundtrip for GlobalLineageMetadata { } })?; - // Convert the bytes into an u64 with the needed +2 offset - let value = u64::from_le_bytes(value_bytes) + 2; - - // Create the internal `NonZeroOneU64` representation of the reference - let value_inner = - NonZeroOneU64::try_from(value).map_err(|err| MetadataError::RoundtripError { - value: Box::new(io::Error::new(io::ErrorKind::InvalidData, err.to_string())), - })?; - + // Convert the bytes into an u64 and a GlobalLineageReference Ok(Self(unsafe { - GlobalLineageReference::from_inner(value_inner) + GlobalLineageReference::from_inner(u64::from_le_bytes(value_bytes)) })) } } diff --git a/rustcoalescence/algorithms/cuda/src/parallelisation/monolithic.rs b/rustcoalescence/algorithms/cuda/src/parallelisation/monolithic.rs index d687ff99a..66e1ff479 100644 --- a/rustcoalescence/algorithms/cuda/src/parallelisation/monolithic.rs +++ b/rustcoalescence/algorithms/cuda/src/parallelisation/monolithic.rs @@ -314,7 +314,7 @@ pub fn simulate< } } - event_buffer.report_events(&mut proxy); + event_buffer.report_events_unordered(&mut proxy); proxy.local_partition().get_reporter().report_progress( &((slow_lineages.len() as u64) + (fast_lineages.len() as u64)).into(),