fix bugs

LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2024 Dongdong Kong
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

README.md

@@ -0,0 +1,15 @@
+# ClimateStats.jl
+
+<!-- [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://jl-pkgs.github.io/ClimateStats.jl/stable) -->
+[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://jl-pkgs.github.io/ClimateStats.jl/dev)
+[![CI](https://github.com/jl-pkgs/ClimateStats.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/jl-pkgs/ClimateStats.jl/actions/workflows/CI.yml)
+[![Codecov](https://codecov.io/gh/jl-pkgs/ClimateStats.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/jl-pkgs/ClimateStats.jl)
+
+> Dongdong Kong
+
+# Installation
+
+```
+using Pkg
+Pkg.add(url="https://github.com/jl-pkgs/ClimateStats.jl")
+```

docs/make.jl

@@ -15,11 +15,7 @@ makedocs(
   ),
   pages=[
     "Introduction" => "index.md"
-    "R Base"       => "RBase.md"
-    "Statistics"   => "Statistics.md"
-    "Slope"        => "Slope.md"
     "Climate"      => "Climate.md"
-    "Parallel"     => "Parallel.md"
   ],
   sitename="ClimateStats.jl",
   warnonly=true,

docs/src/Climate.md

@@ -25,8 +25,8 @@ cal_climatology_full
 
 
 ```@docs
-Ipaper.cal_mTRS_base!
-Ipaper.cal_mTRS_base3!
+ClimateStats.cal_mTRS_base!
+ClimateStats.cal_mTRS_base3!
 ```
 
 

src/ClimateStats.jl

@@ -1,12 +1,14 @@
 module ClimateStats
 
-using Ipaper
-using Printf
-using Dates
-using NaNStatistics
 
 using Reexport
 @reexport using DocStringExtensions: TYPEDSIGNATURES, METHODLIST
+using Printf
+using Dates
+
+using NaNStatistics
+using Ipaper
+using Ipaper: movmean
 
 include("base.jl")
 include("threshold.jl")
@@ -15,13 +17,11 @@ include("anomaly.jl")
 include("warming_level.jl")
 
 # export cal_mTRS_base, cal_mTRS_season, cal_mTRS_full
+export filter_mds
 export format_md
 export cal_climatology_base, cal_climatology_full, cal_climatology_full, 
   _cal_anomaly, 
   cal_anomaly_quantile, cal_anomaly_clim, cal_threshold,
   cal_warming_level, cal_yearly_Tair
 
-
-export filter_mds
-
 end # module ClimateStats

test/Statistics/test-apply.qmd

@@ -0,0 +1,13 @@
+```{julia}
+using Ipaper
+
+## example 02
+dates = make_date(2010):Day(1):make_date(2013, 12, 31)
+n = 100
+ntime = length(dates)
+x = rand(n, n, ntime)
+
+years = year.(dates)
+res = apply(x, 3; by=years, fun=nanquantile, combine=true, probs=[0.05, 0.95])
+obj_size(res)
+```

test/Statistics/test-quantile.qmd

@@ -0,0 +1,26 @@
+```{julia}
+using BenchmarkTools
+
+ntime = 365*30
+n = 100
+x = rand(Float32, n, n, ntime);
+probs = [0.9, 0.99, 0.9999]
+
+# @time r = _nanquantile(x; probs, dims=3);
+@time r0 = nanquantile(x, probs[1]; dims=3);
+
+@time r1 = NanQuantile(x; probs, dims=3);
+@time r2 = NanQuantile_low(x; probs, dims=3);
+@time r3 = NanQuantile_3d(x; probs, dims=3)
+
+
+# r == r1
+r1 == r2
+```
+
+
+```{julia}
+@benchmark r1 = NanQuantile($x; probs, dims=3) evals = 5
+@benchmark r2 = NanQuantile_low($x; probs, dims=3) evals = 5
+@benchmark r3 = NanQuantile_3d($x; probs, dims=3) evals = 5
+```

test/Statistics/test_linreg.qmd

@@ -0,0 +1,69 @@
+```{julia}
+using Ipaper
+using BenchmarkTools
+set_seed(1)
+
+n = 30
+y = collect(1:n)
+x = rand(n)
+```
+
+```{julia}
+function linreg_simple(y::AbstractVector{T}, x::AbstractVector{T}) where {T<:Real}
+  (N = length(x)) == length(y) || throw(DimensionMismatch())
+
+  x̄ = mean(x)
+  ȳ = mean(y)
+
+  sum_a = T(0.0)
+  sum_b = T(0.0)
+  @inbounds for i = 1:N
+    sum_a += (x[i] - x̄) * (y[i] - ȳ)
+    sum_b += (x[i] - x̄) * (x[i] - x̄)
+  end
+  β1 = sum_a / sum_b
+  β0 = ȳ - β1 * x̄
+  # β1 = sum((x .- x̄) .* (y .- ȳ)) / sum((x .- x̄) .^ 2)
+  [β0, β1]
+end
+
+function linreg_simple(y::AbstractVector{T}) where {T<:Real}
+  x = Float32.(1:length(y))
+  linreg_simple(y, x)
+end
+
+function linreg_fast(y::AbstractVector{T}, x::AbstractVector{T}) where {T<:Real}
+  (N = length(x)) == length(y) || throw(DimensionMismatch())
+  ldiv!(
+    cholesky!(Symmetric([T(N) sum(x); zero(T) sum(abs2, x)], :U)),
+    [sum(y), dot(x, y)])
+end
+
+function linreg_fast(y::AbstractVector{T}) where {T<:Real}
+  x = Float32.(1:length(y))
+  linreg_fast(y, x)
+end
+# lm = linreg_fast
+```
+
+
+```{julia}
+# 比R语言的版本快80倍左右
+@benchmark R"""
+rtrend::slope_p($y, $x)
+"""
+```
+
+```{julia}
+using Random
+
+N = 10000
+x = rand(N)
+X = [ones(N) x]
+y = 10 .+ x .* 0.3
+
+# @benchmark linreg1(y, X)
+# @benchmark linreg2(y, X)
+@benchmark linreg_fast(y, x)
+@benchmark linreg_simple(y, x)
+```

test/Statistics/test_vogel2020.jl

@@ -0,0 +1,41 @@
+@time using Ipaper
+
+using nctools
+using JLD2
+
+function _ncwrite(mTRS_full, outfile="OBS_mTRS_full_V2.nc")
+  _dims = nc_dims(f)
+  probs = [0.9, 0.95, 0.99, 0.999, 0.9999]
+  dim_prob = NcDim("prob", probs)
+  dims = [_dims; dim_prob]
+
+  @time nc_write(outfile, "mTRS", mTRS_full, dims; compress=1, overwrite=true)
+  # jldsave("OBS_mTRS_full_V1.jld2"; mTRS_full, lon, lat, dates, compress=true)
+end
+
+f = "/share/Data/CN0.5.1_ChinaDaily_025x025/HI-Tmax_CN05.1_1961_2021_daily_025x025.nc"
+lon = nc_read(f, "lon")
+lat = nc_read(f, "lat")
+
+dates = nc_date(f)
+@time arr = nc_read(f);
+r_range(dates)
+# size(arr)
+# length(dates)
+# dates[[1, end]]
+
+# T_year = cal_mTRS_seasonal(arr, dates)
+
+# `mTRS_full`最耗时
+Threads.nthreads()
+
+if false
+  @time mTRS_full = cal_mTRS_full(arr, dates; use_mov=true);
+  _ncwrite(mTRS_full, "OBS_mTRS_full_V2.nc")
+end
+# 1675.866281 seconds (8.29 G allocations: 3.986 TiB, 24.24% gc time, 0.45% compilation time)
+# 0.465 hours, 16 threads
+
+# @profview_allocs mTRS_full = cal_mTRS_full(arr, dates; probs=[0.9]);
+@time mTRS_full_simple = cal_mTRS_full(arr, dates; use_mov=false);
+_ncwrite(mTRS_full_simple, "OBS_mTRS_full_simple_V2.nc")

test/Statistics/test_vogel2020.qmd

@@ -0,0 +1,70 @@
+```{julia}
+@time using Ipaper
+
+using nctools
+using JLD2
+import Ipaper: cal_mTRS_base, cal_mTRS_season, cal_mTRS_full
+
+function getDataType(data)
+  x = eltype(data)
+  typeof(x) == Union ? x.b : x
+end
+
+function replaceMiss(data; nan = NaN)
+  nan = getDataType(data)(nan)
+  replace(data, missing => nan)
+end
+
+function _ncwrite(mTRS_full, outfile="OBS_mTRS_full_V2.nc"; 
+  compress=1, overwrite=true)
+  
+  _dims = nc_dims(f)
+  # 需要修复nc_dims
+  probs = [0.9, 0.95, 0.99, 0.999, 0.9999]
+  dim_prob = NcDim("prob", probs)
+  dims = [
+    _dims["longitude"];
+    _dims["latitude"];
+    _dims["time"]; dim_prob]
+
+  @time nc_write(outfile, "mTRS", mTRS_full, dims; compress, overwrite)
+  # jldsave("OBS_mTRS_full_V1.jld2"; mTRS_full, lon, lat, dates, compress=true)
+end
+```
+
+```{julia}
+f = "Z:/DATA/China/ChinaMet_HI_RHtestsV4_v20230409/data-raw/ChinaMet_HI_RHtestsV4_v20230409_HI_max_G100.nc"
+nc = nc_open(f)
+
+# f = "/share/Data/CN0.5.1_ChinaDaily_025x025/HI-Tmax_CN05.1_1961_2021_daily_025x025.nc"
+lon = nc_read(f, "longitude")
+lat = nc_read(f, "latitude")
+
+dates = nc_date(f)
+@time arr = nc_read(f);
+arr2 = replaceMiss(arr)
+```
+
+```{julia}
+r_range(dates)
+# size(arr)
+# length(dates)
+# dates[[1, end]]
+
+# T_year = cal_mTRS_seasonal(arr, dates)
+# `mTRS_full`最耗时
+Threads.nthreads()
+
+if false
+  @time mTRS_full = cal_mTRS_full(arr, dates; use_mov=true);
+  _ncwrite(mTRS_full, "OBS_mTRS_full_V2.nc")
+end
+# 1675.866281 seconds (8.29 G allocations: 3.986 TiB, 24.24% gc time, 0.45% compilation time)
+# 0.465 hours, 16 threads
+```
+
+```{julia}
+# @profview_allocs mTRS_full = cal_mTRS_full(arr, dates; probs=[0.9]);
+@time mTRS_full_simple = cal_mTRS_full(arr2, dates; use_mov=false, na_rm=true);
+_ncwrite(mTRS_full_simple, "OBS_mTRS_full_simple_V2.nc")
+```

test/runtests.jl

@@ -1,5 +1,6 @@
 using Test
 using ClimateStats
+import ClimateStats: cal_mTRS_base3!
 using NaNStatistics
 using Dates
 using Ipaper: set_seed, obj_size

test/test-stat_threshold.jl

@@ -7,7 +7,7 @@
   A = rand(Float32, 4, 4, n)
 
   probs = [0.90, 0.95, 0.99, 0.999, 0.9999]
-  kw = (; probs, na_rm=true, parallel=false, (fun!)=Ipaper.cal_mTRS_base3!)
+  kw = (; probs, na_rm=true, parallel=false, (fun!)=cal_mTRS_base3!)
 
   @time r1 = cal_mTRS_base(A, dates; kw..., method_q="mapslices")
   @time r2 = cal_mTRS_base(A, dates; kw..., method_q="base")