Merge pull request #54 from ArndtLab/histogram-bug

Histogram bug

Author: TommasoStentella

src/HetDister.jl

@@ -102,23 +102,24 @@ function CustomEdgeVector(; lo = 1, hi = 10, nbins::Integer)
 end
 
 """
-    adapt_histogram(segments::AbstractVector{<:Integer}; lo::Int=1, hi::Int=50_000_000, nbins::Int=200)
+    adapt_histogram(segments::AbstractVector{<:Integer}; lo::Int=1, hi::Int=50_000_000, nbins::Int=800)
 
 Build an histogram from `segments` logbinned between `lo` and `hi`
 with `nbins` bins.
 
 The upper limit is adapted to ensure logspacing with the requested `nbins`.
 """
-function adapt_histogram(segments::AbstractVector{<:Integer}; lo::Int=1, hi::Int=50_000_000, nbins::Int=200)
+function adapt_histogram(segments::AbstractVector{<:Integer}; lo::Int=1, hi::Int=50_000_000, nbins::Int=800)
     h_obs = Histogram(CustomEdgeVector(;lo, hi, nbins))
     @assert !isempty(segments)
     append!(h_obs, segments)
-    l = findlast(h_obs.weights .> 1) + 1
-    while h_obs.edges[1].edges[l+1] + 1 < hi
+    l = findlast(h_obs.weights .> 1)
+    isnothing(l) && return h_obs
+    while h_obs.edges[1].edges[l+1] < hi
         hi = h_obs.edges[1].edges[l+1]
         h_obs = Histogram(CustomEdgeVector(;lo, hi, nbins))
         append!(h_obs, segments)
-        l = findlast(h_obs.weights .> 1) + 1
+        l = findlast(h_obs.weights .> 1)
     end
     return h_obs
 end

src/Spectra/SMCpIntegrals.jl

@@ -202,7 +202,7 @@ function prordn!(res::AbstractMatrix{<:Real}, jprt::AbstractMatrix{<:Real},
                     s += temp[k,j] * exp(-2rate * (rs[i]-edges[k+1]) * ts[j]) * (- expm1(-2rate * w * ts[j])) / 2ts[j]
                 end
                 w = edges[i+1] - edges[i]
-                if w <= 1
+                if w < 1
                     s += temp[i,j] * (- expm1(-2rate * w * ts[j])) / 2ts[j]
                 else
                     w = rs[i] - edges[i]

src/corrections.jl

@@ -65,7 +65,9 @@ function demoinfer(h_obs::Histogram{T,1,E}, epochrange::AbstractRange{<:Integer}
         fits = map(r->r.f, results),
         chains = map(r->r.chain, results),
         corrections = map(r->r.corrections, results),
-        h_obs = results[1].h_obs
+        h_obs = results[1].h_obs,
+        yth = map(r->r.yth, results),
+        deltas = map(r->r.deltas, results)
     )
 end
 
@@ -83,10 +85,12 @@ function demoinfer(h_obs::Histogram{T,1,E}, epochs::Int, fop_::FitOptions;
 
     chain = []
     corrections = []
+    deltas = []
 
     ho_mod.weights .= h_obs.weights
     corr = zeros(eltype(h_obs.weights), length(h_obs.weights))
     f = nothing
+    yth = nothing
     for iter in 1:iters
         fits = pre_fit!(fop, ho_mod, epochs; require_convergence = false)
         f = fits[end]
@@ -109,7 +113,9 @@ function demoinfer(h_obs::Histogram{T,1,E}, epochs::Int, fop_::FitOptions;
         if iter > 1
             deltacorr = (corrections[iter] .- corrections[iter-1]) ./ corrections[iter-1]
             deltacorr[isnan.(deltacorr)] .= 0.
-            if maximum(abs.(deltacorr)) < reltol
+            delta = maximum(abs.(deltacorr))
+            push!(deltas, delta)
+            if delta < reltol
                 break
             end
         end
@@ -119,7 +125,9 @@ function demoinfer(h_obs::Histogram{T,1,E}, epochs::Int, fop_::FitOptions;
         f,
         chain,
         corrections,
-        h_obs
+        h_obs,
+        yth,
+        deltas
     )
 end
 

test/runtests.jl

@@ -93,9 +93,10 @@ end
 
 @testset "Test core functionality" for (mu,rho,TN) in zip(mus, rhos, TNs)
 
-    h = Histogram(LogEdgeVector(lo = 1, hi = 1_000_000, nbins = 200))
     ibs_segments = get_sim(TN, mu, rho)
-    append!(h, ibs_segments)
+    h = adapt_histogram(ibs_segments; nbins = 200)
+    @test length(h.weights) == 200
+    @test h.weights[end] > 1
 
     stat = pre_fit(h, 2, mu, rho, 10, 100, sum(ibs_segments); require_convergence = false)
     @test isassigned(stat, 1)
@@ -113,7 +114,11 @@ end
         iters = 1
     )
     @test length(res.chains) == length(TN)÷2
+    @test length(res.yth) == length(TN)÷2
     @test all(length.(res.chains) .== 1)
+    @test all(length.(res.corrections) .== 1)
+    @test all(length.(res.deltas) .== 0)
+    @test all(length.(res.yth) .== length(res.h_obs.weights))
     @test !any(isinf.(evd.(res.fits)))
     best = compare_models(res.fits)
     @test !isnothing(best)
@@ -139,7 +144,7 @@ end
 @testset "fitting procedure" begin
     @testset "exhaustive pre-fit $(length(TN)÷2) epochs,  mu $mu, rho $rho" for (mu,rho,TN) in itr
         ibs_segments = get_sim(TN, mu, rho)
-        h = adapt_histogram(ibs_segments)
+        h = adapt_histogram(ibs_segments; nbins = 200)
         Ltot = sum(ibs_segments)
         fop = FitOptions(Ltot, mu, rho; maxnts = 8, force = false)
         fits = pre_fit!(fop, h, 8; require_convergence = false)
@@ -153,7 +158,7 @@ end
     @testset "Iterative fit" begin
         mu, rho, TN = mus[1], rhos[1], TNs[3]
         ibs_segments = get_sim(TN, mu, rho)
-        h = adapt_histogram(ibs_segments)
+        h = adapt_histogram(ibs_segments; nbins = 200)
         Ltot = sum(ibs_segments)
         pfits = pre_fit(h, 5, mu, rho, 10, 100, Ltot; require_convergence = false)
         fop = FitOptions(Ltot, mu, rho; order = 10, ndt = 800)