Add resamplingprobability function

Author: brenhinkeller

src/utilities/Resampling.jl

@@ -1077,7 +1077,7 @@
             @batch for g ∈ eachindex(spatialscaleᵣ)
                 for h ∈ eachindex(agescale)
                     kᵢ = 0.0
-                    @inbounds for j ∈ 1:N
+                    @inbounds @simd ivdep for j ∈ 1:N
                         kᵢⱼ = 1.0 / ((spatialdistᵣ[j]/spatialscaleᵣ[g])^lp + 1.0) + 1.0 / ((abs(age[i] - age[j])/agescale[h])^lp + 1.0)
                         kᵢ += kᵢⱼ
                     end
@@ -1175,18 +1175,12 @@
         k = Array{Float64}(undef, N)
         p = Progress(N÷10, desc="Calculating weights: ")
         @inbounds @batch for i ∈ 1:N
-            if isnan(x[i])
-                # If there is missing data, set k=inf for weight=0
-                k[i] = Inf
-            else
-                # Otherwise, calculate weight
-                kᵢ = 0.0
-                @turbo for j ∈ 1:N
-                    kᵢⱼ = 1.0 / ( (abs(x[i] - x[j])/xscale)^lp + 1.0)
-                    kᵢ += kᵢⱼ
-                end
-                k[i] = kᵢ
+            kᵢ = 0.0
+            @turbo for j ∈ 1:N
+                kᵢⱼ = 1.0 / ( (abs(x[i] - x[j])/xscale)^lp + 1.0)
+                kᵢ += kᵢⱼ
             end
+            k[i] = kᵢ
             (i % 10 == 0) && next!(p)
         end
         return k
@@ -1196,15 +1190,33 @@
 
     """
     ```julia
-    k = invweight_age(age::AbstractArray; lp::Number=2, agescale::Number=38.0)
+    k = invweight_age(age::AbstractArray;
+        \tlp::Number=2, 
+        \tagescale::Number=38.0
+    )
     ```
 
     Find the inverse weights `k` (proportional to temporal sample density) for
     a set of geological samples with specified `age` (of crystallization, deposition, etc.).
     """
-    function invweight_age(age::AbstractArray; lp::Number=2, agescale::Number=38.0)
-        return invweight(age, agescale, lp=lp)
-    end
+    invweight_age(age::AbstractArray; lp::Number=2, agescale::Number=38.0) = invweight(age, agescale; lp)
     export invweight_age
 
+
+    """
+    ```julia
+    resamplingprobability(k; median_probability::Number=1/6)
+    ```
+
+    Calculate scaled resampling probabilities `p` given a vector of inverse weights `k` 
+    and a desired median resampling probability (1/6 by default).
+    """
+    function resamplingprobability(k; median_probability::Number=1/6)
+        median_k = 1/median_probability - 1
+        f = median_k / nanmedian(filter(isfinite, k))
+        p = @. 1 / ((k * f) + 1)
+        return p
+    end
+    export resamplingprobability
+
 ## --- End of file

test/testResampling.jl

@@ -103,6 +103,9 @@
     @test invweight_location(lat, lon) ≈ [2.3478642777118957, 2.950680392065609, 3.2200556525889006, 3.348975353894235, 3.4103053158720016, 3.4297605864726126, 3.413776296513463, 3.3555688532471923, 3.2289225181937002, 2.9601916238626513, 2.3566734930529947, Inf]
     @test invweight_age(age) ≈ [10.744908100409873, 10.808982898265619, 10.859335221424754, 10.895581458515855, 10.91744238026025, 10.924748324734335, 10.91744238026025, 10.895581458515855, 10.859335221424754, 10.808982898265619, 10.744908100409873, Inf]
 
+    @test resamplingprobability(invweight(lat, lon, age)) ≈ [0.17709468716990867, 0.16997032242046156, 0.1667541264064975, 0.16513990970091127, 0.16433863731114293, 0.1640825308282905, 0.16430537469842987, 0.16507611788372856, 0.16666666666666666, 0.16987285789479709, 0.17699668840192856, 0.0]
+    @test resamplingprobability(invweight_location(lat, lon)) ≈ [0.2157181923483138, 0.1795609083825818, 0.16704894086177344, 0.1616579035210186, 0.15921356093554181, 0.158453529370978, 0.1590774311384848, 0.16139152082816496, 0.16666666666666666, 0.17908729309480173, 0.21508527491483678, 0.0]
+    @test resamplingprobability(invweight_age(age)) ≈ [0.16814313324769609, 0.1673131616109112, 0.16666666666666666, 0.16620436987522427, 0.16592678603901323, 0.16583422381194865, 0.16592678603901323, 0.16620436987522427, 0.16666666666666666, 0.1673131616109112, 0.16814313324769609, 0.0]
 
 ## --- bin_bsr