From be05c3b510baac3d5216652fc8662f058806c8b8 Mon Sep 17 00:00:00 2001 From: PharmCat Date: Mon, 2 Jan 2023 03:50:23 +0300 Subject: [PATCH] optim --- src/gmat.jl | 47 ++++++++++++----------- src/rmat.jl | 107 ++++++++++++++++++++++------------------------------ 2 files changed, 71 insertions(+), 83 deletions(-) diff --git a/src/gmat.jl b/src/gmat.jl index 62df0880..4e2786ce 100644 --- a/src/gmat.jl +++ b/src/gmat.jl @@ -57,9 +57,10 @@ function gmat!(mx, θ, ::AR_) mx[i, i] = de end if s > 1 - for m = 1:s - 1 - @inbounds @simd for n = m + 1:s - mx[m, n] = de * θ[2] ^ (n - m) + @inbounds θ2 = θ[2] + for n = 2:s + @inbounds @simd for m = 1:n-1 + mx[m, n] = de * θ2 ^ (n - m) end end end @@ -72,14 +73,16 @@ function gmat!(mx, θ, ::ARH_) mx[m, m] = θ[m] end if s > 1 - for m = 1:s - 1 - @inbounds @simd for n = m + 1:s - mx[m, n] = mx[m, m] * mx[n, n] * θ[end] ^ (n - m) + θe = last(θ) + for n = 2:s + mxnn = mx[n, n] + @inbounds @simd for m = 1:n-1 + mx[m, n] = mx[m, m] * mxnn * θe ^ (n - m) end end end @inbounds @simd for m = 1:s - mx[m, m] = mx[m, m] * mx[m, m] + mx[m, m] *= mx[m, m] end mx end @@ -105,14 +108,14 @@ function gmat!(mx, θ, ::CSH_) end if s > 1 for n = 2:s - @inbounds mxnθe = mx[n, n] * θ[end] + @inbounds mxnθe = mx[n, n] * last(θ) @inbounds @simd for m = 1:n-1 mx[m, n] = mx[m, m] * mxnθe end end end @inbounds @simd for m = 1:s - mx[m, m] = mx[m, m] * mx[m, m] + mx[m, m] *= mx[m, m] end mx end @@ -125,9 +128,10 @@ function gmat!(mx, θ, ::ARMA_) mx[i, i] = de end if s > 1 - for m = 1:s - 1 - @inbounds @simd for n = m + 1:s - mx[m, n] = de * θ[2] * θ[3] ^ (n - m - 1) + deθ2 = de * θ[2] + for n = 2:s + @inbounds @simd for m = 1:n-1 + mx[m, n] = deθ2 * θ[3] ^ (n - m - 1) end end end @@ -141,8 +145,8 @@ function gmat!(mx, θ, ::TOEP_) mx[i, i] = de end if s > 1 - for m = 1:s - 1 - @inbounds @simd for n = m + 1:s + for n = 2:s + @inbounds @simd for m = 1:n-1 mx[m, n] = de * θ[n-m+1] end end @@ -171,14 +175,14 @@ function gmat!(mx, θ, ::TOEPH_) mx[m, m] = θ[m] end if s > 1 - for m = 1:s - 1 - @inbounds @simd for n = m + 1:s + for n = 2:s + @inbounds @simd for m = 1:n-1 mx[m, n] = mx[m, m] * mx[n, n] * θ[n-m+s] end end end @inbounds @simd for m = 1:s - mx[m, m] = mx[m, m] * mx[m, m] + mx[m, m] *= mx[m, m] end mx end @@ -196,7 +200,7 @@ function gmat!(mx, θ, ct::TOEPHP_) end end @inbounds @simd for m = 1:s - mx[m, m] = mx[m, m] * mx[m, m] + mx[m, m] *= mx[m, m] end mx end @@ -207,15 +211,14 @@ function gmat!(mx, θ, ::UN_) mx[m, m] = θ[m] end if s > 1 - for m = 1:s - 1 - @inbounds @simd for n = m + 1:s + for n = 2:s + @inbounds @simd for m = 1:n-1 mx[m, n] = mx[m, m] * mx[n, n] * θ[s+tpnum(m, n, s)] end end end @inbounds @simd for m = 1:s - v = mx[m, m] - mx[m, m] = v * v + mx[m, m] *= mx[m, m] end mx end diff --git a/src/rmat.jl b/src/rmat.jl index 852aba34..1966305a 100644 --- a/src/rmat.jl +++ b/src/rmat.jl @@ -3,22 +3,6 @@ ################################################################################ ################################################################################ -#= -function rmat_base_inc_b!(mx, θ, zrv, covstr) - rmat!(mx, θ, zrv, covstr.repeated.covtype.s) -end -=# -################################################################################ -################################################################################ -#= -function rmat_base_inc!(mx, θ, covstr, block, sblock) - zblock = view(covstr.rz, block, :) - @simd for i ∈ axes(sblock[end], 1) - @inbounds rmat_base_inc_b!(view(mx, sblock[end][i], sblock[end][i]), θ, view(zblock, sblock[end][i], :), covstr) - end - mx -end -=# @noinline function rmat_base_inc!(mx, θ, covstr, bi) en = covstr.rn + 1 block = covstr.vcovblock[bi] @@ -43,23 +27,23 @@ Base.@propagate_inbounds function rmat!(mx, θ, ::AbstractMatrix, ::SI_) end #DIAG function rmat!(mx, θ, rz, ::DIAG_) - #=@turbo=# for i ∈ axes(mx, 1) + for i ∈ axes(mx, 1) @inbounds @simd for c ∈ axes(θ, 1) - mx[i, i] += rz[i, c] * θ[c] * θ[c] + mx[i, i] += rz[i, c] * θ[c] ^ 2 end end mx end #AR function rmat!(mx, θ, ::AbstractMatrix, ::AR_) - rn = size(mx, 1) + s = size(mx, 1) de = θ[1] ^ 2 - @inbounds @simd for m = 1:rn + @inbounds @simd for m = 1:s mx[m, m] += de end - if rn > 1 - for m = 1:rn - 1 - @inbounds @simd for n = m + 1:rn + if s > 1 + for n = 2:s + @inbounds @simd for m = 1:n-1 mx[m, n] += de * θ[2] ^ (n - m) end end @@ -69,10 +53,10 @@ end #ARH function rmat!(mx, θ, rz, ::ARH_) vec = tmul_unsafe(rz, θ) - rn = size(mx, 1) - if rn > 1 - for m = 1:rn - 1 - @inbounds @simd for n = m + 1:rn + s = size(mx, 1) + if s > 1 + for n = 2:s + @inbounds @simd for m = 1:n-1 mx[m, n] += vec[m] * vec[n] * θ[end] ^ (n - m) end end @@ -84,15 +68,15 @@ function rmat!(mx, θ, rz, ::ARH_) end #CS function rmat!(mx, θ, ::AbstractMatrix, ::CS_) - rn = size(mx, 1) + s = size(mx, 1) θsq = θ[1]*θ[1] θsqp = θsq*θ[2] @inbounds @simd for i = 1:size(mx, 1) mx[i, i] += θsq end - if rn > 1 - for m = 1:rn - 1 - @inbounds @simd for n = m + 1:rn + if s > 1 + for n = 2:s + @inbounds @simd for m = 1:n-1 mx[m, n] += θsqp end end @@ -120,15 +104,16 @@ end ################################################################################ #ARMA function rmat!(mx, θ, ::AbstractMatrix, ::ARMA_) - rn = size(mx, 1) + s = size(mx, 1) de = θ[1] ^ 2 - @inbounds @simd for m = 1:rn + @inbounds @simd for m = 1:s mx[m, m] += de end - if rn > 1 - for m = 1:rn - 1 - @inbounds @simd for n = m + 1:rn - mx[m, n] += de * θ[2] * θ[3] ^ (n - m - 1) + if s > 1 + deθ2 = de * θ[2] + for n = 2:s + @inbounds @simd for m = 1:n-1 + mx[m, n] += deθ2 * θ[3] ^ (n - m - 1) end end end @@ -195,13 +180,13 @@ function rmat!(mx, θ, rz, ::SPEXP_) θe = θ[2] θe = abs(θe) < eps() ? sqrt(eps()) : abs(θe) #θe = abs(θ[2]) - rn = size(mx, 1) + s = size(mx, 1) @simd for i = 1:size(mx, 1) @inbounds mx[i, i] += σ² end - if rn > 1 - for m = 1:rn - 1 - @simd for n = m + 1:rn + if s > 1 + for m = 1:s - 1 + @simd for n = m + 1:s mx[m, n] += σ² * exp(-edistance(rz, m, n) / θe) end end @@ -213,13 +198,13 @@ end function rmat!(mx, θ, rz, ::SPPOW_) σ² = θ[1]^2 ρ = θ[2] - rn = size(mx, 1) + s = size(mx, 1) @simd for i = 1:size(mx, 1) @inbounds mx[i, i] += σ² end - if rn > 1 - for m = 1:rn - 1 - @simd for n = m + 1:rn + if s > 1 + for m = 1:s - 1 + @simd for n = m + 1:s ed = edistance(rz, m, n) mx[m, n] += ρ > 0 ? σ² * ρ^ed : σ² * ρ^ed * cos(π * ed) #mx[m, n] += ρ > 0 ? σ² * ρ^ed : σ² * abs(ρ)^ed * sign(ρ) @@ -237,13 +222,13 @@ function rmat!(mx, θ, rz, ::SPGAU_) θe = abs(θe) < eps() ? sqrt(eps()) : θe^2 #θe = θ[2] #θe = abs(θ[2]) - rn = size(mx, 1) + s = size(mx, 1) @simd for i = 1:size(mx, 1) @inbounds mx[i, i] += σ² end - if rn > 1 - for m = 1:rn - 1 - @simd for n = m + 1:rn + if s > 1 + for m = 1:s - 1 + @simd for n = m + 1:s mx[m, n] += σ² * exp(- (edistance(rz, m, n)^2) / θe) end end @@ -257,13 +242,13 @@ function rmat!(mx, θ, rz, ::SPEXPD_) σ²d = θ[1]^2 + σ² θe = θ[3] θe = iszero(θe) ? 1e-16 : abs(θe) - rn = size(mx, 1) + s = size(mx, 1) @simd for i = 1:size(mx, 1) @inbounds mx[i, i] += σ²d end - if rn > 1 - for m = 1:rn - 1 - @simd for n = m + 1:rn + if s > 1 + for m = 1:s - 1 + @simd for n = m + 1:s mx[m, n] += σ² * exp(-edistance(rz, m, n) / θe) end end @@ -275,13 +260,13 @@ function rmat!(mx, θ, rz, ::SPPOWD_) σ² = θ[2]^2 σ²d = θ[1]^2 + σ² ρ = θ[3] - rn = size(mx, 1) + s = size(mx, 1) @simd for i = 1:size(mx, 1) @inbounds mx[i, i] += σ²d end - if rn > 1 - for m = 1:rn - 1 - @simd for n = m + 1:rn + if s > 1 + for m = 1:s - 1 + @simd for n = m + 1:s mx[m, n] += σ² * ρ^edistance(rz, m, n) end end @@ -294,13 +279,13 @@ function rmat!(mx, θ, rz, ::SPGAUD_) σ²d = θ[1]^2 + σ² θe = θ[3] θe = iszero(θe) ? 1e-16 : θe^2 - rn = size(mx, 1) + s = size(mx, 1) @simd for i = 1:size(mx, 1) @inbounds mx[i, i] += σ²d end - if rn > 1 - for m = 1:rn - 1 - @simd for n = m + 1:rn + if s > 1 + for m = 1:s - 1 + @simd for n = m + 1:s mx[m, n] += σ² * exp(- (edistance(rz, m, n)^2) / θe) end end