Allow D_init and nothing tolerances

.github/workflows/CI.yml

@@ -22,16 +22,7 @@ jobs:
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}
-      - uses: actions/cache@v3
-        env:
-          cache-name: cache-artifacts
-        with:
-          path: ~/.julia/artifacts
-          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
-          restore-keys: |
-            ${{ runner.os }}-test-${{ env.cache-name }}-
-            ${{ runner.os }}-test-
-            ${{ runner.os }}-
+      - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1

src/KSVD.jl

@@ -88,11 +88,12 @@ function ksvd(Y::AbstractMatrix{T}, n_atoms::Int, max_nnz=n_atoms÷10;
               ksvd_update_method = BatchedParallelKSVD{false, T}(; shuffle_indices=true, batch_size_per_thread=1),
               sparse_coding_method = ParallelMatchingPursuit(; max_nnz, rtol=5e-2),
               minibatch_size=nothing,
+              D_init::Union{Nothing, <:AbstractMatrix{T}} = nothing,
               # termination conditions
               maxiters::Int=100, #: The maximum number of iterations to perform. Defaults to 100.
               maxtime::Union{Nothing, <:Real}=nothing,# : The maximum time for solving the nonlinear system of equations. Defaults to nothing which means no time limit. Note that setting a time limit does have a small overhead.
-              abstol::Number=real(oneunit(T)) * (eps(real(one(T))))^(4 // 5), #: The absolute tolerance. Defaults to real(oneunit(T)) * (eps(real(one(T))))^(1 // 2).
-              reltol::Number=real(oneunit(T)) * (eps(real(one(T))))^(4 // 5), #: The relative tolerance. Defaults to real(oneunit(T)) * (eps(real(one(T))))^(1 // 2).
+              abstol::Union{Nothing, <:Real}=real(oneunit(T)) * (eps(real(one(T))))^(4 // 5), #: The absolute tolerance. Defaults to real(oneunit(T)) * (eps(real(one(T))))^(1 // 2).
+              reltol::Union{Nothing, <:Real}=real(oneunit(T)) * (eps(real(one(T))))^(4 // 5), #: The relative tolerance. Defaults to real(oneunit(T)) * (eps(real(one(T))))^(1 // 2).
               nnz_per_col_target::Number=0.0,
               # tracing options
               show_trace::Bool=false,
@@ -104,7 +105,7 @@ function ksvd(Y::AbstractMatrix{T}, n_atoms::Int, max_nnz=n_atoms÷10;
 
     # D is a dictionary matrix that contains atoms for columns.
     @timeit_debug timer "Init dict" begin
-        D = init_dictionary(T, emb_dim, n_atoms)  # size(D) == (n, K)
+        D = (isnothing(D_init) ? init_dictionary(T, emb_dim, n_atoms) : D_init)  # size(D) == (n, K)
         @assert all(≈(1.0), norm.(eachcol(D)))
     end
     X = sparse_coding(sparse_coding_method, Y, D; timer)
@@ -152,7 +153,7 @@ function ksvd(Y::AbstractMatrix{T}, n_atoms::Int, max_nnz=n_atoms÷10;
             termination_condition = :maxiter; break
         elseif !isnothing(maxtime) && (time() - tic) > maxtime
             termination_condition = :maxtime; break
-        elseif length(norm_results) > 1 && isapprox(norm_results[end], norm_results[end-1]; atol=abstol, rtol=reltol)
+        elseif (!isnothing(abstol) && !isnothing(reltol)) && length(norm_results) > 1 && isapprox(norm_results[end], norm_results[end-1]; atol=abstol, rtol=reltol)
             termination_condition = :converged; break
         elseif !isempty(nnz_per_col_results) && last(nnz_per_col_results) <= nnz_per_col_target
             termination_condition = :nnz_per_col_target; break

src/matching_pursuit.jl

@@ -245,8 +245,8 @@ products[t+1] = dictionary' * (residual[t] - dictionary[idx] * a)
     products_abs = abs.(products)  # prealloc
 
     for i in 1:max_iter
-        !isfinite(norm(residual)) && @show norm(residual), residual
-        !isfinite(norm_data) && @show norm_data
+        # @assert(norm(residual)) && @show norm(residual), residual
+        # @assert(isfinite(norm_data))
         if norm(residual)/norm_data < rtol
             return sparsevec(xdict, n_atoms)
         end
@@ -261,7 +261,7 @@ products[t+1] = dictionary' * (residual[t] - dictionary[idx] * a)
 
         a = products[maxindex]
         atom = @view dictionary[:, maxindex]
-        @assert norm(atom) ≈ 1. norm(atom)
+        # @assert norm(atom) ≈ 1. norm(atom)
 
         residual .-= a .* atom
         products .-= a .* @view DtD[:, maxindex]