From a90647273f4489d7ca9363060fa8bc01b94acfe1 Mon Sep 17 00:00:00 2001
From: hopeyen <hopey1118@gmail.com>
Date: Mon, 17 Apr 2023 19:21:58 +0200
Subject: [PATCH] fix: force exact numeric sum - simplex projection

Signed-off-by: hopeyen <hopey1118@gmail.com>
---
 src/project.jl  | 1 +
 test/project.jl | 4 ++++
 2 files changed, 5 insertions(+)

diff --git a/src/project.jl b/src/project.jl
index 192badc..21ab37b 100644
--- a/src/project.jl
+++ b/src/project.jl
@@ -16,6 +16,7 @@ function σsimplex(x::AbstractVector{T}, σ::Real) where {T<:Real}
     ρ = maximum((1:n)[μ - (cumsum(μ) .- σ) ./ (1:n) .> zero(T)])
     θ = (sum(μ[1:ρ]) - σ) / ρ
     w = max.(x .- θ, zero(T))
+    w[end] += σ - sum(w)
     return w
 end
 
diff --git a/test/project.jl b/test/project.jl
index f2735a5..18ae941 100644
--- a/test/project.jl
+++ b/test/project.jl
@@ -14,6 +14,10 @@
         x = Float64[-1, 0, 0]
         σ = 5
         @test σsimplex(x, σ) ≈ [1, 2, 2]  # Scale up
+        
+        x = Float64[23133337391432116]
+        σ = 652174.7265297174
+        @test sum(σsimplex(x, σ)) ≈ σ  # exact
     end
 
     @testset "gssp" begin