add phenofit_CA-NS6

data/data_CA-NS6.jl

@@ -0,0 +1,29 @@
+using JLD2
+using RData
+using UnPack
+using VegCurveFit
+using BenchmarkTools
+using Test
+
+l = load("data/phenofit-CA_NS6.rda")
+@unpack INPUT, brks2 = l
+
+dt = brks2["dt"]
+nptperyear = INPUT["nptperyear"] |> Int
+y = INPUT["y"]
+w = INPUT["w"]
+t = INPUT["t"]
+
+serialize("data/phenofit_CA-NS6", (;y, w, t))
+
+# QC_flag = INPUT["QC_flag"]
+n = length(y)
+w = ones(n)
+t = 1:n
+
+# ylu = INPUT["ylu"]
+# ylu = [ylu[1], ylu[2]]
+# using RTableTools
+# d = DT(; y, w)
+# fwrite(d, "data.csv")
+# jldsave("data/phenofit-CA_NS6.jld2", true; y, t, w, QC_flag, ylu, nptperyear, dt)

data/test_CA-NS6.jl

@@ -0,0 +1,41 @@
+using VegCurveFit
+using BenchmarkTools
+using Test
+
+y, w, t = deserialize("data/phenofit_CA-NS6")
+
+lambda = 2.0
+include_cve = false
+
+@testset "whit3" begin
+  @time z, cve = whit3(y, w; lambda, include_cve)
+  @time z2, cve2 = WHIT(y, w; d=3, lambda=2, include_cve)
+
+  @test maximum(z - z2) <= 1e-10
+  @test abs(cve - cve2) <= 1e-10
+end
+
+@testset "whit2" begin
+  @time z, cve = whit2(y, w; lambda, include_cve)
+  @time z2, cve2 = WHIT(y, w; d=2, lambda=2, include_cve)
+
+  @test maximum(z - z2) <= 1e-10
+  @test abs(cve - cve2) <= 1e-10
+end
+
+## whit3
+include_cve = true
+@btime z, cve = whit3($y, $w; lambda, include_cve);
+@btime z2, cve2 = WHIT($y, $w; d=3, lambda=2, include_cve);
+
+## whit2
+@btime z, cve = whit2($y, $w; lambda, include_cve);
+@btime z2, cve2 = WHIT($y, $w; d=2, lambda=2, include_cve);
+
+begin
+  using Plots
+  plot()
+  # plot(y, label="y", color="grey")
+  plot!(z, label="fast")
+  plot!(z2, label="matrix")
+end

src/Smooth/Whittaker/WHIT.jl

@@ -33,8 +33,8 @@ function WHIT(y::AbstractVector, w::AbstractVector, x::AbstractVector;
   L = cholesky(A, perm=1:n).L # sparse
   z = L' \ (L \ (w.*y))
   # display(sparse(L)[1:7, 1:7])
-  
-  cve = 999.0
+
+  cve = -999.0
   if include_cve
     inv_L = Matrix(sparse(L))^-1 # 这一步可能消耗了较多的时间
     H = inv_L' * inv_L * W # 借力cholesky