refactor loops to eachindex for robustness

Project.toml

@@ -3,6 +3,7 @@ uuid = "3630a16b-0f2f-4d88-afbf-c7d59eccf553"
 version = "0.1.0"
 
 [deps]
+JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
 ManifoldsBase = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb"
@@ -13,6 +14,7 @@ RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 
 [compat]
+JuliaFormatter = "2.3.1"
 Manifolds = "0.11.20"
 ManifoldsBase = "2.3.5"
 Manopt = "0.5.37"

src/api/cpd.jl

@@ -49,14 +49,6 @@ Base.@kwdef mutable struct _CPDComponentTraceHistory
     ambient_velocity_top3_share_history::Vector{Float64} = Float64[]
     ambient_velocity_effective_components_history::Vector{Float64} = Float64[]
     ambient_velocity_argmax_component_history::Vector{Int} = Int[]
-    # Legacy aliases (coordinate / block-norm rgrad diagnostics).
-    rgrad_energy_history::Vector{Vector{Float64}} = Vector{Float64}[]
-    rgrad_share_history::Vector{Vector{Float64}} = Vector{Float64}[]
-    rgrad_top1_share_history::Vector{Float64} = Float64[]
-    rgrad_top2_share_history::Vector{Float64} = Float64[]
-    rgrad_top3_share_history::Vector{Float64} = Float64[]
-    rgrad_effective_components_history::Vector{Float64} = Float64[]
-    rgrad_argmax_component_history::Vector{Int} = Int[]
     rgrad_failed_count::Int = 0
 end
 
@@ -93,7 +85,7 @@ function _cpd_component_deltas(prev::CPDPoint, curr::CPDPoint)
     U_curr = factors(curr)
     r = length(λ_curr)
     deltas = Vector{Float64}(undef, r)
-    @inbounds for k = 1:r
+    @inbounds for k in eachindex(deltas)
         n_prev = _rankone_norm2(λ_prev, U_prev, k)
         n_curr = _rankone_norm2(λ_curr, U_curr, k)
         cross = _rankone_inner(λ_prev, U_prev, λ_curr, U_curr, k)
@@ -142,7 +134,7 @@ end
 function _cpd_coordinate_rgrad_energy_from_blocks(gλ, gU)
     r = length(gλ)
     energies = zeros(Float64, r)
-    @inbounds for k = 1:r
+    @inbounds for k in eachindex(energies)
         energies[k] += abs2(Float64(gλ[k]))
         for m in eachindex(gU)
             energies[k] += sum(abs2, @view gU[m][:, k])
@@ -191,7 +183,7 @@ end
 function _cpd_join_component_metric_energy(M::ProductManifold, p, g, k::Int, d::Int)
     base = (k - 1) * (d + 1)
     energy = 0.0
-    @inbounds for b = 1:(d+1)
+    @inbounds for b in eachindex(Base.OneTo(d + 1))
         energy += _cpd_product_block_metric_energy(M, p, g, base + b)
     end
     return energy
@@ -209,7 +201,7 @@ function _cpd_canonical_component_metric_energy(M::ProductManifold, p, g, k::Int
     pp = point_parts(p)
     gp = point_parts(g)
     energy = abs2(Float64(gp[1][k]))
-    @inbounds for m = 1:d
+    @inbounds for m in eachindex(Base.OneTo(d))
         mode_M = factors[m+1]
         energy += Float64(
             ManifoldsBase.inner(mode_M.manifolds[k], pp[m+1][k], gp[m+1][k], gp[m+1][k]),
@@ -232,11 +224,20 @@ function _cpd_component_metric_rgrad_energies(model::RankRCPDModel, p, g)
     if M isa ProductManifold
         nf = length(M.manifolds)
         if nf == r * (d + 1)
-            return [_cpd_join_component_metric_energy(M, p, g, k, d) for k = 1:r]
+            return [
+                _cpd_join_component_metric_energy(M, p, g, k, d) for
+                k in eachindex(Base.OneTo(r))
+            ]
         elseif nf == r
-            return [_cpd_native_component_metric_energy(M, p, g, k) for k = 1:r]
+            return [
+                _cpd_native_component_metric_energy(M, p, g, k) for
+                k in eachindex(Base.OneTo(r))
+            ]
         elseif nf == d + 1
-            return [_cpd_canonical_component_metric_energy(M, p, g, k, d) for k = 1:r]
+            return [
+                _cpd_canonical_component_metric_energy(M, p, g, k, d) for
+                k in eachindex(Base.OneTo(r))
+            ]
         end
     end
     return _cpd_component_coordinate_rgrad_energies(model, g)
@@ -314,7 +315,7 @@ function _cpd_component_ambient_velocities(model::RankRCPDModel, p, g)
     U = factors(q)
     r = model.r
     energies = Vector{Float64}(undef, r)
-    @inbounds for k = 1:r
+    @inbounds for k in eachindex(energies)
         δλk, δu = _decoded_tangent_for_component_k(model.geometry, λ̃, Ũ, gλ, gU, k)
         u_cols = [@view U[m][:, k] for m in eachindex(U)]
         energies[k] = _ambient_rankone_tangent_norm2(λ[k], u_cols, δλk, δu)
@@ -367,13 +368,6 @@ function _push_component_energy_trace!(
         push!(hist.coordinate_rgrad_top3_share_history, summary.top3)
         push!(hist.coordinate_rgrad_effective_components_history, summary.effective)
         push!(hist.coordinate_rgrad_argmax_component_history, summary.argmax_component)
-        push!(hist.rgrad_energy_history, energies_f)
-        push!(hist.rgrad_share_history, summary.shares)
-        push!(hist.rgrad_top1_share_history, summary.top1)
-        push!(hist.rgrad_top2_share_history, summary.top2)
-        push!(hist.rgrad_top3_share_history, summary.top3)
-        push!(hist.rgrad_effective_components_history, summary.effective)
-        push!(hist.rgrad_argmax_component_history, summary.argmax_component)
     elseif kind == :metric
         push!(hist.metric_rgrad_energy_history, energies_f)
         push!(hist.metric_rgrad_share_history, summary.shares)
@@ -492,34 +486,23 @@ function _cpd_component_trace_info(rec::_CPDComponentTraceRecorder)
         component_trace_ambient_velocity_top3_share_history = hist.ambient_velocity_top3_share_history,
         component_trace_ambient_velocity_effective_components_history = hist.ambient_velocity_effective_components_history,
         component_trace_ambient_velocity_argmax_component_history = hist.ambient_velocity_argmax_component_history,
-        component_trace_rgrad_energy_history = hist.rgrad_energy_history,
-        component_trace_rgrad_share_history = hist.rgrad_share_history,
-        component_trace_rgrad_top1_share_history = hist.rgrad_top1_share_history,
-        component_trace_rgrad_top2_share_history = hist.rgrad_top2_share_history,
-        component_trace_rgrad_top3_share_history = hist.rgrad_top3_share_history,
-        component_trace_rgrad_effective_components_history = hist.rgrad_effective_components_history,
-        component_trace_rgrad_argmax_component_history = hist.rgrad_argmax_component_history,
         component_trace_rgrad_failed_count = hist.rgrad_failed_count,
-        component_trace_rgrad_top1_share_final = _trace_history_final(
-            hist.rgrad_top1_share_history,
+        component_trace_coordinate_rgrad_top1_share_final = _trace_history_final(
+            hist.coordinate_rgrad_top1_share_history,
             NaN,
         ),
-        component_trace_rgrad_top2_share_final = _trace_history_final(
-            hist.rgrad_top2_share_history,
+        component_trace_coordinate_rgrad_top2_share_final = _trace_history_final(
+            hist.coordinate_rgrad_top2_share_history,
             NaN,
         ),
-        component_trace_rgrad_top3_share_final = _trace_history_final(
-            hist.rgrad_top3_share_history,
+        component_trace_coordinate_rgrad_top3_share_final = _trace_history_final(
+            hist.coordinate_rgrad_top3_share_history,
             NaN,
         ),
-        component_trace_rgrad_effective_components_final = _trace_history_final(
-            hist.rgrad_effective_components_history,
+        component_trace_coordinate_rgrad_effective_components_final = _trace_history_final(
+            hist.coordinate_rgrad_effective_components_history,
             NaN,
         ),
-        component_trace_rgrad_argmax_component_final = _trace_history_final(
-            hist.rgrad_argmax_component_history,
-            0,
-        ),
         component_trace_coordinate_rgrad_argmax_component_final = _trace_history_final(
             hist.coordinate_rgrad_argmax_component_history,
             0,

src/btd/core/btd_grad.jl

@@ -6,7 +6,7 @@ function _btd_core_grad_block(backend::BTDBackend, parts, b::Int)
 
     grad_core = copy(Gb)
     grad_core .-= _tucker_project_target(backend, pb, backend.target)
-    @inbounds for c = 1:backend.r
+    @inbounds for c in eachindex(parts)
         c == b && continue
         grad_core .+= _tucker_cross_core(backend, pb, parts[c])
     end
@@ -19,12 +19,13 @@ function _btd_factor_grad_block(backend::BTDBackend, parts, b::Int)
     T = eltype(Gb)
     N = length(Ub)
     grad_factors = Vector{Matrix{T}}(undef, N)
-    core_unfolds = [unfold_mode(Gb, m) for m = 1:N]
-    core_gram_factors = [core_unfolds[m] * transpose(core_unfolds[m]) for m = 1:N]
+    core_unfolds = [unfold_mode(Gb, m) for m in eachindex(Ub)]
+    core_gram_factors =
+        [core_unfolds[m] * transpose(core_unfolds[m]) for m in eachindex(Ub)]
     oneT = one(T)
     zeroT = zero(T)
 
-    @inbounds for m = 1:N
+    @inbounds for m in eachindex(Ub)
         Gbm = core_unfolds[m]
         GbmT = transpose(Gbm)
         nmode = size(Ub[m], 1)
@@ -37,7 +38,7 @@ function _btd_factor_grad_block(backend::BTDBackend, parts, b::Int)
 
         mul!(grad_m, Ub[m], core_gram_factors[m], oneT, oneT)
 
-        for c = 1:backend.r
+        for c in eachindex(parts)
             c == b && continue
             _, Uc = _tucker_data(parts[c])
             cross_unfold =

src/btd/core/inner_prod.jl

@@ -186,7 +186,7 @@ function _tucker_tucker_inner(p::Manifolds.TuckerPoint, q::Manifolds.TuckerPoint
         throw(DimensionMismatch("Tucker core/factor count mismatch for second point."))
     N == length(Uq) ||
         throw(DimensionMismatch("Tucker mode count mismatch between points."))
-    @inbounds for k = 1:N
+    @inbounds for k in eachindex(Up)
         size(Up[k], 2) == size(Gp, k) ||
             throw(DimensionMismatch("Tucker factor dimension mismatch for mode $k."))
         size(Uq[k], 2) == size(Gq, k) ||
@@ -196,7 +196,7 @@ function _tucker_tucker_inner(p::Manifolds.TuckerPoint, q::Manifolds.TuckerPoint
     end
 
     Ht = Gq
-    @inbounds for k = 1:N
+    @inbounds for k in eachindex(Up)
         Ht = mode_n_product(Ht, Up[k]' * Uq[k], k)
     end
     return sum(Gp .* Ht)
@@ -226,7 +226,7 @@ function _target_tucker_inner(A::AbstractArray, p::Manifolds.TuckerPoint)
     G, U = _tucker_data(p)
     N = length(U)
     ndims(G) == N || throw(DimensionMismatch("Tucker core/factor count mismatch."))
-    @inbounds for k = 1:N
+    @inbounds for k in eachindex(U)
         size(U[k], 2) == size(G, k) ||
             throw(DimensionMismatch("Tucker factor dimension mismatch for mode $k."))
         size(U[k], 1) == size(A, k) ||

src/btd/model.jl

@@ -70,7 +70,7 @@ function _btd_sequential_tucker_init(model::JoinModel{<:AbstractFloat,<:BTDBacke
     # same full-tensor fit into every block.
     residual = copy(backend.target)
     parts = Vector{Manifolds.TuckerPoint{eltype(backend.target)}}(undef, backend.r)
-    for k = 1:backend.r
+    for k in eachindex(parts)
         pk = _manifold_init(backend.manifolds[k], residual, init)
         parts[k] = pk
         _subtract_ambient_tensor!(
@@ -85,7 +85,7 @@ end
 
 function _btd_block_ranks_by_mode(backend::BTDBackend{T,N}) where {T,N}
     ranks = Vector{NTuple{N,Int}}(undef, backend.r)
-    for b = 1:backend.r
+    for b in eachindex(backend.manifolds)
         M = backend.manifolds[b]
         M isa Manifolds.Tucker || throw(
             ArgumentError(
@@ -131,7 +131,7 @@ end
 
 function _btd_project_core(A::AbstractArray{T,N}, factors) where {T<:AbstractFloat,N}
     core = A
-    for mode = 1:N
+    for mode in eachindex(factors)
         core = mode_n_product(core, factors[mode]', mode)
     end
     return core
@@ -145,14 +145,14 @@ function _btd_hosvd_split_candidate(
     candidate::Int,
 ) where {T<:AbstractFloat,N}
     columns_by_mode = ntuple(N) do mode
-        ranks_m = [ranks_by_block[b][mode] for b = 1:backend.r]
+        ranks_m = [ranks_by_block[b][mode] for b in eachindex(ranks_by_block)]
         _btd_split_columns(rng, size(subspaces[mode], 2), ranks_m, candidate)
     end
 
     residual = copy(backend.target)
     # Candidate blocks are always Tucker points here
     parts = Vector{Manifolds.TuckerPoint{T}}(undef, backend.r)
-    for b = 1:backend.r
+    for b in eachindex(parts)
         factors =
             ntuple(mode -> Matrix(@view subspaces[mode][:, columns_by_mode[mode][b]]), N)
         core = _btd_project_core(residual, factors)
@@ -234,7 +234,7 @@ function initial_point(
     best_cost = Inf
     split_candidate = 1
 
-    for c = 1:init.candidates
+    for c in eachindex(Base.OneTo(init.candidates))
         p_candidate = if init.include_sequential && c == 1
             _btd_sequential_tucker_init(model, :sthosvd)
         else

src/core/initialization.jl

@@ -198,7 +198,7 @@ end
 
 function random_unit_matrix(n::Int, r::Int, ::Type{T} = Float64) where {T<:AbstractFloat}
     U = zeros(T, n, r)
-    for k = 1:r
+    for k in axes(U, 2)
         @views U[:, k] .= random_unit_vector(n, T)
     end
     return U
@@ -209,7 +209,7 @@ function _tucker_diag(core::AbstractArray{T}, r::Int) where {T<:AbstractFloat}
     out = zeros(T, r)
     maxr = minimum(size(core))
     rr = min(r, maxr)
-    for k = 1:rr
+    for k in eachindex(Base.OneTo(rr))
         idx = ntuple(_ -> k, N)
         out[k] = core[idx...]
     end

src/core/pack_points.jl

@@ -36,16 +36,16 @@ function pack_rankr_native(
 ) where {T<:AbstractFloat}
     d = length(U)
     length(λ) == r || throw(DimensionMismatch("length(λ)=$(length(λ)) must equal r=$r"))
-    for m = 1:d
+    for m in eachindex(U)
         size(U[m], 2) == r ||
             throw(DimensionMismatch("U[$m] has $(size(U[m],2)) columns, expected r=$r"))
     end
 
     comps = Vector{Vector{Vector{T}}}(undef, r)
-    @inbounds for k = 1:r
+    @inbounds for k in eachindex(λ)
         λk = λ[k]
         Uk = Vector{Vector{T}}(undef, d)
-        for m = 1:d
+        for m in eachindex(U)
             u = Vector{T}(@view U[m][:, k])
             λk = _normalize_column_into_lambda!(u, λk)
             Uk[m] = u
@@ -91,16 +91,16 @@ function pack_rankr_canonical_tuple(
 ) where {T<:AbstractFloat}
     d = length(U)
     length(λ) == r || throw(DimensionMismatch("length(λ)=$(length(λ)) must equal r=$r"))
-    for m = 1:d
+    for m in eachindex(U)
         size(U[m], 2) == r ||
             throw(DimensionMismatch("U[$m] has $(size(U[m],2)) columns, expected r=$r"))
     end
 
     λn = Vector{T}(undef, r)
-    mode_cols = [Vector{Vector{T}}(undef, r) for _ = 1:d]
-    @inbounds for k = 1:r
+    mode_cols = [Vector{Vector{T}}(undef, r) for _ in eachindex(U)]
+    @inbounds for k in eachindex(λ)
         λk = λ[k]
-        for m = 1:d
+        for m in eachindex(U)
             u = Vector{T}(@view U[m][:, k])
             λk = _normalize_column_into_lambda!(u, λk)
             mode_cols[m][k] = u
@@ -118,16 +118,16 @@ function pack_rankr_join_tuple(
 ) where {T<:AbstractFloat}
     d = length(U)
     length(λ) == r || throw(DimensionMismatch("length(λ)=$(length(λ)) must equal r=$r"))
-    for m = 1:d
+    for m in eachindex(U)
         size(U[m], 2) == r ||
             throw(DimensionMismatch("U[$m] has $(size(U[m],2)) columns, expected r=$r"))
     end
     parts = Vector{Vector{T}}(undef, r * (d + 1))
     idx = 1
-    @inbounds for k = 1:r
+    @inbounds for k in eachindex(λ)
         parts[idx] = T[λ[k]]
         idx += 1
-        for m = 1:d
+        for m in eachindex(U)
             parts[idx] = Vector{T}(@view U[m][:, k])
             idx += 1
         end