revise WHIT ddmat

docs/formula/main_symb.jl

@@ -1,13 +1,6 @@
 using Symbolics
 import Symbolics: scalarize, variables
 
-function tri_lower(name, n)
-  x = variables(name, 1:n, 1:n)
-  for i = 1:n, j = i+1:n
-    x[i, j] = 0
-  end
-  x
-end
 
 Mat(name, n) = variables(name, 1:n, 1:n)
 Vec(name, n) = variables(name, 1:n)
@@ -33,6 +26,25 @@ function tri_upper(name, n)
   x
 end
 
+function tri_lower(name, n)
+  x = variables(name, 1:n, 1:n)
+  for i = 1:n, j = i+1:n
+    x[i, j] = 0
+  end
+  x
+end
+
+# speye(n) = SparseArrays.sparse(I, n, n)
+# function ddmat(x::AbstractVector, d::Integer=2)
+#   m = length(x)
+#   if d == 0
+#     return speye(m)
+#   else
+#     # dx = x[(d+1):m] - x[1:(m-d)] # bug may here
+#     return diff(ddmat(x, d - 1))
+#   end
+# end
+
 @variables λ
 
 function diag_m(x)
@@ -48,14 +60,6 @@ function diag_m(x)
   M
 end
 
-
-
-
-function Base.diff(x::AbstractMatrix, d::Integer=1)
-  D = x[2:end, :] .- x[1:end-1, :]
-  d >= 2 ? diff(D, d - 1) : D
-end
-
 # 代数余子式; algebraic complement
 function complement(A::AbstractArray, i=1, j=1; verbose=false)
   i, j = j, i

docs/formula/main_whit.jl

@@ -1,90 +1,13 @@
 using LinearAlgebra, SparseArrays
 
 include("main_symb.jl")
-speye(n) = SparseArrays.sparse(I, n, n)
 
-
-function ddmat(x::AbstractVector, d::Integer=2)
-  m = length(x)
-  if d == 0
-    return speye(m)
-  else
-    # dx = x[(d+1):m] - x[1:(m-d)] # bug may here
-    return diff(ddmat(x, d - 1))
-  end
-end
-
-# 高斯消元法
-function LU_gauss(A)
-  n = size(A, 1)
-  T = typeof(A)
-  L = T(diagm(ones(n)))
-
-  for i = 1:n-1
-    r1 = A[i, :]
-    # U[i, :] = r1
-    for j = i+1:n
-      f = A[j, i] / A[i, i]
-      L[j, i] = f
-      A[j, :] .= A[j, :] .- (f * r1)
-      # 为啥要引入U，这样已经求解完成了
-      # L[:, 1] = A₁[j, :]
-    end
-  end
-  (; L, U=A)
-end
-
-# Doolittle's method
-# > 可避免修改矩阵A
-function LU_full(a)
-  n = size(a, 1)
-  u = tri_upper(:u, n)
-  l = tri_lower(:l, n)
-  for i = 1:n
-    l[i, i] = 1
-    for j = i:n
-      u[i, j] = a[i, j] - sum(l[i, 1:i-1] .* u[1:i-1, j])
-    end
-
-    for i2 = i+1:n
-      l[i2, i] = (a[i2, i] - sum(l[i2, 1:i-1] .* u[1:i-1, i])) / u[i, i]
-    end
-  end
-  l, u
-end
-
-# 带状矩阵的LU分解, Doolittle's method
-# - [x] 有效的减少循环次数
-# - [x] 压缩数据存储
-function LU_band(a; p=2, q=1)
-  n = size(a, 1)
-  u = tri_upper(:u, n)
-  # l = tri_lower(:l, n)
-  l = variables(:l, 1:p, 1:n)   # [i+k, i] -> [k, i]  ;  [i, j] -> [i-j, j]
-  u = variables(:l, 1:n, 1:q+1) # [i, i+k] -> [i, k+1];  [i, j] -> [i, j-i+1]
-
-  fill!(u, 0)
-  fill!(l, 0)
-  for i = 1:n
-    ## U矩阵同时压缩数据
-    for j = i:min(i + q, n)
-      u[i, j-i+1] = a[i, j]
-      for k = max(i - p, j - q, 1):min(i - 1, j - 1)
-        # 1 <= j - k <= q  ==> j - q <= k <= j - 1
-        # 1 <= i - k <= p  ==> i - p <= k <= i - 1
-        u[i, j-i+1] -= l[i-k, k] * u[k, j-k+1]
-      end
-    end
-
-    for i2 = i+1:min(i + p, n)
-      l[i2-i, i] = a[i2, i]
-      for k = max(i2 - p, i - q, 1):min(i - 1, i2 - 1)
-        # 1 <= i - k <= q  ==> i - q <= k <= i - 1
-        # 1 <= i2 - k <= p  ==> i2 - p <= k <= i2 - 1
-        l[i2-i, i] -= l[i2-k, k] * u[k, i-k+1]
-      end
-      l[i2-i, i] /= u[i, 1]
-    end
+function zip_U(U::AbstractMatrix, q)
+  U2 = variables(:U, 1:n, 1:q+1)
+  fill!(U2, 0)
+  for i = 1:n, j = i:min(n, q + i)
+    # 0 <= j-i <= q ==> i <= j <= q+i
+    U2[i, j-i+1] = U[i, j]
   end
-  l, u
+  U2
 end

src/Smooth/Whittaker/WHIT.jl

@@ -1,22 +1,33 @@
 using SparseArrays
 
-speye(n) = SparseArrays.sparse(I, n, n)
+function Base.diff(x::AbstractMatrix, d::Integer=1)
+  D = x[2:end, :] .- x[1:end-1, :]
+  d >= 2 ? diffn(D, d - 1) : D
+end
 
 function Base.diff(x::SparseMatrixCSC, d::Integer=1)
   D = x[2:end, :] .- x[1:end-1, :]
   d >= 2 ? diff(D, d - 1) : D
 end
 
+# ddmat(x::AbstractVector, d::Integer=2) = diff(spdiagm(x), d)
+# function Base.diff(x::SparseMatrixCSC, d::Integer=1)
+#   D = x[2:end, :] .- x[1:end-1, :]
+#   d >= 2 ? diff(D, d - 1) : D
+# end
+
+# 相比于原版进行了微调
 function ddmat(x::AbstractVector, d::Integer=2)
   n = length(x)
   if d == 0
-    return speye(n)
+    return sparse(I, n, n)
   else
-    return diff(ddmat(x, d - 1))
+    dx = ( x[d+1:n] - x[1:(n-d)] ) ./ d
+    V = spdiagm( 1 ./ dx)
+    return V * diff(ddmat(x, d - 1))
   end
 end
 
-
 function WHIT(y::AbstractVector, w::AbstractVector; kw...)
   WHIT(y, w, 1:length(y); kw...)
 end

test/test-smooth_whit.jl

@@ -1,5 +1,6 @@
 using Test
 using UnPack
+using VegCurveFit
 
 @testset "whittaker smoother" begin
   y = [5.0, 8, 9, 10, 12, 10, 15, 10, 9, 19, 19, 17, 13, 14, 18, 19, 18, 12, 18,
@@ -18,7 +19,7 @@ using UnPack
 
   lamb_cv = lambda_cv(y, w, is_plot=true)
   lamb_vcurve = lambda_vcurve(y, w, is_plot=true)
-  
+
   z1, cve_cv = whit2(y, w, lambda=lamb_cv)
   z2, cve_vcurve = whit2(y, w, lambda=lamb_vcurve)
   @test cve_cv < cve
@@ -28,23 +29,26 @@ end
 
 
 @testset "whit2" begin
-  d = deserialize("../data/Tumbarumba_EVI2")
+  d = deserialize("data/Tumbarumba_EVI2")
   @unpack y, t, w = d
 
-  z, cve1 = whit2_cv(y, w, lambda=2)
-  z, cve2 = whit2(y, w, lambda=2)
+  z1, cve1 = whit2_cv(y, w, lambda=2)
+  z2, cve2 = whit2(y, w, lambda=2)
   @test cve1 ≈ cve2
   @test_nowarn r = smooth_whit(y, w)
 
   # whit3
+  x = 1:length(y)
   z1, cve1 = whit3(y, w; lambda=2)
-  z2, cve2 = WHIT(y, w; lambda=2, d=3)
+  z2, cve2 = WHIT(y, w, x; lambda=2, d=3)
+
   @test cve1 ≈ cve2
   @test maximum(z1 - z2) <= 1e-10
 
   # whit2
   z1, cve1 = whit2(y, w; lambda=2)
-  z2, cve2 = WHIT(y, w; lambda=2, d=2)
+  z2, cve2 = WHIT(y, w, x; lambda=2, d=2)
+
   @test cve1 ≈ cve2
   @test maximum(z1 - z2) <= 1e-10
 end