diff --git a/Project.toml b/Project.toml
index bd4414c..09e5246 100644
--- a/Project.toml
+++ b/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"
diff --git a/src/api/cpd.jl b/src/api/cpd.jl
index 12f5094..3aa8fa3 100644
--- a/src/api/cpd.jl
+++ b/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,
diff --git a/src/btd/core/btd_grad.jl b/src/btd/core/btd_grad.jl
index 1bfafa5..4fb741f 100644
--- a/src/btd/core/btd_grad.jl
+++ b/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 =
diff --git a/src/btd/core/inner_prod.jl b/src/btd/core/inner_prod.jl
index bed16b7..30cf757 100644
--- a/src/btd/core/inner_prod.jl
+++ b/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) ||
diff --git a/src/btd/model.jl b/src/btd/model.jl
index 509905d..0c36e88 100644
--- a/src/btd/model.jl
+++ b/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
diff --git a/src/core/initialization.jl b/src/core/initialization.jl
index e295da5..7bcd7b4 100644
--- a/src/core/initialization.jl
+++ b/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
diff --git a/src/core/pack_points.jl b/src/core/pack_points.jl
index c51475d..cb8e919 100644
--- a/src/core/pack_points.jl
+++ b/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
diff --git a/src/core/tensor_ops.jl b/src/core/tensor_ops.jl
index 5ab2f89..1746eaa 100644
--- a/src/core/tensor_ops.jl
+++ b/src/core/tensor_ops.jl
@@ -20,7 +20,7 @@ using TensorOperations
     perm = Vector{Int}(undef, N)
     perm[1] = mode
     j = 2
-    @inbounds for m = 1:N
+    @inbounds for m in eachindex(perm)
         m == mode && continue
         perm[j] = m
         j += 1
@@ -32,7 +32,7 @@ end
     N = length(U)
     mats = Vector{eltype(U)}(undef, N - 1)
     j = 1
-    @inbounds for m = N:-1:1
+    @inbounds for m in Iterators.reverse(eachindex(U))
         m == mode && continue
         mats[j] = U[m]
         j += 1
@@ -95,7 +95,7 @@ end
     U::Vector{<:AbstractVector{T}},
 ) where {T<:AbstractFloat,N}
     prod_val = one(T)
-    @inbounds for m = 1:N
+    @inbounds for m in eachindex(U)
         prod_val *= U[m][I[m]]
     end
     return prod_val
@@ -104,7 +104,7 @@ end
 @inline function _rank1_entry_product_parts(I::CartesianIndex{N}, parts) where {N}
     T = eltype(parts[2])
     prod_val = one(T)
-    @inbounds for m = 1:N
+    @inbounds for m in eachindex(Base.OneTo(N))
         prod_val *= parts[m+1][I[m]]
     end
     return prod_val
@@ -116,7 +116,7 @@ end
     mode::Int,
 ) where {T<:AbstractFloat,N}
     prod_val = one(T)
-    @inbounds for m = 1:N
+    @inbounds for m in eachindex(U)
         m == mode && continue
         prod_val *= U[m][I[m]]
     end
@@ -130,7 +130,7 @@ end
 ) where {N}
     T = eltype(parts[2])
     prod_val = one(T)
-    @inbounds for m = 1:N
+    @inbounds for m in eachindex(Base.OneTo(N))
         m == mode && continue
         prod_val *= parts[m+1][I[m]]
     end
@@ -144,7 +144,7 @@ end
     k::Int,
 ) where {T<:AbstractFloat,N}
     prod_val = one(T)
-    @inbounds for m = 1:N
+    @inbounds for m in eachindex(U)
         m == mode && continue
         prod_val *= U[m][I[m], k]
     end
@@ -316,7 +316,7 @@ function rank1_mode_contract_column!(
     k::Int,
 ) where {T<:AbstractFloat,N}
     fill!(out, zero(T))
-    @inbounds for m = 1:N
+    @inbounds for m in eachindex(U)
         size(U[m], 1) == size(A, m) ||
             throw(DimensionMismatch("mode $m factor row count mismatch"))
     end
@@ -399,8 +399,8 @@ end
 function build_cross_matrix(components::Vector{RankOneTensor{T}}) where {T<:AbstractFloat}
     r = length(components)
     cross_mat = zeros(T, r, r)
-    for k = 1:r
-        for l = 1:r
+    for k in eachindex(components)
+        for l in eachindex(components)
             cross_mat[k, l] = cross_component(components[k], components[l])
         end
     end
@@ -412,8 +412,8 @@ function build_cross_matrix_unit(
 ) where {T<:AbstractFloat}
     r = length(components)
     cross_mat = zeros(T, r, r)
-    for k = 1:r
-        for l = 1:r
+    for k in eachindex(components)
+        for l in eachindex(components)
             cross_mat[k, l] = prod(
                 dot(components[k].vectors[m], components[l].vectors[m]) for
                 m in eachindex(components[1].vectors)
@@ -538,7 +538,10 @@ function factors_from_components(
     isempty(components) && return Vector{Matrix{T}}()
     r = length(components)
     N = length(components[1].vectors)
-    return [hcat((components[k].vectors[m] for k = 1:r)...) for m = 1:N]
+    return [
+        hcat((components[k].vectors[m] for k in eachindex(components))...) for
+        m in eachindex(components[1].vectors)
+    ]
 end
 
 function components_from_factors(
@@ -547,5 +550,8 @@ function components_from_factors(
 ) where {T<:AbstractFloat}
     r = length(λ)
     N = length(U)
-    return [RankOneTensor(λ[k], [Vector(@view U[m][:, k]) for m = 1:N]) for k = 1:r]
+    return [
+        RankOneTensor(λ[k], [Vector(@view U[m][:, k]) for m in eachindex(U)]) for
+        k in eachindex(λ)
+    ]
 end
diff --git a/src/core/unpack_points.jl b/src/core/unpack_points.jl
index 8828e86..ad417d0 100644
--- a/src/core/unpack_points.jl
+++ b/src/core/unpack_points.jl
@@ -9,9 +9,9 @@ function unpack_rankr_native(p, dims::NTuple{N,Int}, r::Int) where {N}
     d = length(dims)
     proto = parts[1][2]
     λ = similar(proto, T, r)
-    U = [similar(proto, T, dims[m], r) for m = 1:d]
+    U = [similar(proto, T, dims[m], r) for m in eachindex(dims)]
 
-    @inbounds for k = 1:r
+    @inbounds for k in eachindex(λ)
         comp = parts[k]
         length(comp) == d + 1 || throw(
             DimensionMismatch(
@@ -19,7 +19,7 @@ function unpack_rankr_native(p, dims::NTuple{N,Int}, r::Int) where {N}
             ),
         )
         λ[k] = comp[1][1]
-        for m = 1:d
+        for m in eachindex(dims)
             @views U[m][:, k] .= comp[m+1]
         end
     end
@@ -36,11 +36,11 @@ function unpack_rankr_canonical(p, dims::NTuple{N,Int}, r::Int) where {N}
     λ = parts[1]
     T = eltype(λ)
     U = Vector{Matrix{T}}(undef, d)
-    @inbounds for m = 1:d
+    @inbounds for m in eachindex(dims)
         mode_m = parts[m+1]
         proto = mode_m[1]
         Um = similar(proto, T, dims[m], r)
-        for k = 1:r
+        for k in eachindex(λ)
             uk = mode_m[k]
             length(uk) == dims[m] || throw(
                 DimensionMismatch(
@@ -59,12 +59,12 @@ function unpack_rankr_join(p, dims::NTuple{N,Int}, r::Int) where {N}
     T = eltype(parts[1])
     proto = parts[2]
     λ = similar(proto, T, r)
-    U = [similar(proto, T, dims[m], r) for m = 1:N]
+    U = [similar(proto, T, dims[m], r) for m in eachindex(dims)]
     idx = 1
-    @inbounds for k = 1:r
+    @inbounds for k in eachindex(λ)
         λ[k] = parts[idx][1]
         idx += 1
-        for m = 1:N
+        for m in eachindex(dims)
             @views U[m][:, k] .= parts[idx]
             idx += 1
         end
@@ -167,12 +167,12 @@ function unpack_point_rankr(
     d = length(dims)
     block = 1 + sum(dims)
     λ = similar(p, T, r)
-    U = [similar(p, T, dims[m], r) for m = 1:d]
-    for k = 1:r
+    U = [similar(p, T, dims[m], r) for m in eachindex(dims)]
+    for k in eachindex(λ)
         b = (k - 1) * block
         λ[k] = p[b+1]
         idx = b + 2
-        for m = 1:d
+        for m in eachindex(dims)
             n = dims[m]
             src = @view p[idx:(idx+n-1)]
             @views U[m][:, k] .= src
@@ -187,5 +187,9 @@ function pack_point_rank1_to_vector(λ̃::T, U::Vector{Vector{T}}) where {T}
 end
 
 function pack_point_rankr_to_vector(λ::Vector{T}, U::Vector{Matrix{T}}, r::Int) where {T}
-    return vcat((vcat(λ[k], ((@view U[m][:, k]) for m in eachindex(U))...) for k = 1:r)...)
+    return vcat(
+        (
+            vcat(λ[k], ((@view U[m][:, k]) for m in eachindex(U))...) for k in eachindex(λ)
+        )...,
+    )
 end
diff --git a/src/cpd/core/cp_cost.jl b/src/cpd/core/cp_cost.jl
index c99547b..4d4d32a 100644
--- a/src/cpd/core/cp_cost.jl
+++ b/src/cpd/core/cp_cost.jl
@@ -161,7 +161,7 @@ function _inner_from_mttkrp_first_mode(
 ) where {T<:AbstractFloat}
     r = size(M1, 2)
     inner = Vector{T}(undef, r)
-    @inbounds for k = 1:r
+    @inbounds for k in eachindex(inner)
         inner[k] = dot(@view(U[1][:, k]), @view(M1[:, k]))
     end
     return inner
@@ -181,8 +181,8 @@ function _CPRankrEvalCache(
     dims::NTuple{N,Int},
     r::Int,
 ) where {T<:AbstractFloat,N}
-    contracts = [Matrix{T}(undef, dims[m], r) for m = 1:N]
-    grams = [Matrix{T}(undef, r, r) for _ = 1:N]
+    contracts = [Matrix{T}(undef, dims[m], r) for m in eachindex(dims)]
+    grams = [Matrix{T}(undef, r, r) for _ in eachindex(dims)]
     return _CPRankrEvalCache{T}(
         nothing,
         false,
@@ -226,13 +226,13 @@ function _cp_rankr_refresh_cache!(
         return cache
     end
 
-    @inbounds for m = 1:N
+    @inbounds for m in eachindex(U)
         copyto!(cache.contracts[m], mttkrp(A, U, m; method))
         mul!(cache.grams[m], transpose(U[m]), U[m])
     end
 
     fill!(cache.cross_mat, one(T))
-    @inbounds for m = 1:N
+    @inbounds for m in eachindex(U)
         cache.cross_mat .*= cache.grams[m]
     end
 
@@ -265,9 +265,9 @@ function _rankr_gradU_from_terms(
     r = length(λ)
     gradU = Vector{Matrix{T}}(undef, Nmodes)
     λouter = λ * transpose(λ)
-    @inbounds for m = 1:Nmodes
+    @inbounds for m in eachindex(U)
         prod_except = ones(T, r, r)
-        for j = 1:Nmodes
+        for j in eachindex(U)
             j == m && continue
             prod_except .*= grams[j]
         end
@@ -321,7 +321,7 @@ function egrad_secant_rankr(
         Nmodes = length(U)
 
         contracts = Vector{Matrix{T}}(undef, Nmodes)
-        for m = 1:Nmodes
+        for m in eachindex(U)
             contracts[m] = mttkrp(A, U, m; method = :auto)
         end
 
@@ -333,9 +333,9 @@ function egrad_secant_rankr(
         gradU = _rankr_gradU_from_terms(U, λ, contracts, grams)
 
         if scale_by_lambda
-            for k = 1:r
+            for k in eachindex(λ)
                 λ_abs = max(abs(λ[k]), lambda_eps_T)
-                for m = 1:Nmodes
+                for m in eachindex(gradU)
                     gradU[m][:, k] ./= λ_abs
                 end
             end
@@ -357,7 +357,7 @@ end
         throw(DimensionMismatch("expected $(N+1) tuple parts, got $(length(parts))"))
     λ = parts[1]
     length(λ) == r || throw(DimensionMismatch("expected λ length $r, got $(length(λ))"))
-    @inbounds for m = 1:N
+    @inbounds for m in eachindex(Ubuf)
         mode_m = parts[m+1]
         length(mode_m) == r ||
             throw(DimensionMismatch("mode $m has $(length(mode_m)) vectors, expected $r"))
@@ -367,7 +367,7 @@ end
         else
             Um = Ubuf[m]
         end
-        for k = 1:r
+        for k in eachindex(Base.OneTo(r))
             uk = mode_m[k]
             length(uk) == dims[m] || throw(
                 DimensionMismatch(
@@ -429,9 +429,9 @@ function egrad_rankr_canonical(
         gradU = _rankr_gradU_from_terms(Ubuf, λ, cache.contracts, cache.grams)
 
         if scale_by_lambda
-            for k = 1:r
+            for k in eachindex(λ)
                 λ_abs = max(abs(λ[k]), lambda_eps_T)
-                for m = 1:Nmodes
+                for m in eachindex(gradU)
                     gradU[m][:, k] ./= λ_abs
                 end
             end
@@ -491,17 +491,17 @@ function egrad_rankr_native(
         gradU = _rankr_gradU_from_terms(U, λ, cache.contracts, cache.grams)
 
         if scale_by_lambda
-            for k = 1:r
+            for k in eachindex(λ)
                 λ_abs = max(abs(λ[k]), lambda_eps_T)
-                for m = 1:Nmodes
+                for m in eachindex(gradU)
                     gradU[m][:, k] ./= λ_abs
                 end
             end
         end
 
         grad_parts = Vector{Vector{Vector{T}}}(undef, r)
-        @inbounds for k = 1:r
-            grad_Uk = [Vector(@view gradU[m][:, k]) for m = 1:Nmodes]
+        @inbounds for k in eachindex(grad_parts)
+            grad_Uk = [Vector(@view gradU[m][:, k]) for m in eachindex(gradU)]
             grad_parts[k] = pack_tangent_rank1_segre(grad_λ[k], grad_Uk)
         end
 
diff --git a/src/cpd/core/cp_normalization.jl b/src/cpd/core/cp_normalization.jl
index 0c0f61e..3c4e5c9 100644
--- a/src/cpd/core/cp_normalization.jl
+++ b/src/cpd/core/cp_normalization.jl
@@ -52,7 +52,7 @@ function normalize_components!(
     isempty(factors) && return factors
     r = length(lambda)
     d = length(factors)
-    @inbounds for m = 1:d
+    @inbounds for m in eachindex(factors)
         size(factors[m], 2) == r || throw(
             DimensionMismatch("factor $m has $(size(factors[m], 2)) columns, expected $r"),
         )
@@ -84,9 +84,9 @@ function _normalize_components_policy!(
 ) where {T<:AbstractFloat}
     r = length(lambda)
     d = length(factors)
-    @inbounds for k = 1:r
+    @inbounds for k in eachindex(lambda)
         lambda[k] = max(lambda[k], zero(T))
-        for m = 1:d
+        for m in eachindex(factors)
             col = @view factors[m][:, k]
             col .= max.(col, zero(T))
         end
@@ -101,9 +101,9 @@ function _normalize_components_policy!(
 ) where {T<:AbstractFloat}
     r = length(lambda)
     d = length(factors)
-    @inbounds for k = 1:r
+    @inbounds for k in eachindex(lambda)
         scale = lambda[k]
-        for m = 1:d
+        for m in eachindex(factors)
             col = @view factors[m][:, k]
             scale = _normalize_column_into_lambda!(col, scale)
         end
diff --git a/src/cpd/core/cp_points.jl b/src/cpd/core/cp_points.jl
index 1815cd2..9ef2db8 100644
--- a/src/cpd/core/cp_points.jl
+++ b/src/cpd/core/cp_points.jl
@@ -91,7 +91,7 @@ function normalize_rank1_segre_point(p, dims::NTuple{N,Int}) where {N}
     λ = isfinite(λpart0[1]) ? λpart0[1] : zero(T)
     flip_first = λ < 0
     out[1] = T[abs(λ)]
-    @inbounds for m = 1:N
+    @inbounds for m in eachindex(dims)
         u0 = _unwrap_part(parts[m+1])
         u0 isa AbstractVector || throw(
             DimensionMismatch("Segre mode $m part must be a vector, got $(typeof(u0))"),
@@ -119,7 +119,7 @@ function normalize_rankr_native_point(p, dims::NTuple{N,Int}, r::Int) where {N}
     parts = parts_tuple(p)
     length(parts) == r ||
         throw(DimensionMismatch("expected $r Segre components, got $(length(parts))"))
-    comps = [normalize_rank1_segre_point(parts[k], dims) for k = 1:r]
+    comps = [normalize_rank1_segre_point(parts[k], dims) for k in eachindex(parts)]
     return (comps...,)
 end
 
@@ -150,7 +150,7 @@ function normalize_rankr_canonical_point(p, dims::NTuple{N,Int}, r::Int) where {
                 DimensionMismatch("mode $m has $(length(mode_m)) vectors, expected $r"),
             )
             cols = Vector{Vector{T}}(undef, r)
-            for k = 1:r
+            for k in eachindex(cols)
                 uk0 = _unwrap_part(mode_m[k])
                 uk0 isa AbstractVector || throw(
                     DimensionMismatch(
@@ -192,7 +192,7 @@ function normalize_rankr_join_point(p, dims::NTuple{N,Int}, r::Int) where {N}
 
     out = Vector{Vector{T}}(undef, expected)
     idx = 1
-    @inbounds for k = 1:r
+    @inbounds for k in eachindex(Base.OneTo(r))
         λ0 = _unwrap_part(parts[idx])
         λ0 isa AbstractVector ||
             throw(DimensionMismatch("join λ[$k] must be a vector, got $(typeof(λ0))"))
@@ -200,7 +200,7 @@ function normalize_rankr_join_point(p, dims::NTuple{N,Int}, r::Int) where {N}
             throw(DimensionMismatch("join λ[$k] must have length 1, got $(length(λ0))"))
         out[idx] = Vector{T}(λ0)
         idx += 1
-        for m = 1:N
+        for m in eachindex(dims)
             u0 = _unwrap_part(parts[idx])
             u0 isa AbstractVector || throw(
                 DimensionMismatch(
diff --git a/src/cpd/core/cpd_init.jl b/src/cpd/core/cpd_init.jl
index b6964e7..9843d7e 100644
--- a/src/cpd/core/cpd_init.jl
+++ b/src/cpd/core/cpd_init.jl
@@ -36,15 +36,15 @@ function _cp_core_diag_init(core::AbstractArray{T,N}, r::Int) where {T<:Abstract
     core_dims = size(core)
     λ0 = _tucker_diag(core, r)
     U0 = Vector{Matrix{T}}(undef, N)
-    for m = 1:N
+    for m in eachindex(core_dims)
         rm = core_dims[m]
         Um = fill!(similar(core, rm, r), zero(T))
         n_eye = min(rm, r)
-        for k = 1:n_eye
+        for k in eachindex(Base.OneTo(n_eye))
             Um[k, k] = one(T)
         end
         if rm > 0 && r > n_eye
-            Um[:, n_eye+1:r] .= random_unit_matrix(rm, r - n_eye, T)
+            Um[:, (n_eye+1):r] .= random_unit_matrix(rm, r - n_eye, T)
         end
         U0[m] = Um
     end
@@ -77,7 +77,7 @@ function _cp_least_squares_lambda(
 ) where {T<:AbstractFloat,N}
     dlen = length(A)
     Z = similar(A, T, dlen, r)
-    @inbounds for k = 1:r
+    @inbounds for k in axes(Z, 2)
         vecs = [collect(@view U0[m][:, k]) for m in eachindex(U0)]
         Z[:, k] .= vec(reconstruct_cp_rank1(one(T), vecs))
     end
@@ -90,7 +90,8 @@ function init_cpd_factors(
     init::Symbol = :random,
 ) where {T<:AbstractFloat,N}
     dims = size(A)
-    init == :random && return (ones(T, r), [random_unit_matrix(dims[m], r, T) for m = 1:N])
+    init == :random &&
+        return (ones(T, r), [random_unit_matrix(dims[m], r, T) for m in eachindex(dims)])
     init in (:hosvd, :tucker_diag, :tucker) || throw(
         ArgumentError("Unknown init=$init. Use :random, :hosvd, :tucker_diag, or :tucker."),
     )
@@ -101,7 +102,7 @@ function init_cpd_factors(
         hcat(
             factors[m],
             random_unit_matrix(size(factors[m], 1), r - size(factors[m], 2), T),
-        ) : factors[m] for m = 1:N
+        ) : factors[m] for m in eachindex(factors)
     ]
     λ0 = if init == :tucker_diag
         _tucker_diag(core, r)
@@ -118,12 +119,12 @@ function init_cp_rank1(
     init::Symbol = :random,
 ) where {T<:AbstractFloat,N}
     dims = size(A)
-    init == :random && return [random_unit_vector(dims[m], T) for m = 1:N]
+    init == :random && return [random_unit_vector(dims[m], T) for m in eachindex(dims)]
     init in (:hosvd, :tucker_diag, :tucker) || throw(
         ArgumentError("Unknown init=$init. Use :random, :hosvd, :tucker_diag, or :tucker."),
     )
     _, factors = tucker_hosvd(A, ntuple(_ -> 1, N))
-    return [vec(factors[m][:, 1]) for m = 1:N]
+    return [vec(factors[m][:, 1]) for m in eachindex(factors)]
 end
 
 function cp_init_tucker(
diff --git a/src/cpd/core/mttkrp.jl b/src/cpd/core/mttkrp.jl
index df8d2b2..45ee5c2 100644
--- a/src/cpd/core/mttkrp.jl
+++ b/src/cpd/core/mttkrp.jl
@@ -1,5 +1,5 @@
 # cpd/core/mttkrp.jl — CPD-specific MTTKRP kernels and dispatch
-# Improvement of resolving mttkrp bottleneck still in progress
+# Improvement of resolving mttkrp bottleneck
 export mttkrp, mttkrp!, khatri_rao, khatri_rao!
 @inline _mttkrp_needs_kr_workspace(method::Symbol) = method == :khatri_rao
 @inline _mttkrp_needs_tmp_workspace(method::Symbol) = method in (:direct3, :direct4)
@@ -63,16 +63,16 @@ function khatri_rao(mats::AbstractVector{<:AbstractMatrix{T}}) where {T<:Abstrac
         size(m, 2) == r || throw(ArgumentError("khatri_rao: column counts must match"))
     end
     out = copy(mats[1])
-    for i = 2:length(mats)
+    for i in Iterators.drop(eachindex(mats), 1)
         A = mats[i]
         new = similar(out, T, size(out, 1) * size(A, 1), r)
         rows_out = size(out, 1)
         rows_A = size(A, 1)
-        @inbounds for k = 1:r
+        @inbounds for k in axes(out, 2)
             idx = 1
-            for io = 1:rows_out
+            for io in axes(out, 1)
                 scale = out[io, k]
-                for ia = 1:rows_A
+                for ia in axes(A, 1)
                     new[idx, k] = scale * A[ia, k]
                     idx += 1
                 end
@@ -111,15 +111,15 @@ function khatri_rao!(
     copyto!(view(out, 1:rows_cur, :), mats[1])
     src = out
     dst = work
-    for i = 2:length(mats)
+    for i in Iterators.drop(eachindex(mats), 1)
         A = mats[i]
         rows_A = size(A, 1)
         new_rows = rows_cur * rows_A
-        @inbounds for k = 1:r
+        @inbounds for k in axes(out, 2)
             idx = 1
-            for io = 1:rows_cur
+            for io in eachindex(Base.OneTo(rows_cur))
                 scale = src[io, k]
-                for ia = 1:rows_A
+                for ia in axes(A, 1)
                     dst[idx, k] = scale * A[ia, k]
                     idx += 1
                 end
@@ -209,19 +209,24 @@ function _mttkrp_direct3!(
     mode::Int,
 ) where {T<:AbstractFloat}
     I, J, K = size(A)
+    k_axis = axes(U[3], 1)
     fill!(out, zero(T))
     if mode == 1
-        @inbounds for k = 1:K
+        @inbounds for kk in eachindex(k_axis)
+            k = k_axis[kk]
             mul!(tmp, @view(A[:, :, k]), U[2])
             _accumulate_scaled_columns!(out, tmp, @view(U[3][k, :]))
         end
     elseif mode == 2
-        @inbounds for k = 1:K
+        @inbounds for kk in eachindex(k_axis)
+            k = k_axis[kk]
             mul!(tmp, transpose(@view(A[:, :, k])), U[1])
             _accumulate_scaled_columns!(out, tmp, @view(U[3][k, :]))
         end
     elseif mode == 3
-        @inbounds for j = 1:J
+        j_axis = axes(U[2], 1)
+        @inbounds for jj in eachindex(j_axis)
+            j = j_axis[jj]
             mul!(tmp, transpose(reshape(@view(A[:, j, :]), I, K)), U[1])
             _accumulate_scaled_columns!(out, tmp, @view(U[2][j, :]))
         end
@@ -249,32 +254,43 @@ function _mttkrp_direct4!(
     U::AbstractVector{<:AbstractMatrix{T}},
     mode::Int,
 ) where {T<:AbstractFloat}
-    I, J, K, L = size(A)
+    I, _, K, L = size(A)
+    j_axis = axes(U[2], 1)
+    k_axis = axes(U[3], 1)
+    l_axis = axes(U[4], 1)
     fill!(out, zero(T))
     if mode == 1
-        @inbounds for l = 1:L
-            for k = 1:K
+        @inbounds for ll in eachindex(l_axis)
+            l = l_axis[ll]
+            for kk in eachindex(k_axis)
+                k = k_axis[kk]
                 mul!(tmp, @view(A[:, :, k, l]), U[2])
                 _accumulate_scaled_columns!(out, tmp, @view(U[3][k, :]), @view(U[4][l, :]))
             end
         end
     elseif mode == 2
-        @inbounds for l = 1:L
-            for k = 1:K
+        @inbounds for ll in eachindex(l_axis)
+            l = l_axis[ll]
+            for kk in eachindex(k_axis)
+                k = k_axis[kk]
                 mul!(tmp, transpose(@view(A[:, :, k, l])), U[1])
                 _accumulate_scaled_columns!(out, tmp, @view(U[3][k, :]), @view(U[4][l, :]))
             end
         end
     elseif mode == 3
-        @inbounds for l = 1:L
-            for j = 1:J
+        @inbounds for ll in eachindex(l_axis)
+            l = l_axis[ll]
+            for jj in eachindex(j_axis)
+                j = j_axis[jj]
                 mul!(tmp, transpose(reshape(@view(A[:, j, :, l]), I, K)), U[1])
                 _accumulate_scaled_columns!(out, tmp, @view(U[2][j, :]), @view(U[4][l, :]))
             end
         end
     elseif mode == 4
-        @inbounds for k = 1:K
-            for j = 1:J
+        @inbounds for kk in eachindex(k_axis)
+            k = k_axis[kk]
+            for jj in eachindex(j_axis)
+                j = j_axis[jj]
                 mul!(tmp, transpose(reshape(@view(A[:, j, k, :]), I, L)), U[1])
                 _accumulate_scaled_columns!(out, tmp, @view(U[2][j, :]), @view(U[3][k, :]))
             end
@@ -294,7 +310,7 @@ function _mttkrp_contract(
     r = size(U[1], 2)
     out = similar(U[1], T, dims[mode], r)
     fill!(out, zero(T))
-    for k = 1:r
+    for k in eachindex(axes(out, 2))
         @views out[:, k] .= rank1_mode_contract_column(A, U, mode, k)
     end
     return out
@@ -307,7 +323,7 @@ function _mttkrp_contract!(
     mode::Int,
 ) where {T<:AbstractFloat,N}
     r = size(U[1], 2)
-    @inbounds for q = 1:r
+    @inbounds for q in eachindex(axes(out, 2))
         @views out[:, q] .= rank1_mode_contract_column(A, U, mode, q)
     end
     return out
@@ -320,7 +336,7 @@ function _mttkrp_contract(
 ) where {T<:AbstractFloat}
     r = size(U[1], 2)
     out = similar(U[1], T, size(A, mode), r)
-    @inbounds for q = 1:r
+    @inbounds for q in eachindex(axes(out, 2))
         @views out[:, q] .= rank1_mode_contract_column(A, U, mode, q)
     end
     return out
@@ -333,7 +349,7 @@ function _mttkrp_contract(
 ) where {T<:AbstractFloat}
     r = size(U[1], 2)
     out = similar(U[1], T, size(A, mode), r)
-    @inbounds for q = 1:r
+    @inbounds for q in eachindex(axes(out, 2))
         @views out[:, q] .= rank1_mode_contract_column(A, U, mode, q)
     end
     return out
@@ -363,7 +379,7 @@ function mttkrp(
         throw(DimensionMismatch("mttkrp: expected $N factor matrices, got $(length(U))"))
 
     r = size(U[1], 2)
-    for m = 1:N
+    for m in eachindex(U)
         size(U[m], 1) == dims[m] || throw(
             DimensionMismatch(
                 "mttkrp: U[$m] has $(size(U[m], 1)) rows, expected $(dims[m])",
@@ -414,7 +430,7 @@ function mttkrp!(
         throw(DimensionMismatch("mttkrp: expected $N factor matrices, got $(length(U))"))
 
     r = size(U[1], 2)
-    @inbounds for m = 1:N
+    @inbounds for m in eachindex(U)
         size(U[m], 1) == dims[m] || throw(
             DimensionMismatch(
                 "mttkrp: U[$m] has $(size(U[m], 1)) rows, expected $(dims[m])",
diff --git a/src/cpd/core/reconstruct.jl b/src/cpd/core/reconstruct.jl
index 8dcaff2..2e31f4d 100644
--- a/src/cpd/core/reconstruct.jl
+++ b/src/cpd/core/reconstruct.jl
@@ -18,8 +18,8 @@ function reconstruct_cp_rank1(
     isempty(U) && throw(ArgumentError("reconstruct_cp_rank1: empty factor vectors"))
     N = length(U)
     # A[i1,...,iN] = λ * u1[i1] * u2[i2] * ... * uN[iN] — outer product scaled by λ
-    tensors = [vec(U[m]) for m = 1:N]
-    return T(λ) .* ncon(tensors, [[-m] for m = 1:N])
+    tensors = [vec(U[m]) for m in eachindex(U)]
+    return T(λ) .* ncon(tensors, [[-m] for m in eachindex(U)])
 end
 
 """Rank-r CP from (λ, U). Batched broadcast + sum over r."""
@@ -27,7 +27,7 @@ function reconstruct_cpd_rankr(λ::Vector{T}, U::Vector{Matrix{T}}) where {T<:Ab
     (isempty(U) || isempty(λ)) &&
         throw(ArgumentError("reconstruct_cpd_rankr: empty factors or λ"))
     r, N = length(λ), length(U)
-    for m = 1:N
+    for m in eachindex(U)
         size(U[m], 2) == r || throw(
             DimensionMismatch(
                 "reconstruct_cpd_rankr: U[$m] has $(size(U[m],2)) cols, expected r=$r",
@@ -56,7 +56,7 @@ function reconstruct_cpd_rankr(
                 "reconstruct_cpd_rankr: component $k has $(length(components[k].vectors)) modes",
             ),
         )
-        for m = 1:N
+        for m in eachindex(components[1].vectors)
             length(components[k].vectors[m]) == length(components[1].vectors[m]) || throw(
                 DimensionMismatch(
                     "reconstruct_cpd_rankr: component $k mode $m incompatible",
@@ -65,7 +65,10 @@ function reconstruct_cpd_rankr(
         end
     end
     λ = [c.λ for c in components]
-    U = [hcat((c.vectors[m] for c in components)...) for m = 1:N]
+    U = [
+        hcat((c.vectors[m] for c in components)...) for
+        m in eachindex(components[1].vectors)
+    ]
     return reconstruct_cpd_rankr(λ, U)
 end
 
diff --git a/src/cpd/model/rank1.jl b/src/cpd/model/rank1.jl
index 8903818..9af24f6 100644
--- a/src/cpd/model/rank1.jl
+++ b/src/cpd/model/rank1.jl
@@ -38,7 +38,7 @@ function Rank1CPDModel(
         factors[1] =
             use_softplus_metric ? SoftplusEuclidean(1; ε = pullback_eps) :
             (use_pullback_metric ? SqEuclidean(1; ε = pullback_eps) : Euclidean(1))
-        @inbounds for m = 1:N
+        @inbounds for m in eachindex(dims)
             factors[m+1] =
                 use_softplus_metric ? SoftplusEuclidean(dims[m]; ε = pullback_eps) :
                 (
diff --git a/src/cpd/model/rankr.jl b/src/cpd/model/rankr.jl
index 8df2904..e684912 100644
--- a/src/cpd/model/rankr.jl
+++ b/src/cpd/model/rankr.jl
@@ -61,12 +61,12 @@ function RankRCPDModel(
         n_factors = r * (Nd + 1)
         manifolds = Vector{AbstractManifold}(undef, n_factors)
         idx = 1
-        @inbounds for k = 1:r
+        @inbounds for k in eachindex(Base.OneTo(r))
             manifolds[idx] =
                 use_pullback_softplus ? SoftplusEuclidean(1; ε = pullback_eps) :
                 (use_pullback_sq ? SqEuclidean(1; ε = pullback_eps) : Euclidean(1))
             idx += 1
-            for m = 1:Nd
+            for m in eachindex(dims)
                 manifolds[idx] =
                     use_pullback_softplus ?
                     SoftplusEuclidean(dims[m]; ε = pullback_eps) :
@@ -185,7 +185,7 @@ function model_cost_egrad_functions(model::RankRCPDModel{T,N}) where {T,N}
             grad_λ = grad_lambda_cp(λ, cache.inner, cache.cross_mat)
             gradU = _rankr_gradU_from_terms(U, λ, cache.contracts, cache.grams)
             if model.scale_by_lambda
-                for k = 1:model.r
+                for k in eachindex(λ)
                     λ_abs = max(abs(λ[k]), model.lambda_eps)
                     for m in eachindex(U)
                         gradU[m][:, k] ./= λ_abs
@@ -193,7 +193,7 @@ function model_cost_egrad_functions(model::RankRCPDModel{T,N}) where {T,N}
                 end
             end
             grad_parts = Vector{Vector{Vector{T}}}(undef, model.r)
-            @inbounds for k = 1:model.r
+            @inbounds for k in eachindex(grad_parts)
                 grad_Uk = [Vector(@view gradU[m][:, k]) for m in eachindex(U)]
                 grad_parts[k] = pack_tangent_rank1_segre(grad_λ[k], grad_Uk)
             end
@@ -218,9 +218,9 @@ function model_cost_egrad_functions(model::RankRCPDModel{T,N}) where {T,N}
             grad_λ = grad_lambda_cp(λ, cache.inner, cache.cross_mat)
             gradU = _rankr_gradU_from_terms(Ubuf, λ, cache.contracts, cache.grams)
             if model.scale_by_lambda
-                for k = 1:model.r
+                for k in eachindex(λ)
                     λ_abs = max(abs(λ[k]), model.lambda_eps)
-                    for m = 1:Nmodes
+                    for m in eachindex(gradU)
                         gradU[m][:, k] ./= λ_abs
                     end
                 end
@@ -294,7 +294,7 @@ cp_als_data(model::RankRCPDModel) = (model.A, model.r)
 @inline function _segre_component_tensorvec(comp)
     parts = parts_tuple(comp)
     λ = parts[1][1]
-    U = [parts[m+1] for m = 1:(length(parts)-1)]
+    U = [parts[m+1] for m in eachindex(Base.OneTo(length(parts) - 1))]
     return vec(reconstruct_cp_rank1(λ, U))
 end
 
@@ -307,12 +307,12 @@ function _segre_tangent_tensorvec(comp, Xcomp)
 
     λ = pparts[1][1]
     ν = xparts[1][1]
-    U = [pparts[m+1] for m = 1:d]
-    Udot = [xparts[m+1] for m = 1:d]
+    U = [pparts[m+1] for m in eachindex(Base.OneTo(d))]
+    Udot = [xparts[m+1] for m in eachindex(Base.OneTo(d))]
 
     v = vec(reconstruct_cp_rank1(ν, U))
-    @inbounds for m = 1:d
-        Um = [j == m ? Udot[j] : U[j] for j = 1:d]
+    @inbounds for m in eachindex(Base.OneTo(d))
+        Um = [j == m ? Udot[j] : U[j] for j in eachindex(Base.OneTo(d))]
         v .+= vec(reconstruct_cp_rank1(λ, Um))
     end
     return v
@@ -384,7 +384,7 @@ function rgrad(model::RankRCPDModel{T,N}, p) where {T,N} # Riemannian gradient f
     r = model.r
 
     contracts = Vector{Matrix{T}}(undef, Nmodes)
-    for m = 1:Nmodes
+    for m in eachindex(U)
         contracts[m] = mttkrp(model.A, U, m; method = :auto)
     end
 
@@ -394,16 +394,16 @@ function rgrad(model::RankRCPDModel{T,N}, p) where {T,N} # Riemannian gradient f
     grad_λ = grad_lambda_cp(λ, inner, cross_mat)
 
     gradU = _rankr_gradU_from_terms(U, λ, contracts, grams)
-    for m = 1:Nmodes
+    for m in eachindex(U)
         Gm = gradU[m]
         if model.scale_by_lambda
-            for k = 1:r
+            for k in eachindex(λ)
                 Gm[:, k] ./= max(abs(λ[k]), model.lambda_eps)
             end
         end
 
         # Tangent projection on each sphere factor (u_{m,k}^T g_{m,k} = 0).
-        for k = 1:r
+        for k in eachindex(λ)
             uk = @view U[m][:, k] # view the k-th column of the m-th mode
             gk = @view Gm[:, k] # view the k-th column of the gradient for the m-th mode
             Gm[:, k] .-= dot(uk, gk) .* uk
diff --git a/src/join/join_backend.jl b/src/join/join_backend.jl
index e8d88d7..1939994 100644
--- a/src/join/join_backend.jl
+++ b/src/join/join_backend.jl
@@ -205,7 +205,8 @@ function _sum_backend_parts(
     tflat = vec(tgt)
     tgt_len = length(tgt)
     # One ambient buffer per component lets reconstruction reuse storage across iterations.
-    component_bufs = [_join_vector_workspace_like(tgt, tgt_len) for _ = 1:r]
+    component_bufs =
+        [_join_vector_workspace_like(tgt, tgt_len) for _ in eachindex(manifolds)]
     work_rec = _join_vector_workspace_like(tgt, tgt_len)
     work_residual = _join_vector_workspace_like(tgt, tgt_len)
 
@@ -418,12 +419,12 @@ function extract_components(
     N = length(backend.target_shape)
     point_type = Union{}
     manifold_type = Union{}
-    @inbounds for k = 1:backend.r
+    @inbounds for k in eachindex(parts)
         point_type = typejoin(point_type, typeof(parts[k]))
         manifold_type = typejoin(manifold_type, typeof(backend.manifolds[k]))
     end
     comps = Vector{DecompositionComponent{T,N,point_type,manifold_type}}(undef, backend.r)
-    @inbounds for k = 1:backend.r
+    @inbounds for k in eachindex(comps)
         # Components keep only point/manifold metadata and reconstruct derived tensors on demand.
         comps[k] = DecompositionComponent{T,N,point_type,manifold_type}(
             parts[k],
@@ -529,7 +530,7 @@ function _join_reconstruct!(out::AbstractArray, backend::Union{JoinBackend,BTDBa
 
     fill!(out, zero(eltype(out)))
 
-    @inbounds for k = 1:r
+    @inbounds for k in eachindex(parts)
         # Reconstruct each component into its preallocated workspace.
         _ambient_vector!(bufs[k], manifolds[k], parts[k])
 
diff --git a/src/join/join_init.jl b/src/join/join_init.jl
index 10a83b3..ba35c11 100644
--- a/src/join/join_init.jl
+++ b/src/join/join_init.jl
@@ -102,10 +102,10 @@ function _tucker_all_except_mode_products(
 ) where {T<:AbstractFloat,N}
     N == 1 && return Array{T,N}[Array(core)]
     products = Vector{Array{T,N}}(undef, N)
-    @inbounds for m = 1:N
+    @inbounds for m in eachindex(products)
         # Precompute all "except mode m" contractions once so factor gradients can reuse them.
         B = core
-        for j = 1:N
+        for j in eachindex(products)
             j == m && continue
             B = mode_n_product(B, factors[j], j)
         end
@@ -126,15 +126,15 @@ function _tucker_egrad(M::Manifolds.Tucker, p, R)
     N = length(factors)
 
     grad_core = copy(R)
-    for m = 1:N
+    for m in eachindex(factors)
         grad_core = mode_n_product(grad_core, factors[m]', m)
     end
 
-    residual_unfolds = [unfold_mode(R, m) for m = 1:N]
+    residual_unfolds = [unfold_mode(R, m) for m in eachindex(factors)]
     # Reuse these intermediates across all factor-gradient blocks.
     all_except_mode = _tucker_all_except_mode_products(core, factors)
     grad_factors = Vector{Matrix{eltype(core)}}(undef, N)
-    for m = 1:N
+    for m in eachindex(factors)
         Bm = unfold_mode(all_except_mode[m], m)
         grad_factors[m] = residual_unfolds[m] * transpose(Bm)
     end
diff --git a/src/manifolds/gradients.jl b/src/manifolds/gradients.jl
index 6eefd42..d81665c 100644
--- a/src/manifolds/gradients.jl
+++ b/src/manifolds/gradients.jl
@@ -70,7 +70,7 @@ function egrad_to_rgrad(M::Manifolds.Segre, p, egrad)
     out = Vector{Vector{T}}(undef, d + 1)
     out[1] = _grad_vec(gparts[1], T)
     length(out[1]) == 1 || throw(DimensionMismatch("Segre λ-gradient must have length 1."))
-    @inbounds for m = 1:d
+    @inbounds for m in eachindex(Base.OneTo(d))
         u = _grad_vec(parts[m+1], T)
         gm = _grad_vec(gparts[m+1], T)
         length(u) == length(gm) ||
@@ -141,7 +141,7 @@ function _embedded_euclidean_basis_fallback(
     coeff = zeros(T, d)
     e_j = zeros(T, d)
     u = similar(egrad, T, length(egrad))
-    @inbounds for j = 1:d
+    @inbounds for j in eachindex(Base.OneTo(d))
         fill!(e_j, zero(T))
         e_j[j] = one(T)
         ξj = ManifoldsBase.get_vector(M, p, e_j, basis)
diff --git a/src/manifolds/join.jl b/src/manifolds/join.jl
index f67f03d..65c8e6c 100644
--- a/src/manifolds/join.jl
+++ b/src/manifolds/join.jl
@@ -41,12 +41,16 @@ end
 
 function _join_product(base::Manifolds.Segre, r::Int)
     factors = _segre_flat_factors(base)
-    return ProductManifold([deepcopy(f) for _ = 1:r for f in factors]...)
+    return ProductManifold(
+        [deepcopy(f) for _ in eachindex(Base.OneTo(r)) for f in factors]...,
+    )
 end
 
 function _join_product(base::ProductManifold, r::Int)
     factors = base.manifolds
-    return ProductManifold([deepcopy(f) for _ = 1:r for f in factors]...)
+    return ProductManifold(
+        [deepcopy(f) for _ in eachindex(Base.OneTo(r)) for f in factors]...,
+    )
 end
 
 function _join_product(base::AbstractManifold, r::Int)
diff --git a/src/solvers/btd_als.jl b/src/solvers/btd_als.jl
index 2852167..1dc34c7 100644
--- a/src/solvers/btd_als.jl
+++ b/src/solvers/btd_als.jl
@@ -25,7 +25,12 @@ function _btd_tucker_point_to_result(p::Manifolds.TuckerPoint{T}) where {T<:Abst
     N = length(factors)
     d = ndims(core)
     d == N || throw(DimensionMismatch("_btd_tucker_point_to_result: core ndims=$d, N=$N"))
-    return TuckerResult{T,N}(core, collect(factors), collect(1:N), [T[] for _ = 1:N])
+    return TuckerResult{T,N}(
+        core,
+        collect(factors),
+        collect(eachindex(factors)),
+        [T[] for _ in eachindex(factors)],
+    )
 end
 
 function _btd_block_fit_tucker(
@@ -99,10 +104,10 @@ function fit_btd_als(
         parts0 = point_parts(p_start)
         _check_parts_len(parts0, backend.r, "BTD ALS init")
 
-        points = [parts0[k] for k = 1:backend.r]
-        block_tensors = [_btd_block_tensor(points[k]) for k = 1:backend.r]
+        points = [parts0[k] for k in eachindex(parts0)]
+        block_tensors = [_btd_block_tensor(points[k]) for k in eachindex(points)]
         residual = copy(A)
-        @inbounds for k = 1:backend.r
+        @inbounds for k in eachindex(points)
             residual .-= block_tensors[k]
         end
 
@@ -121,7 +126,7 @@ function fit_btd_als(
             ) : NoMethodProgress()
 
         for iter = 1:maxiter
-            @inbounds for b = 1:backend.r
+            @inbounds for b in eachindex(points)
                 residual .+= block_tensors[b]
 
                 ranks_b = _btd_block_ranks(backend, b)
diff --git a/src/solvers/cp_als.jl b/src/solvers/cp_als.jl
index 3135f84..e2bf0ad 100644
--- a/src/solvers/cp_als.jl
+++ b/src/solvers/cp_als.jl
@@ -62,27 +62,28 @@ function CPALSWorkspace(
     r::Int;
     mttkrp_method::Symbol = :auto,
 ) where {T<:AbstractFloat,N}
-    Gs = [_cp_als_matrix_workspace_like(A, r, r) for _ = 1:N]
+    Gs = [_cp_als_matrix_workspace_like(A, r, r) for _ in eachindex(dims)]
     V = _cp_als_matrix_workspace_like(A, r, r)
-    transposed_work = [_cp_als_matrix_workspace_like(A, r, dims[n]) for n = 1:N]
-    denom_work = [_cp_als_matrix_workspace_like(A, dims[n], r) for n = 1:N]
-    mttkrp_bufs = [_cp_als_matrix_workspace_like(A, dims[n], r) for n = 1:N]
+    transposed_work =
+        [_cp_als_matrix_workspace_like(A, r, dims[n]) for n in eachindex(dims)]
+    denom_work = [_cp_als_matrix_workspace_like(A, dims[n], r) for n in eachindex(dims)]
+    mttkrp_bufs = [_cp_als_matrix_workspace_like(A, dims[n], r) for n in eachindex(dims)]
     total_dim_prod = prod(dims)
     resolved_mttkrp_methods =
-        [_mttkrp_resolve_method(mttkrp_method, dims, r, n) for n = 1:N]
+        [_mttkrp_resolve_method(mttkrp_method, dims, r, n) for n in eachindex(dims)]
     mttkrp_tmp_work = Any[
         _mttkrp_needs_tmp_workspace(resolved_mttkrp_methods[n]) ?
-        _cp_als_matrix_workspace_like(A, dims[n], r) : nothing for n = 1:N
+        _cp_als_matrix_workspace_like(A, dims[n], r) : nothing for n in eachindex(dims)
     ]
     mttkrp_kr_work = Any[
         _mttkrp_needs_kr_workspace(resolved_mttkrp_methods[n]) ?
         _cp_als_matrix_workspace_like(A, div(total_dim_prod, dims[n]), r) : nothing for
-        n = 1:N
+        n in eachindex(dims)
     ]
     mttkrp_kr_work2 = Any[
         _mttkrp_needs_kr_workspace(resolved_mttkrp_methods[n]) ?
         _cp_als_matrix_workspace_like(A, div(total_dim_prod, dims[n]), r) : nothing for
-        n = 1:N
+        n in eachindex(dims)
     ]
     cross_buf = _cp_als_matrix_workspace_like(A, r, r)
     return CPALSWorkspace(
@@ -189,7 +190,7 @@ function _cp_als_stats(
         copyto!(cross, Gs[1])
     end
 
-    @inbounds for m = 2:length(Gs)
+    @inbounds for m in Iterators.drop(eachindex(Gs), 1)
         cross .*= Gs[m]
     end
     normX2 = zero(T)
@@ -287,7 +288,7 @@ function fit_cp_als(
 
     workspace = CPALSWorkspace(A, dims, r; mttkrp_method = mttkrp_method)
     Gs = workspace.Gs
-    @inbounds for n = 1:N
+    @inbounds for n in eachindex(U)
         _update_G!(Gs[n], U[n])
     end
     V = workspace.V
@@ -307,7 +308,7 @@ function fit_cp_als(
     for iter = 1:maxiter
         pg_sq = zero(T)
         u_sq = zero(T)
-        for n = 1:N
+        for n in eachindex(U)
             _hadamard_G_except!(V, Gs, n)
             M_mttkrp = mttkrp!(
                 mttkrp_bufs[n],
diff --git a/src/solvers/nncp_updates.jl b/src/solvers/nncp_updates.jl
index b497ed4..ffdf964 100644
--- a/src/solvers/nncp_updates.jl
+++ b/src/solvers/nncp_updates.jl
@@ -104,7 +104,7 @@ function _nncp_nnls_row_update!(
 ) where {T<:AbstractFloat}
     floor = sqrt(eps(T))
     mul!(work, V, x)
-    @inbounds for _ = 1:max_cd_sweeps
+    @inbounds for _ in eachindex(Base.OneTo(max_cd_sweeps))
         max_delta = zero(T)
         max_x = maximum(x)
         for k in eachindex(x)
@@ -135,7 +135,7 @@ function _nncp_nnls_row_update!(
 ) where {T<:AbstractFloat}
     floor = sqrt(eps(T))
     d = max.(diag(V), floor)
-    @inbounds for _ = 1:max_cd_sweeps
+    @inbounds for _ in eachindex(Base.OneTo(max_cd_sweeps))
         mul!(work, V, x)
         x_new = max.(x .- (work .- g) ./ d, floor)
         max_delta = maximum(abs.(x_new .- x))
@@ -177,7 +177,7 @@ function _nncp_nnls_mode_update!(
     floor = sqrt(eps(T))
     d = max.(diag(V), floor)
     d_row = reshape(d, 1, :)
-    @inbounds for _ = 1:max_cd_sweeps
+    @inbounds for _ in eachindex(Base.OneTo(max_cd_sweeps))
         mul!(work, U, V)
         U_new = max.(U .- (work .- M_mttkrp) ./ d_row, floor)
         max_delta = maximum(abs.(U_new .- U))
@@ -278,7 +278,7 @@ function _projected_grad_norm_nonnegative!(
 ) where {T<:AbstractFloat,N}
     sq = zero(T)
     u_sq = zero(T)
-    for n = 1:N
+    for n in eachindex(U)
         _hadamard_G_except!(V, grams, n)
         M_mttkrp = mttkrp!(
             mttkrp_bufs[n],
diff --git a/src/solvers/rals.jl b/src/solvers/rals.jl
index 4cd9861..1792019 100644
--- a/src/solvers/rals.jl
+++ b/src/solvers/rals.jl
@@ -31,7 +31,7 @@ function fit_cp_rals(
     else
         init_sym = _builtin_initializer_symbol(init)
         if init_sym == :random
-            (ones(T, r), [Matrix(qr(randn(T, dims[n], r)).Q) for n = 1:N])
+            (ones(T, r), [Matrix(qr(randn(T, dims[n], r)).Q) for n in eachindex(dims)])
         elseif init_sym in (:hosvd, :tucker, :tucker_diag)
             init_cpd_factors(A, r; init = init_sym)
         else
@@ -43,7 +43,7 @@ function fit_cp_rals(
     A_work = A
     if mix
         Q_mix = [Matrix(qr(randn(T, d, d)).Q) for d in dims]
-        for n = 1:N
+        for n in eachindex(dims)
             A_work = mode_n_product(A_work, Q_mix[n], n)
         end
     end
@@ -51,7 +51,9 @@ function fit_cp_rals(
     fit_old = 0.0
     converged = false
     iter_final = maxiter
-    others_indices = [CartesianIndices(Tuple(dims[setdiff(1:N, n)])) for n = 1:N]
+    others_indices = [
+        CartesianIndices(Tuple(dims[setdiff(eachindex(dims), n)])) for n in eachindex(dims)
+    ]
     progress =
         maxiter > 0 ?
         (
@@ -71,8 +73,8 @@ function fit_cp_rals(
         ) : NoMethodProgress()
 
     for iter = 1:maxiter
-        for n = 1:N
-            others = setdiff(1:N, n)
+        for n in eachindex(dims)
+            others = setdiff(eachindex(dims), n)
             sample_idx = rand(1:length(others_indices[n]), samples)
             coords = others_indices[n][sample_idx]
 
@@ -87,7 +89,7 @@ function fit_cp_rals(
             # We need to extract A[i, :, k, ...] where (i, k, ...) comes from coords
             # Construct indices for A
             X_s = zeros(T, samples, dims[n])
-            for s = 1:samples
+            for s in eachindex(coords)
                 idx = ntuple(d -> d == n ? Colon() : coords[s][findfirst(==(d), others)], N)
                 X_s[s, :] = A_work[idx...]
             end
@@ -136,7 +138,7 @@ function fit_cp_rals(
 
     # Un-mix factors if needed
     if mix
-        for n = 1:N
+        for n in eachindex(dims)
             # U_orig = Q * U_mixed
             U[n] = Q_mix[n] * U[n]
         end
diff --git a/src/tucker/hooi.jl b/src/tucker/hooi.jl
index 9307b79..555d519 100644
--- a/src/tucker/hooi.jl
+++ b/src/tucker/hooi.jl
@@ -19,7 +19,7 @@ function hooi(
 
     factors0, singular_vals = _hooi_initial_factors(A, ranks, init)
 
-    factors = [copy(factors0[m]) for m = 1:d]
+    factors = [copy(factors0[m]) for m in eachindex(factors0)]
     prev_rel_error = T(Inf)
     converged = false
     iter_final = maxiter
@@ -35,7 +35,7 @@ function hooi(
     for iter = 1:maxiter
         factors_prev = copy(factors)
         prefix = A
-        for k = 1:d
+        for k in eachindex(factors)
             Y = prefix
             for j = (k+1):d
                 Y = mode_n_product(Y, factors_prev[j]', j)
@@ -69,7 +69,7 @@ function hooi(
 
     if maxiter == 0
         S = copy(A)
-        for k = 1:d
+        for k in eachindex(factors)
             S = mode_n_product(S, factors[k]', k)
         end
     end
@@ -102,7 +102,7 @@ function _hooi_initial_factors(
             "hooi: TuckerResult.core has size $(size(init.core)), expected core size $ranks",
         ),
     )
-    @inbounds for m = 1:N
+    @inbounds for m in eachindex(ranks)
         size(init.factors[m], 1) == dims[m] || throw(
             DimensionMismatch(
                 "hooi: TuckerResult factor $m has $(size(init.factors[m], 1)) rows, expected $(dims[m])",
@@ -114,7 +114,7 @@ function _hooi_initial_factors(
             ),
         )
     end
-    return init.factors, [copy(init.singular_values[m]) for m = 1:N]
+    return init.factors, [copy(init.singular_values[m]) for m in eachindex(ranks)]
 end
 
 function _hooi_initial_factors(
diff --git a/src/tucker/hosvd.jl b/src/tucker/hosvd.jl
index 3bf72b9..dd64a9a 100644
--- a/src/tucker/hosvd.jl
+++ b/src/tucker/hosvd.jl
@@ -20,14 +20,14 @@ the original tensor onto these factor spaces.
 function tucker_hosvd(A::AbstractArray{T}, ranks::NTuple{N,Int}) where {T<:AbstractFloat,N}
     dims = size(A)
     factors = Vector{Matrix{T}}(undef, N)
-    for mode = 1:N
+    for mode in eachindex(factors)
         A_mode = unfold_mode(A, mode)
         U, _, _ = svd(A_mode)
         r = min(ranks[mode], size(U, 2))
         factors[mode] = Matrix(@view U[:, 1:r])
     end
     core = A
-    for mode = 1:N
+    for mode in eachindex(factors)
         core = mode_n_product(core, factors[mode]', mode)
     end
     return core, factors
@@ -65,7 +65,7 @@ function reconstruct_tucker!(
             "reconstruct_tucker!: got $(length(factors)) factors for an order-$N core.",
         ),
     )
-    @inbounds for mode = 1:N
+    @inbounds for mode in eachindex(factors)
         size(out, mode) == size(factors[mode], 1) || throw(
             DimensionMismatch(
                 "reconstruct_tucker!: output mode $mode has size $(size(out, mode)), " *
@@ -84,7 +84,7 @@ function reconstruct_tucker!(
         acc = zero(T)
         for J in CartesianIndices(core)
             val = core[J]
-            for mode = 1:N
+            for mode in eachindex(factors)
                 val *= factors[mode][I[mode], J[mode]]
             end
             acc += val
diff --git a/src/tucker/sthosvd.jl b/src/tucker/sthosvd.jl
index 2b0e575..1904349 100644
--- a/src/tucker/sthosvd.jl
+++ b/src/tucker/sthosvd.jl
@@ -27,7 +27,7 @@ For a Tucker result with core `S` and factors `U_1, ..., U_d`, this returns
 """
 function reconstruct(td::TuckerResult{T,N}) where {T,N}
     A = td.core
-    for k = 1:N
+    for k in eachindex(td.factors)
         A = mode_n_product(A, td.factors[k], k)
     end
     return A
@@ -48,7 +48,7 @@ dimension is a reasonable default.
 """
 function optimal_mode_order(dims::NTuple{N,Int}, ranks::NTuple{N,Int}) where {N}
     # Sort modes by compression ratio nk/rk (ascending = most compressible first)
-    ratios = [dims[k] / ranks[k] for k = 1:N]
+    ratios = [dims[k] / ranks[k] for k in eachindex(dims)]
     return sortperm(ratios)  # most compressible first
 end
 
@@ -105,7 +105,7 @@ function sthosvd(
     d = N
 
     # Validate inputs
-    for k = 1:d
+    for k in eachindex(ranks)
         @assert 1 <= ranks[k] <= dims[k] "Rank r[$k]=$(ranks[k]) must be in [1, $(dims[k])]"
     end
     @assert sort(processing_order) == 1:d "processing_order must be a permutation of 1:$d"
@@ -213,7 +213,7 @@ function sthosvd(
         rk = max(rk, 1)  # keep at least rank 1
 
         if verbose
-            discarded = rk < length(sigma) ? sqrt(sum(sigma[rk+1:end] .^ 2)) : 0.0
+            discarded = rk < length(sigma) ? sqrt(sum(sigma[(rk+1):end] .^ 2)) : 0.0
             update_progress!(
                 progress,
                 step;
@@ -269,7 +269,7 @@ function thosvd(
     dims = size(A)
     d = N
 
-    for k = 1:d
+    for k in eachindex(ranks)
         @assert 1 <= ranks[k] <= dims[k] "Rank r[$k]=$(ranks[k]) must be in [1, $(dims[k])]"
     end
 
@@ -281,7 +281,7 @@ function thosvd(
     singular_vals = Vector{Vector{T}}(undef, d)
 
     # Step 1: Compute all factor matrices from the ORIGINAL tensor
-    for k = 1:d
+    for k in eachindex(ranks)
         Ak = unfold_mode(A, k)
         F = svd(Ak)
         factors[k] = Matrix(@view F.U[:, 1:ranks[k]])
@@ -296,7 +296,7 @@ function thosvd(
 
     # Step 2: Compute core tensor  S = A ×₁ U₁ᵀ ×₂ U₂ᵀ ⋯ ×_d Udᵀ
     S = copy(A)
-    for k = 1:d
+    for k in eachindex(ranks)
         S = mode_n_product(S, factors[k]', k)
     end
     verbose && finish_progress!(
@@ -318,11 +318,11 @@ end
 """
 function error_bound(td::TuckerResult{T,N}) where {T,N}
     sq_error = zero(T)
-    for k = 1:N
+    for k in eachindex(td.singular_values)
         rk = size(td.core, k)
         sigma = td.singular_values[k]
         if rk < length(sigma)
-            sq_error += sum(sigma[rk+1:end] .^ 2)
+            sq_error += sum(sigma[(rk+1):end] .^ 2)
         end
     end
     return sqrt(sq_error)
diff --git a/test/basic_tests.jl b/test/basic_tests.jl
index f225ccf..7facd9c 100644
--- a/test/basic_tests.jl
+++ b/test/basic_tests.jl
@@ -501,9 +501,12 @@ end
     @test hasproperty(trace_info, :component_trace_iterations)
     @test hasproperty(trace_info, :component_trace_max_delta_history)
     @test hasproperty(trace_info, :component_trace_delta_history)
-    @test hasproperty(trace_info, :component_trace_rgrad_top1_share_history)
-    @test hasproperty(trace_info, :component_trace_rgrad_top3_share_history)
-    @test hasproperty(trace_info, :component_trace_rgrad_effective_components_history)
+    @test hasproperty(trace_info, :component_trace_coordinate_rgrad_top1_share_history)
+    @test hasproperty(trace_info, :component_trace_coordinate_rgrad_top3_share_history)
+    @test hasproperty(
+        trace_info,
+        :component_trace_coordinate_rgrad_effective_components_history,
+    )
     @test hasproperty(trace_info, :component_trace_coordinate_rgrad_energy_history)
     @test hasproperty(trace_info, :component_trace_metric_rgrad_energy_history)
     @test hasproperty(trace_info, :component_trace_ambient_component_velocity_history)
@@ -524,12 +527,13 @@ end
           length(trace_info.component_trace_max_delta_history)
     @test length(trace_info.component_trace_delta_history) ==
           length(trace_info.component_trace_max_delta_history)
-    @test length(trace_info.component_trace_rgrad_top1_share_history) ==
-          length(trace_info.component_trace_iterations)
-    @test length(trace_info.component_trace_rgrad_top3_share_history) ==
+    @test length(trace_info.component_trace_coordinate_rgrad_top1_share_history) ==
           length(trace_info.component_trace_iterations)
-    @test length(trace_info.component_trace_rgrad_effective_components_history) ==
+    @test length(trace_info.component_trace_coordinate_rgrad_top3_share_history) ==
           length(trace_info.component_trace_iterations)
+    @test length(
+        trace_info.component_trace_coordinate_rgrad_effective_components_history,
+    ) == length(trace_info.component_trace_iterations)
     @test length(trace_info.component_trace_metric_rgrad_top1_share_history) ==
           length(trace_info.component_trace_iterations)
     @test length(trace_info.component_trace_ambient_velocity_top1_share_history) ==
@@ -546,26 +550,26 @@ end
     @test all(isfinite, trace_info.component_trace_max_delta_history)
     @test all(
         x -> isnan(x) || -1e-12 <= x <= 1 + 1e-12,
-        trace_info.component_trace_rgrad_top1_share_history,
+        trace_info.component_trace_coordinate_rgrad_top1_share_history,
     )
     @test all(
         x -> isnan(x) || -1e-12 <= x <= 1 + 1e-12,
-        trace_info.component_trace_rgrad_top3_share_history,
+        trace_info.component_trace_coordinate_rgrad_top3_share_history,
     )
     @test all(
         zip(
-            trace_info.component_trace_rgrad_top1_share_history,
-            trace_info.component_trace_rgrad_top3_share_history,
+            trace_info.component_trace_coordinate_rgrad_top1_share_history,
+            trace_info.component_trace_coordinate_rgrad_top3_share_history,
         ),
     ) do (top1, top3)
         isnan(top1) || isnan(top3) || top1 <= top3 + 1e-12
     end
     @test all(
         x -> isnan(x) || 1 <= x <= r,
-        trace_info.component_trace_rgrad_effective_components_history,
+        trace_info.component_trace_coordinate_rgrad_effective_components_history,
     )
-    @test trace_info.component_trace_rgrad_argmax_component_final ==
-          trace_info.component_trace_rgrad_argmax_component_history[end]
+    @test trace_info.component_trace_coordinate_rgrad_argmax_component_final ==
+          trace_info.component_trace_coordinate_rgrad_argmax_component_history[end]
     @test trace_info.component_trace_metric_rgrad_argmax_component_final ==
           trace_info.component_trace_metric_rgrad_argmax_component_history[end]
     @test trace_info.component_trace_ambient_velocity_argmax_component_final ==