From b6b6e8115434ab90fa0069c22a7aa302cb91121b Mon Sep 17 00:00:00 2001
From: Dongdong Kong <kongdd@users.noreply.github.com>
Date: Mon, 8 Apr 2024 23:13:07 +0800
Subject: [PATCH] fix pow for negative value

---
 data/x_ref              | Bin 0 -> 446 bytes
 src/Dist/logLogistic.jl |   2 +-
 src/SPEI.jl             |   9 ++++++++-
 src/main_spei.jl        |   4 ++++
 test/test-spei.jl       |  26 +++++++++++++-------------
 5 files changed, 26 insertions(+), 15 deletions(-)
 create mode 100644 data/x_ref

diff --git a/data/x_ref b/data/x_ref
new file mode 100644
index 0000000000000000000000000000000000000000..16ba63151f526980c1f30a14c5942f250a005d9c
GIT binary patch
literal 446
zcmXr_@{wR+U|=v2VC3>k%uP)RDJ{rJmG_(N!UPl*VX9?ht%xs5O%n))i3tS2SPZ-m
zbeI`L8F&rtfNDV?DzD4VxxC%(&<8=gb9~$EwB9ba`_IW@Ki^j0Zql9wc5ciLc8g9b
z*>CnaXZK-;f}KHzrmfY}e7gl+I(AzwoVLr~ImNCkklTKZ-YL6Ny_$B`E3NE8HYeJp
zrz+Wn^{lpw$X;Zpvc<#x%|9#qTiUU9zl4w5@p*IDA9=9J&O%t*-b(eJ-I3>$?LzpX
z?WRff*t3Nu+rR31VyE{~!rq%<vi;12J9e8&cG^i+H`*<YaIn{#dBpCr{9L<hjdShN
zw43d`5}w$zJ&&^gV~}q5Zoh+lP2wRtsn@-BnYGOJ)8C2OzhPQycWu)ZJ41^tc0XKJ
Q+I6va*a7|jcc9z~02)fC0{{R3

literal 0
HcmV?d00001

diff --git a/src/Dist/logLogistic.jl b/src/Dist/logLogistic.jl
index b400dd2..ce9381b 100644
--- a/src/Dist/logLogistic.jl
+++ b/src/Dist/logLogistic.jl
@@ -26,7 +26,7 @@ end
 
 function cdf_logLogistic(x::Real, params)
   β, α, γ = params[1:3]
-  1 / (1 + ((α / (x - β))^γ))
+  1 / (1 + (pow(α / (x - β), γ)))
 end
 
 
diff --git a/src/SPEI.jl b/src/SPEI.jl
index f6e0256..4782a7e 100644
--- a/src/SPEI.jl
+++ b/src/SPEI.jl
@@ -8,7 +8,14 @@ using SpecialFunctions: gamma_inc, lgamma, loggamma
 
 
 qnorm(p::Real) = quantile(Normal(), p)
-pow = ^
+
+export pow
+# pow = ^
+# import Base:^
+pow(x, y) = x^y
+function pow(x::T, y::AbstractFloat) where {T<:AbstractFloat}
+  x < 0 ? T(NaN) : x^y
+end
 
 
 include("PWM.jl")
diff --git a/src/main_spei.jl b/src/main_spei.jl
index 3c6a3bd..d7e0828 100644
--- a/src/main_spei.jl
+++ b/src/main_spei.jl
@@ -1,6 +1,8 @@
 function param_spei(x_ref::AbstractVector)
   x2 = x_ref[.!isnan.(x_ref)] |> sort # have to sort
   beta = pwm(x2, 0.0, 0.0, 0)
+  @show beta
+
   params = _fit_logLogistic(beta)
   params
 end
@@ -21,6 +23,8 @@ function param_spi(x_ref::AbstractVector; fit="lmom")
   params(D), q
 end
 
+export param_spei, param_spi
+
 
 """
     spei(x::AbstractVector)
diff --git a/test/test-spei.jl b/test/test-spei.jl
index 6f4132f..ea901eb 100644
--- a/test/test-spei.jl
+++ b/test/test-spei.jl
@@ -2,24 +2,27 @@ using Test
 using RCall
 using SPEI
 
+R"""
+library(SPEI2)
+"""
+
 function nanmax(x)
   x2 = x[.!isnan.(x)]
   maximum(x2)
 end
 
-R"""
-library(SPEI2)
-"""
-wb = R"wb" |> rcopy |> x -> Float32.(x)
+# wb = R"wb" |> rcopy |> x -> Float32.(x)
+wb = deserialize("data/wb")
+x, x_ref = deserialize("data/x_ref")
 
 @testset "spei" begin
-  @time r_R = R"cal_spei(wb)" |> rcopy
-  @time r_jl = spei(wb)
+  @time r_R = R"cal_spei($x, $x_ref)$z" |> rcopy
+  @time r_jl = spei(x, x_ref).z
 
-  z_r = r_R[:z]
-  z_jl = r_jl.z
-  e = z_r - z_jl
-  @test maximum(abs.(e)) <= 1e-3
+  # z_r = r_R[:z]
+  # z_jl = r_jl.z
+  e = r_R - r_jl
+  @test nanmax(e) <= 1e-3
 end
 
 @testset "spi" begin
@@ -33,11 +36,8 @@ end
 end
 
 # > 计算结果完全一致
-
-# ## 误差还挺大的
 # begin
 #   using Plots
-
 #   plot(z_r, label="R")
 #   plot!(z_jl, label="Julia")
 #   plot!(e, label="R - Julia")