Updated default_algtype to properly select the most appropriate algorithm

src/NCMProblem.jl

@@ -36,30 +36,6 @@ struct NCMProblem{T, M, P, K}
         require_matrix(A)
         require_square(A)
         require_real(A)
-
-        mask = normalize_mask(mask)
-
         return new{typeof(A), typeof(mask), typeof(p), typeof(kwargs)}(A, mask, p, kwargs)
     end
 end
-
-
-normalize_mask(::Nothing) = nothing
-
-function normalize_mask(B::BitMatrix)
-    n = require_square(B)
-
-    for i in 1:(n - 1), j in (i + 1):n
-        b = B[i, j] || B[j, i]
-        B[i, j] = b
-        B[j, i] = b
-    end
-
-    for ii in diagind(B)
-        B[ii] = false
-    end
-
-    return B
-end
-
-normalize_mask(B::AbstractMatrix{T}) where {T} = normalize_mask(BitMatrix(B))

src/NCMSolution.jl

@@ -89,7 +89,7 @@ function CommonSolve.solve!(solver::NCMSolver, args...; kwargs...)
         # definiteness perturbs every entry, so re-apply the mask afterwards. The fixed elements
         # (and unit diagonal) take precedence - the result is PD up to O(√eps).
         if sol.solver.mask !== nothing
-            project_fixed!(sol.X, sol.solver.A_orig, sol.solver.mask)
+            project_fixed!(sol.X, solver.A_orig, solver.mask)
             project_unit!(sol.X)
         else
             cov2cor!(sol.X)

src/NCMSolver.jl

@@ -35,13 +35,6 @@ mutable struct NCMSolver{TA, P, Talg, Tc, Ttol, Tm}
     A_orig::TA      # a copy of A, or an alias of A if no mask is given
 end
 
-"""
-    default_algtype(prob)
-
-Get the default algorithm type for a given input matrix.
-"""
-default_algtype(prob::NCMProblem) = prob.mask === nothing ? Newton : AcceleratedAP
-
 """
     init(prob, alg, args...; kwargs...)
 
@@ -94,17 +87,20 @@ function CommonSolve.init(
         kwargs...
     )
     # Resolve the effective mask first: an explicit `mask=` kwarg overrides any mask set on the
-    # problem; otherwise fall back to the problem's (already-normalized) mask, then to nothing.
+    # problem; otherwise fall back to the problem's mask, then finally default to nothing.
     mask = if mask !== nothing
         verbose && println("Using fixed-element mask")
-        normalize_mask(mask)
+        mask
     elseif prob.mask !== nothing
         verbose && println("Using fixed-element mask from the problem")
         prob.mask
     else
         nothing
     end
 
+    # Ensure that the mask is normalized to a BitMatrix or Nothing
+    mask = normalize_mask(mask)
+
     # A non-empty effective mask must be enforced by an algorithm that supports it. Check the
     # *effective* mask (not just the kwarg) so a mask set on the problem is not silently ignored
     # by an unsupported algorithm.
@@ -235,10 +231,9 @@ end
 
 Initialize the solver, and autotune the algorithm to the problem.
 """
-function CommonSolve.init(
-        prob::NCMProblem, algtype::Type{<:NCMAlgorithm}, args...; kwargs...
-    )
-    return init(prob, autotune(algtype, prob), args...; kwargs...)
+function CommonSolve.init(prob::NCMProblem, algtype::Type{<:NCMAlgorithm}, args...; kwargs...)
+    alg = autotune(algtype, prob)
+    return init(prob, alg, args...; kwargs...)
 end
 
 """
@@ -256,5 +251,50 @@ end
 Initialize the solver with the default algorithm autotuned to the problem.
 """
 function CommonSolve.init(prob::NCMProblem, ::Nothing, args...; kwargs...)
-    return init(prob, default_algtype(prob), args...; kwargs...)
+    algtype = default_algtype(prob; kwargs...)
+    return init(prob, algtype, args...; kwargs...)
+end
+
+"""
+    default_algtype(prob; kwargs...)
+
+Get the default algorithm type for a given input matrix.
+"""
+function default_algtype(prob::NCMProblem; mask = nothing, kwargs...)
+    mask = if mask !== nothing
+        mask
+    elseif prob.mask !== nothing
+        prob.mask
+    else
+        nothing
+    end
+
+    return mask === nothing ? Newton : AcceleratedAP
+end
+
+"""
+    normalize_mask(mask)
+
+Normalizes a mask by converting it to a BitMatrix, or leaves as `nothing`. If a matrix is
+given, the resulting BitMatrix is forced to be symmetric. If A[i,j] or A[j,i] is true, then
+the resulting BitMatrix will have a true value in both positions. The diagonal elements are
+always set to false, since the algorithm should handle setting the diagonal elements to 1.
+"""
+function normalize_mask(B::BitMatrix)
+    n = require_square(B)
+
+    for i in 1:(n - 1), j in (i + 1):n
+        b = B[i, j] || B[j, i]
+        B[i, j] = b
+        B[j, i] = b
+    end
+
+    for ii in diagind(B)
+        B[ii] = false
+    end
+
+    return B
 end
+
+normalize_mask(B::AbstractMatrix{T}) where {T} = normalize_mask(BitMatrix(B))
+normalize_mask(::Nothing) = nothing

src/simple_interface.jl

@@ -12,8 +12,7 @@ positive definite, and corrected if it is not.
 When a fixed-element mask is passed (via the `mask` keyword or on the problem), the result has
 an exact unit diagonal and retains the masked elements exactly. Because repairing positive
 definiteness perturbs every entry, strict PD and exact fixed-element feasibility cannot both be
-guaranteed: the fixed elements take precedence, and the result is positive definite up to
-``O(\\sqrt{\\mathrm{eps}})``.
+guaranteed: the fixed elements take precedence, and the result is positive definite up to `√ϵ`.
 
 # Examples
 
@@ -42,8 +41,11 @@ true
 ```
 """
 function nearest_cor!(A, alg; kwargs...)
-    sol = solve(
-        NCMProblem(A),
+    # (#41) Pass any potential mask to the problem so that the proper algorithm can be selected
+    prob = NCMProblem(A; kwargs...)
+
+    solver = init(
+        prob,
         alg;
         alias_A = true,
         fix_sym = true,
@@ -52,6 +54,8 @@ function nearest_cor!(A, alg; kwargs...)
         kwargs...
     )
 
+    sol = solve!(solver)
+
     copyto!(A, sol.X)
     return A
 end
@@ -71,8 +75,7 @@ positive definite, and corrected if it is not.
 When a fixed-element mask is passed (via the `mask` keyword or on the problem), the result has
 an exact unit diagonal and retains the masked elements exactly. Because repairing positive
 definiteness perturbs every entry, strict PD and exact fixed-element feasibility cannot both be
-guaranteed: the fixed elements take precedence, and the result is positive definite up to
-``O(\\sqrt{\\mathrm{eps}})``.
+guaranteed: the fixed elements take precedence, and the result is positive definite up to `√ϵ`.
 
 # Examples
 
@@ -101,8 +104,15 @@ true
 ```
 """
 function nearest_cor(A, alg; kwargs...)
-    sol = solve(
-        NCMProblem(A),
+    # (#41) Pass any potential mask to the problem so that the proper algorithm can be selected
+    prob = NCMProblem(A; kwargs...)
+
+    # `nearest_cor` cannot simply call `nearest_cor!` with `alias_A=false`. For some reason
+    # the `alias_A` keyword is passed properly, but `A` still ends up getting aliased anyway.
+    # The solution is to call `init` with `alias_A=false` independently.
+
+    solver = init(
+        prob,
         alg;
         alias_A = false,
         fix_sym = true,
@@ -111,7 +121,8 @@ function nearest_cor(A, alg; kwargs...)
         kwargs...
     )
 
+    sol = solve!(solver)
+
     return sol.X
 end
-
 nearest_cor(A; kwargs...) = nearest_cor(A, nothing; kwargs...)

test/test_api.jl

@@ -1,4 +1,5 @@
 using Test
+using LinearAlgebra
 using InteractiveUtils
 using NearestCorrelationMatrix
 using NearestCorrelationMatrix.Internals: default_negdef
@@ -46,3 +47,19 @@ cache = init(prob)
 @test_isdefined solve!
 @test_isimplemented solve!(cache)
 @test solve!(cache) isa NCMSolution
+
+# alias_A must be respected
+A = default_negdef()
+prob = NCMProblem(A)
+@test Base.mightalias(A, prob.A)
+solver = init(prob; alias_A = true)
+@test Base.mightalias(prob.A, solver.A)
+solver = init(prob; alias_A = false)
+@test !Base.mightalias(prob.A, solver.A)
+S = Symmetric(A)
+prob = NCMProblem(S)
+@test Base.mightalias(S, prob.A)
+solver = init(prob; alias_A = true)
+@test Base.mightalias(prob.A, solver.A)
+solver = init(prob; alias_A = false)
+@test !Base.mightalias(prob.A, solver.A)

test/test_simple_api.jl

@@ -40,5 +40,15 @@ A = rand(Float16, 4, 4)
 
 # (#41) uses an algorithm that supports masking when a mask is given
 A, m = default_negdef(; include_mask = true)
-@test_broken nearest_cor(A; mask = m)
-@test_broken nearest_cor!(A; mask = m)
+@test_nothrow nearest_cor(A; mask = m)
+@test_nothrow nearest_cor!(A; mask = m)
+
+# nearest_cor must not modify original matrix UNLESS user passes `alias_A=true`
+A = default_negdef()
+Y = nearest_cor(A, AlternatingProjections)
+@test Base.mightalias(Y, A) == false
+@test !isapprox(Y, A)
+A = default_negdef()
+Y = nearest_cor(A, AlternatingProjections; alias_A = true)
+@test Base.mightalias(Y, A) == true
+@test isapprox(Y, A)