diff --git a/src/NMFkCluster.jl b/src/NMFkCluster.jl index 96741e5..1ee4c84 100644 --- a/src/NMFkCluster.jl +++ b/src/NMFkCluster.jl @@ -89,7 +89,7 @@ function robustkmeans(X::AbstractMatrix, k::Integer, repeats::Integer=1000; maxi Xn = zerostoepsilon(X) for i = 1:repeats local c_new - @Suppressor.suppress begin + Suppressor.@suppress begin c_new = Clustering.kmeans(X, k; maxiter=maxiter, tol=tol, display=display, distance=distance) end Xd = Distances.pairwise(distance, Xn; dims=2) diff --git a/src/NMFkExecute.jl b/src/NMFkExecute.jl index 5869888..5cd1369 100644 --- a/src/NMFkExecute.jl +++ b/src/NMFkExecute.jl @@ -14,7 +14,7 @@ function execute(X::AbstractArray{T,N}, nkrange::AbstractRange{Int}, nNMF::Integ end @info("Results") for nk in nkrange - println("Signals: $(@Printf.sprintf("%2d", nk)) Fit: $(@Printf.sprintf("%12.7g", fitquality[nk])) Silhouette: $(@Printf.sprintf("%12.7g", robustness[nk])) AIC: $(@Printf.sprintf("%12.7g", aic[nk]))") + println("Signals: $(Printf.@sprintf("%2d", nk)) Fit: $(Printf.@sprintf("%12.7g", fitquality[nk])) Silhouette: $(Printf.@sprintf("%12.7g", robustness[nk])) AIC: $(Printf.@sprintf("%12.7g", aic[nk]))") end kopt = getk(nkrange, robustness[nkrange], cutoff) if isnothing(kopt) @@ -72,7 +72,7 @@ function execute(X::Union{AbstractMatrix{T},AbstractArray{T}}, nk::Integer, nNMF W, H, fitquality, robustness, aic = NMFk.execute_run(X, nk, nNMF; clusterWmatrix=clusterWmatrix, veryquiet=veryquiet, resultdir=resultdir, casefilename=casefilename, kw...) end so = signalorder(W, H) - !veryquiet && println("Signals: $(@Printf.sprintf("%2d", nk)) Fit: $(@Printf.sprintf("%12.7g", fitquality)) Silhouette: $(@Printf.sprintf("%12.7g", robustness)) AIC: $(@Printf.sprintf("%12.7g", aic)) Signal order: $(so)") + !veryquiet && println("Signals: $(Printf.@sprintf("%2d", nk)) Fit: $(Printf.@sprintf("%12.7g", fitquality)) Silhouette: $(Printf.@sprintf("%12.7g", robustness)) AIC: $(Printf.@sprintf("%12.7g", aic)) Signal order: $(so)") if save && casefilename != "" filename = joinpathcheck(resultdir, "$casefilename-$nk-$nNMF.jld") recursivemkdir(filename) diff --git a/src/NMFkHelpers.jl b/src/NMFkHelpers.jl index e3c6d06..9f7f51a 100644 --- a/src/NMFkHelpers.jl +++ b/src/NMFkHelpers.jl @@ -387,10 +387,10 @@ function estimateflip_permutedims(X::AbstractMatrix{T}, Y::AbstractMatrix{T}, A: vflip = falses(nparam) for i = 1:nparam local H1 - @Suppressor.suppress W, H1, of, sil, aic = NMFk.execute(permutedims(Y[:,i]), nk, nNNF; Winit=permutedims(X[:,i]), Wfixed=true, save=save, method=method, regularizationweight=regularizationweight, kw...); + Suppressor.@suppress W, H1, of, sil, aic = NMFk.execute(permutedims(Y[:,i]), nk, nNNF; Winit=permutedims(X[:,i]), Wfixed=true, save=save, method=method, regularizationweight=regularizationweight, kw...); a = NMFk.normnan(permutedims(A) .- (permutedims(B) * H1)) local H2 - @Suppressor.suppress W, H2, of, sil, aic = NMFk.execute(permutedims(NMFk.flip(Y[:,i])), nk, nNNF; Winit=permutedims(NMFk.flip(X[:,i])), Wfixed=true, save=save, method=method, regularizationweight=regularizationweight, kw...); + Suppressor.@suppress W, H2, of, sil, aic = NMFk.execute(permutedims(NMFk.flip(Y[:,i])), nk, nNNF; Winit=permutedims(NMFk.flip(X[:,i])), Wfixed=true, save=save, method=method, regularizationweight=regularizationweight, kw...); b = NMFk.normnan(permutedims(A) .- (permutedims(B) * H2)) vflip[i] = a < b ? false : true end @@ -407,10 +407,10 @@ function estimateflip(X::AbstractMatrix{T}, Y::AbstractMatrix{T}, A::AbstractMat vflip = falses(nparam) for i = 1:nparam local H1 - @Suppressor.suppress W, H1, of, sil, aic = NMFk.execute(Y[i:i,:], nk, nNNF; Winit=X[i:i,:], Wfixed=true, save=save, method=method, regularizationweight=regularizationweight, kw...); + Suppressor.@suppress W, H1, of, sil, aic = NMFk.execute(Y[i:i,:], nk, nNNF; Winit=X[i:i,:], Wfixed=true, save=save, method=method, regularizationweight=regularizationweight, kw...); a = NMFk.normnan(B .- (A * H1)) local H2 - @Suppressor.suppress W, H2, of, sil, aic = NMFk.execute(NMFk.flip(Y[i:i,:]), nk, nNNF; Winit=NMFk.flip(X[i:i,:]), Wfixed=true, save=save, method=method, regularizationweight=regularizationweight, kw...); + Suppressor.@suppress W, H2, of, sil, aic = NMFk.execute(NMFk.flip(Y[i:i,:]), nk, nNNF; Winit=NMFk.flip(X[i:i,:]), Wfixed=true, save=save, method=method, regularizationweight=regularizationweight, kw...); b = NMFk.normnan(B .- (A * H2)) vflip[i] = a < b ? false : true end diff --git a/src/NMFkIO.jl b/src/NMFkIO.jl index 8a930d2..ac59b4d 100644 --- a/src/NMFkIO.jl +++ b/src/NMFkIO.jl @@ -42,7 +42,7 @@ function load(nk::Integer, nNMF::Integer=10; type::DataType=Float64, dim::Intege if isfile(filename) W, H, fitquality, robustness, aic = JLD.load(filename, "W", "H", "fit", "robustness", "aic") so = signalorder(W, H) - !quiet && println("Signals: $(@Printf.sprintf("%2d", nk)) Fit: $(@Printf.sprintf("%12.7g", fitquality)) Silhouette: $(@Printf.sprintf("%12.7g", robustness)) AIC: $(@Printf.sprintf("%12.7g", aic)) Signal order: $(so)") + !quiet && println("Signals: $(Printf.@sprintf("%2d", nk)) Fit: $(Printf.@sprintf("%12.7g", fitquality)) Silhouette: $(Printf.@sprintf("%12.7g", robustness)) AIC: $(Printf.@sprintf("%12.7g", aic)) Signal order: $(so)") return W[:,so], H[so,:], fitquality, robustness, aic else !quiet && @warn("File named $filename is missing!") diff --git a/src/NMFkJuMP.jl b/src/NMFkJuMP.jl index 2b2f304..4cd6f25 100644 --- a/src/NMFkJuMP.jl +++ b/src/NMFkJuMP.jl @@ -122,32 +122,32 @@ function jump(X::AbstractArray{T}, nk::Int; method::Symbol=:nlopt, algorithm::Sy W = Winit else constrainW && normalize!(Winit) - @JuMP.variable(m, W[i=1:nummixtures, j=1:nk], start=convert(T, Winit[i, j])) - Wnonneg && @JuMP.constraint(m, W .>= 0) - constrainW && @JuMP.constraint(m, W .<= 1) + JuMP.@variable(m, W[i=1:nummixtures, j=1:nk], start=convert(T, Winit[i, j])) + Wnonneg && JuMP.@constraint(m, W .>= 0) + constrainW && JuMP.@constraint(m, W .<= 1) end if Hfixed H = Hinit else constrainH && normalize!(Hinit) - @JuMP.variable(m, H[i=1:nk, j=1:numconstituents], start=convert(T, Hinit[i, j])) - Hnonneg && @JuMP.constraint(m, H .>= 0) - constrainH && @JuMP.constraint(m, H .<= 1) + JuMP.@variable(m, H[i=1:nk, j=1:numconstituents], start=convert(T, Hinit[i, j])) + Hnonneg && JuMP.@constraint(m, H .>= 0) + constrainH && JuMP.@constraint(m, H .<= 1) end if kullbackleibler smallnumber = eps(Float64) - @JuMP.NLobjective(m, Min, sum(X[i, j] * (log(smallnumber + X[i, j]) - log(smallnumber + sum(W[i, k] * H[k, j] for k = 1:nk))) - X[i, j] + sum(W[i, k] * H[k, j] for k = 1:nk) for i=1:nummixtures, j=1:numconstituents)) + JuMP.@NLobjective(m, Min, sum(X[i, j] * (log(smallnumber + X[i, j]) - log(smallnumber + sum(W[i, k] * H[k, j] for k = 1:nk))) - X[i, j] + sum(W[i, k] * H[k, j] for k = 1:nk) for i=1:nummixtures, j=1:numconstituents)) else if regularizationweight == 0. - @JuMP.NLobjective(m, Min, + JuMP.@NLobjective(m, Min, sum(sum(obsweights[i, j] * (sum(W[i, k] * H[k, j] for k=1:nk) - X[i, j])^2 for i=1:nummixtures) for j=1:numconstituents)) else if Hfixed - @JuMP.NLobjective(m, Min, + JuMP.@NLobjective(m, Min, regularizationweight * sum(sum(log(1. + W[i, j])^2 for i=1:numconstituents) for j=1:nk) / nk + sum(sum(obsweights[i, j] * (sum(W[i, k] * H[k, j] for k=1:nk) - X[i, j])^2 for i=1:nummixtures) for j=1:numconstituents)) else - @JuMP.NLobjective(m, Min, + JuMP.@NLobjective(m, Min, regularizationweight * sum(sum(log(1. + H[i, j])^2 for i=1:nk) for j=1:numconstituents) / nk + sum(sum(obsweights[i, j] * (sum(W[i, k] * H[k, j] for k=1:nk) - X[i, j])^2 for i=1:nummixtures) for j=1:numconstituents)) end @@ -160,7 +160,7 @@ function jump(X::AbstractArray{T}, nk::Int; method::Symbol=:nlopt, algorithm::Sy jumpvariables = JuMP.all_variables(m) jumpvalues = JuMP.start_value.(jumpvariables) if quiet - @Suppressor.suppress JuMP.optimize!(m) + Suppressor.@suppress JuMP.optimize!(m) else JuMP.optimize!(m) end @@ -194,7 +194,7 @@ function jump(X::AbstractArray{T}, nk::Int; method::Symbol=:nlopt, algorithm::Sy frame += 1 end if quiet - @Suppressor.suppress JuMP.optimize!(m) + Suppressor.@suppress JuMP.optimize!(m) else JuMP.optimize!(m) end diff --git a/src/NMFkMap.jl b/src/NMFkMap.jl index 96ff12a..b3d277d 100644 --- a/src/NMFkMap.jl +++ b/src/NMFkMap.jl @@ -82,7 +82,7 @@ function plotmap(X::AbstractMatrix, fips::AbstractVector, dim::Integer=1, signal ltitle = "$(titletext) $(signalidtext)" end end - p = @VegaLite.vlplot( + p = VegaLite.@vlplot( title=ttitle, :geoshape, width=500, height=300, @@ -119,7 +119,7 @@ function plotmap(X::AbstractVector, fips::AbstractVector; us10m=VegaDatasets.dat @assert length(X) == length(fips) nc = length(unique(sort(X))) + 1 df = DataFrames.DataFrame(FIPS=[fips[goodcounties]; fips[.!goodcounties]], Z=[X; zeros(sum(.!goodcounties))]) - p = @VegaLite.vlplot( + p = VegaLite.@vlplot( :geoshape, width=500, height=300, data={ diff --git a/src/NMFkMapping.jl b/src/NMFkMapping.jl index f046578..0bfa229 100644 --- a/src/NMFkMapping.jl +++ b/src/NMFkMapping.jl @@ -9,13 +9,13 @@ function mapping_old(X::AbstractMatrix{T}, Y::AbstractMatrix{T}, A::AbstractMatr np = size(X, 2) local W1, H1, of1, sil1, aic1 @info "Mapping matrix size: $nk x $(size(Y, 1))" - @Suppressor.suppress W1, H1, of1, sil1, aic1 = NMFk.execute(permutedims(Y), nk, nNNF; Winit=permutedims(X), Wfixed=true, save=save, method=method, regularizationweight=regularizationweight, kw...) + Suppressor.@suppress W1, H1, of1, sil1, aic1 = NMFk.execute(permutedims(Y), nk, nNNF; Winit=permutedims(X), Wfixed=true, save=save, method=method, regularizationweight=regularizationweight, kw...) a = NMFk.normnan(permutedims(B) .- (permutedims(A) * H1)) vflip = NMFk.estimateflip_permutedims(X, Y, A, B) Xn = hcat(map(i->vflip[i] ? NMFk.flip(X[:,i]) : X[:,i], 1:np)...) Yn = hcat(map(i->vflip[i] ? NMFk.flip(Y[:,i]) : Y[:,i], 1:np)...) local W2, H2, of2, sil2, aic2 - @Suppressor.suppress W2, H2, of2, sil2, aic2 = NMFk.execute(permutedims(Yn), nk, nNNF; Winit=permutedims(Xn), Wfixed=true, save=save, method=method, regularizationweight=regularizationweight, kw...) + Suppressor.@suppress W2, H2, of2, sil2, aic2 = NMFk.execute(permutedims(Yn), nk, nNNF; Winit=permutedims(Xn), Wfixed=true, save=save, method=method, regularizationweight=regularizationweight, kw...) b = NMFk.normnan(permutedims(B) .- (permutedims(A) * H2)) if a < b return W1, H1, of1, sil1, aic1 @@ -51,7 +51,7 @@ function mapping(X::AbstractMatrix{T}, Y::AbstractMatrix{T}, A::AbstractMatrix{T @info "Mapping matrix size: $nk x $(size(Y, 2))" X[inan] .= 0 local W1, H1, of1, sil1, aic1 - @Suppressor.suppress W1, H1, of1, sil1, aic1 = NMFk.execute(Y, nk, nNNF; Winit=X, Wfixed=true, save=save, method=method, kw..., kwx...) + Suppressor.@suppress W1, H1, of1, sil1, aic1 = NMFk.execute(Y, nk, nNNF; Winit=X, Wfixed=true, save=save, method=method, kw..., kwx...) iz = vec(NMFk.maximumnan(Y; dims=1) .== 0) H1[:, iz] .= 0 if fliptest @@ -60,7 +60,7 @@ function mapping(X::AbstractMatrix{T}, Y::AbstractMatrix{T}, A::AbstractMatrix{T Xn = permutedims(hcat(map(i->vflip[i] ? NMFk.flip(X[i,:]) : X[i,:], 1:np)...)) Yn = permutedims(hcat(map(i->vflip[i] ? NMFk.flip(Y[i,:]) : Y[i,:], 1:np)...)) local W2, H2, of2, sil2, aic2 - @Suppressor.suppress W2, H2, of2, sil2, aic2 = NMFk.execute(Yn, nk, nNNF; Winit=Xn, Wfixed=true, save=save, method=method, kw..., kwx...) + Suppressor.@suppress W2, H2, of2, sil2, aic2 = NMFk.execute(Yn, nk, nNNF; Winit=Xn, Wfixed=true, save=save, method=method, kw..., kwx...) iz = vec(NMFk.maximumnan(Yn; dims=1) .== 0) H2[:, iz] .= 0 b = NMFk.normnan(B .- (A * H2)) diff --git a/src/NMFkMixMatrix.jl b/src/NMFkMixMatrix.jl index 2cb21f3..4a9d823 100644 --- a/src/NMFkMixMatrix.jl +++ b/src/NMFkMixMatrix.jl @@ -93,21 +93,21 @@ function mixmatchdata(concentrations_in::AbstractMatrix{T}, numbuckets::Int; met m = JuMP.Model(NLopt.Optimizer) JuMP.set_optimizer_attributes(m, "algorithm" => algorithm, "maxeval" => maxiter) # "xtol_abs" => tolX, "ftol_abs" => tol end - @JuMP.variable(m, mixer[i=1:nummixtures, j=1:numbuckets], start = convert(T, Winit[i, j])) - @JuMP.variable(m, buckets[i=1:numbuckets, j=1:numconstituents], start = convert(T, Hinit[i, j])) + JuMP.@variable(m, mixer[i=1:nummixtures, j=1:numbuckets], start = convert(T, Winit[i, j])) + JuMP.@variable(m, buckets[i=1:numbuckets, j=1:numconstituents], start = convert(T, Hinit[i, j])) if !normalize - @JuMP.constraint(m, buckets .>= 0) + JuMP.@constraint(m, buckets .>= 0) end - @JuMP.constraint(m, mixer .>= 0) + JuMP.@constraint(m, mixer .>= 0) for i = 1:nummixtures - @JuMP.constraint(m, sum(mixer[i, j] for j=1:numbuckets) == 1.) + JuMP.@constraint(m, sum(mixer[i, j] for j=1:numbuckets) == 1.) end if sizeof(ratios) == 0 - @JuMP.NLobjective(m, Min, + JuMP.@NLobjective(m, Min, regularizationweight * sum(sum(log(1. + buckets[i, j])^2 for i=1:numbuckets) for j=1:numconstituents) / numbuckets + sum(sum(concweights[i, j] * (sum(mixer[i, k] * buckets[k, j] for k=1:numbuckets) - concentrations[i, j])^2 for i=1:nummixtures) for j=1:numconstituents)) else - @JuMP.NLobjective(m, Min, + JuMP.@NLobjective(m, Min, regularizationweight * sum(sum(log(1. + buckets[i, j])^2 for i=1:numbuckets) for j=1:numconstituents) / numbuckets + sum(sum(concweights[i, j] * (sum(mixer[i, k] * buckets[k, j] for k=1:numbuckets) - concentrations[i, j])^2 for i=1:nummixtures) for j=1:numconstituents) + sum(sum(ratiosweightmatrix[i, j] * @@ -122,7 +122,7 @@ function mixmatchdata(concentrations_in::AbstractMatrix{T}, numbuckets::Int; met jumpvariables = JuMP.all_variables(m) jumpvalues = JuMP.start_value.(jumpvariables) if quiet - @Suppressor.suppress JuMP.optimize!(m) + Suppressor.@suppress JuMP.optimize!(m) else JuMP.optimize!(m) end @@ -143,7 +143,7 @@ function mixmatchdata(concentrations_in::AbstractMatrix{T}, numbuckets::Int; met while LinearAlgebra.norm(jumpvalues - JuMP.value.(jumpvariables)) > tolX && ofbest > tol && baditers < maxbaditers && reattempts < maxreattempts jumpvalues = JuMP.value.(jumpvariables) if quiet - @Suppressor.suppress JuMP.optimize!(m) + Suppressor.@suppress JuMP.optimize!(m) else JuMP.optimize!(m) end @@ -308,31 +308,31 @@ function mixmatchdeltas(concentrations_in::AbstractMatrix{T}, deltas_in::Abstrac m = JuMP.Model(NLopt.Optimizer) JuMP.set_optimizer_attributes(m, "algorithm" => algorithm, "maxeval" => maxiter, "xtol_abs" => tolX, "ftol_abs" => tol) end - @JuMP.variable(m, mixer[i=1:nummixtures, j=1:numbuckets], start = convert(T, Winit[i, j])) - @JuMP.variable(m, buckets[i=1:numbuckets, j=1:numconstituents], start = convert(T, Hinit[i, j])) - @JuMP.variable(m, bucketdeltas[i=1:numbuckets, j=1:numdeltas], start = convert(T, Hinitd[i, j])) - @JuMP.constraint(m, buckets .>= 0) - @JuMP.constraint(m, mixer .>= 0) + JuMP.@variable(m, mixer[i=1:nummixtures, j=1:numbuckets], start = convert(T, Winit[i, j])) + JuMP.@variable(m, buckets[i=1:numbuckets, j=1:numconstituents], start = convert(T, Hinit[i, j])) + JuMP.@variable(m, bucketdeltas[i=1:numbuckets, j=1:numdeltas], start = convert(T, Hinitd[i, j])) + JuMP.@constraint(m, buckets .>= 0) + JuMP.@constraint(m, mixer .>= 0) for i = 1:nummixtures - @JuMP.constraint(m, sum(mixer[i, j] for j=1:numbuckets) == 1.) + JuMP.@constraint(m, sum(mixer[i, j] for j=1:numbuckets) == 1.) end #= for i = 1:numbuckets for j = 1:numconstituents if i != 1 || j != 1 - @JuMP.constraint(m, buckets[i, j] == Hinit[i, j]) # Fix buckets for testing + JuMP.@constraint(m, buckets[i, j] == Hinit[i, j]) # Fix buckets for testing end end end for i = 1:numbuckets for j = 1:numdeltas #if i != 1 || j != 1 - @JuMP.constraint(m, bucketdeltas[i, j] == Hinitd[i, j]) # Fix buckets for testing + JuMP.@constraint(m, bucketdeltas[i, j] == Hinitd[i, j]) # Fix buckets for testing #end end end =# - @JuMP.NLobjective(m, Min, + JuMP.@NLobjective(m, Min, regularizationweight * sum(sum(log(1. + buckets[i, j])^2 for i=1:numbuckets) for j=1:numconstituents) / numbuckets + regularizationweight * sum(sum(log(1. + abs(bucketdeltas[i, j]))^2 for i=1:numbuckets) for j=1:numdeltas) / numbuckets + sum(sum(concweights[i, j] * (concentrations[i, j] - (sum(mixer[i, k] * buckets[k, j] for k=1:numbuckets)))^2 for i=1:nummixtures) for j=1:numconstituents) + @@ -382,21 +382,21 @@ function mixmatchwaterdeltas(deltas::AbstractMatrix{T}, numbuckets::Int; method: m = JuMP.Model(Ipopt.Optimizer) JuMP.set_optimizer_attributes(m, "max_iter" => maxiter, "print_level" => verbosity) if random - @JuMP.variable(m, mixer[1:nummixtures, 1:numbuckets] >= 0., start=randn(T)) - @JuMP.variable(m, buckets[1:numbuckets, 1:numconstituents], start=maxdeltaguess * rand(T)) + JuMP.@variable(m, mixer[1:nummixtures, 1:numbuckets] >= 0., start=randn(T)) + JuMP.@variable(m, buckets[1:numbuckets, 1:numconstituents], start=maxdeltaguess * rand(T)) else - @JuMP.variable(m, mixer[1:nummixtures, 1:numbuckets] >= 0.) - @JuMP.variable(m, buckets[1:numbuckets, 1:numconstituents]) + JuMP.@variable(m, mixer[1:nummixtures, 1:numbuckets] >= 0.) + JuMP.@variable(m, buckets[1:numbuckets, 1:numconstituents]) end - @JuMP.constraint(m, mixer .<= 1.) + JuMP.@constraint(m, mixer .<= 1.) for i in axes(deltas, 1) - @JuMP.constraint(m, sum(mixer[i, j] for j=1:numbuckets) == 1.) + JuMP.@constraint(m, sum(mixer[i, j] for j=1:numbuckets) == 1.) end concweights = ones(T, size(deltas)) nans = isnan(deltas) deltas[nans] = 0 concweights[nans] = 0 - @JuMP.NLobjective(m, Min, + JuMP.@NLobjective(m, Min, regularizationweight * sum(sum((buckets[i, j] - bucketmeans[i, j])^2 for i=1:numbuckets) for j=1:numconstituents) / numbuckets + sum(sum(concweights[i, j] * (sum(mixer[i, k] * buckets[k, j] for k=1:numbuckets) - deltas[i, j])^2 for i=1:nummixtures) for j=1:numconstituents)) if method == :ipopt diff --git a/src/NMFkMixTensor.jl b/src/NMFkMixTensor.jl index 80a301e..f53eac4 100644 --- a/src/NMFkMixTensor.jl +++ b/src/NMFkMixTensor.jl @@ -51,24 +51,24 @@ function mixmatchdata(concentrations::AbstractArray{T, 3}, numbuckets::Int; meth m = JuMP.Model(NLopt.Optimizer) JuMP.set_optimizer_attributes(m, "algorithm" => algorithm, "maxeval" => maxiter) # "xtol_abs" => tolX, "ftol_abs" => tol end - @JuMP.variable(m, mixer[i=1:nummixtures, j=1:numbuckets, k=1:ntimes], start=convert(T, Winit[i, j, k])) - @JuMP.variable(m, buckets[i=1:numbuckets, j=1:numconstituents], start=convert(T, Hinit[i, j])) + JuMP.@variable(m, mixer[i=1:nummixtures, j=1:numbuckets, k=1:ntimes], start=convert(T, Winit[i, j, k])) + JuMP.@variable(m, buckets[i=1:numbuckets, j=1:numconstituents], start=convert(T, Hinit[i, j])) if !normalize - @JuMP.constraint(m, buckets .>= 0) + JuMP.@constraint(m, buckets .>= 0) end - @JuMP.constraint(m, mixer .>= 0) + JuMP.@constraint(m, mixer .>= 0) for k = 1:ntimes for i = 1:nummixtures - @JuMP.constraint(m, sum(mixer[i, j, k] for j=1:numbuckets) == 1.) + JuMP.@constraint(m, sum(mixer[i, j, k] for j=1:numbuckets) == 1.) end end - @JuMP.NLobjective(m, Min, + JuMP.@NLobjective(m, Min, regularizationweight * sum(sum(log(1. + buckets[i, j])^2 for i=1:numbuckets) for j=1:numconstituents) / numbuckets + sum(sum(sum(concweights[i, j, t] * (sum(mixer[i, k, t] * buckets[k, j] for k=1:numbuckets) - concentrations[i, j, t])^2 for i=1:nummixtures) for j=1:numconstituents) for t=1:ntimes)) jumpvariables = JuMP.all_variables(m) jumpvalues = JuMP.start_value.(jumpvariables) if quiet - @Suppressor.suppress JuMP.optimize!(m) + Suppressor.@suppress JuMP.optimize!(m) else JuMP.optimize!(m) end @@ -84,7 +84,7 @@ function mixmatchdata(concentrations::AbstractArray{T, 3}, numbuckets::Int; meth while LinearAlgebra.norm(jumpvalues - JuMP.value.(jumpvariables)) > tolX && ofbest > tol && baditers < maxbaditers && reattempts < maxreattempts jumpvalues = JuMP.value.(jumpvariables) if quiet - @Suppressor.suppress JuMP.optimize!(m) + Suppressor.@suppress JuMP.optimize!(m) else JuMP.optimize!(m) end diff --git a/src/NMFkPlot.jl b/src/NMFkPlot.jl index 2a8b255..436d71d 100644 --- a/src/NMFkPlot.jl +++ b/src/NMFkPlot.jl @@ -498,8 +498,8 @@ function plot2dmatrixcomponents(M::AbstractMatrix, dim::Integer=1; quiet::Bool=f return ff end -"Convert `@Printf.sprintf` macro into `sprintf` function" -sprintf(args...) = eval(:@Printf.sprintf($(args...))) +"Convert `Printf.@sprintf` macro into `sprintf` function" +sprintf(args...) = eval(:Printf.@sprintf($(args...))) "Generate Sankey plots" function sankey(c1::AbstractVector, c2::AbstractVector, t1::AbstractString, t2::AbstractString; htmlfile::AbstractString="", pdffile::AbstractString="") diff --git a/src/NMFkRegression.jl b/src/NMFkRegression.jl index 4b96375..aa6e040 100644 --- a/src/NMFkRegression.jl +++ b/src/NMFkRegression.jl @@ -44,7 +44,7 @@ function regression(P::Array{T}, Mtrain::AbstractMatrix{T}, Mpredict::AbstractMa else Wnonneg = false end - @Suppressor.suppress W, H, of, sil, aic = NMFk.execute(permutedims(P[:,:,k]), nk, 1; Wnonneg=Wnonneg, Hinit=[Mtrain ones(size(Mtrain, 1))]', Hfixed=true, method=:ipopt, regularizationweight=0.) + Suppressor.@suppress W, H, of, sil, aic = NMFk.execute(permutedims(P[:,:,k]), nk, 1; Wnonneg=Wnonneg, Hinit=[Mtrain ones(size(Mtrain, 1))]', Hfixed=true, method=:ipopt, regularizationweight=0.) # @show NMFk.normnan((W * [Mtrain ones(size(Mtrain, 1))]')' .- P[:,:,k]) Xe[:,:,k] = permutedims(W * permutedims([Mpredict ones(size(Mpredict, 1))])) end diff --git a/src/NMFkUncertainty.jl b/src/NMFkUncertainty.jl index 971b4d2..9c467a1 100644 --- a/src/NMFkUncertainty.jl +++ b/src/NMFkUncertainty.jl @@ -18,7 +18,7 @@ function uncertainty(X::AbstractArray{T,N}, nk::Integer, nreruns::Integer, nNMF: !quiet && @info("Uncertainty run $(i) out of $(nreruns):") casefilenamemod = save == true ? casefilename * "_$i" : "" if quiet - @Suppressor.suppress W[i], H[i], fitquality[i], robustness[i], aic[i] = NMFk.execute(X[1:window,:], nk, nNMF; kw..., save=save, resultdir=resultdir, casefilename=casefilenamemod) + Suppressor.@suppress W[i], H[i], fitquality[i], robustness[i], aic[i] = NMFk.execute(X[1:window,:], nk, nNMF; kw..., save=save, resultdir=resultdir, casefilename=casefilenamemod) else W[i], H[i], fitquality[i], robustness[i], aic[i] = NMFk.execute(X[1:window,:], nk, nNMF; kw..., save=save, resultdir=resultdir, casefilename=casefilenamemod) end @@ -30,7 +30,7 @@ function uncertainty(X::AbstractArray{T,N}, nk::Integer, nreruns::Integer, nNMF: @info("Uncertainty results stage #1:") end for i in 1:nreruns - println("Run: $(@Printf.sprintf("%5d", i)) Fit: $(@Printf.sprintf("%12.7g", fitquality[i])) Silhouette: $(@Printf.sprintf("%12.7g", robustness[i])) AIC: $(@Printf.sprintf("%12.7g", aic[i]))") + println("Run: $(Printf.@sprintf("%5d", i)) Fit: $(Printf.@sprintf("%12.7g", fitquality[i])) Silhouette: $(Printf.@sprintf("%12.7g", robustness[i])) AIC: $(Printf.@sprintf("%12.7g", aic[i]))") end end if window != size(X,1) @@ -38,7 +38,7 @@ function uncertainty(X::AbstractArray{T,N}, nk::Integer, nreruns::Integer, nNMF: !quiet && @info("Uncertainty run stage #2 $(i) out of $(nreruns):") casefilenamemod = save == true ? casefilename * "_stage2_$i" : "" if quiet - @Suppressor.suppress W[i], H[i], fitquality[i], robustness[i], aic[i] = NMFk.execute(X[1:maxwindow,:], nk; Hinit=convert.(T, H[i]), Hfixed=true, kw..., save=save, resultdir=resultdir, casefilename=casefilenamemod) + Suppressor.@suppress W[i], H[i], fitquality[i], robustness[i], aic[i] = NMFk.execute(X[1:maxwindow,:], nk; Hinit=convert.(T, H[i]), Hfixed=true, kw..., save=save, resultdir=resultdir, casefilename=casefilenamemod) else W[i], H[i], fitquality[i], robustness[i], aic[i] = NMFk.execute(X[1:maxwindow,:], nk; Hinit=convert.(T, H[i]), Hfixed=true, kw..., save=save, resultdir=resultdir, casefilename=casefilenamemod)s end @@ -46,7 +46,7 @@ function uncertainty(X::AbstractArray{T,N}, nk::Integer, nreruns::Integer, nNMF: if !quiet @info("Uncertainty results stage #2:") for i in 1:nreruns - println("Run: $(@Printf.sprintf("%5d", i)) Fit: $(@Printf.sprintf("%12.7g", fitquality[i])) Silhouette: $(@Printf.sprintf("%12.7g", robustness[i])) AIC: $(@Printf.sprintf("%12.7g", aic[i]))") + println("Run: $(Printf.@sprintf("%5d", i)) Fit: $(Printf.@sprintf("%12.7g", fitquality[i])) Silhouette: $(Printf.@sprintf("%12.7g", robustness[i])) AIC: $(Printf.@sprintf("%12.7g", aic[i]))") end end end diff --git a/test/runtests.jl b/test/runtests.jl index 5aeede4..e649a22 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,7 +4,7 @@ import Random import Suppressor import LinearAlgebra -@Test.testset "NMFk" begin +Test.@testset "NMFk" begin function runtest(concs::AbstractMatrix, buckets::AbstractMatrix, ratios::Array{Float32, 2}=Array{Float32, 2}(undef, 0, 0), ratioindices::Union{Array{Int, 1},Array{Int, 2}}=Array{Int, 2}(undef, 0, 0); conccomponents=collect(1:size(concs, 2)), ratiocomponents=Int[], tol::Number=0.0000001) numbuckets = size(buckets, 1) @@ -14,9 +14,9 @@ function runtest(concs::AbstractMatrix, buckets::AbstractMatrix, ratios::Array{F predictedconcs = mixerestimate * bucketestimate predictedconcs[idxnan] .= 0 if length(conccomponents) > 0 - @Test.test LinearAlgebra.norm(predictedconcs[:, conccomponents] - concs[:, conccomponents], 2) / LinearAlgebra.norm(concs[:, conccomponents], 2) < 1 # fit the data within 1% + Test.@test LinearAlgebra.norm(predictedconcs[:, conccomponents] - concs[:, conccomponents], 2) / LinearAlgebra.norm(concs[:, conccomponents], 2) < 1 # fit the data within 1% for j in axes(buckets, 1) - @Test.test minimum(map(i->LinearAlgebra.norm(buckets[i, conccomponents] - bucketestimate[j, conccomponents]) / LinearAlgebra.norm(buckets[i, conccomponents], 2), 1:size(buckets, 1))) < 1 # reproduce the buckets within 30% + Test.@test minimum(map(i->LinearAlgebra.norm(buckets[i, conccomponents] - bucketestimate[j, conccomponents]) / LinearAlgebra.norm(buckets[i, conccomponents], 2), 1:size(buckets, 1))) < 1 # reproduce the buckets within 30% end end checkratios(mixerestimate, bucketestimate, ratios, ratiocomponents) @@ -31,8 +31,8 @@ function checkratios(mixerestimate::AbstractMatrix, bucketestimate::AbstractMatr for i in axes(mixerestimate, 1) for j = 1:numberofratios ratioratio = predictedconcs[i, ratiocomponents[j, 1]] / predictedconcs[i, ratiocomponents[j, 2]] / ratios[i, j] - @Test.test ratioratio > .4 # get the ratio within a factor of 2 - @Test.test ratioratio < 4. + Test.@test ratioratio > .4 # get the ratio within a factor of 2 + Test.@test ratioratio < 4. end end end @@ -141,13 +141,13 @@ W = [a b] H = [.1 1 0 0 .1; 0 0 .1 .5 .2] X = W * H X = [a a*10 b b*5 a+b*2] -@Suppressor.suppress global We, He, p, s = NMFk.execute(X, 2, 10; save=false, method=:ipopt, tolX=1e-3, tol=1e-12) -@Test.test isapprox(p, 0; atol=1e-3) -@Test.test isapprox(s, 1; rtol=1e-1) -@Test.test isapprox(He[1,2] / He[1,1], 10; rtol=1e-3) -@Test.test isapprox(He[2,2] / He[2,1], 10; rtol=1e-3) -@Test.test isapprox(He[2,4] / He[2,3], 5; rtol=1e-3) -@Test.test isapprox(He[1,4] / He[1,3], 5; rtol=1e-3) +Suppressor.@suppress global We, He, p, s = NMFk.execute(X, 2, 10; save=false, method=:ipopt, tolX=1e-3, tol=1e-12) +Test.@test isapprox(p, 0; atol=1e-3) +Test.@test isapprox(s, 1; rtol=1e-1) +Test.@test isapprox(He[1,2] / He[1,1], 10; rtol=1e-3) +Test.@test isapprox(He[2,2] / He[2,1], 10; rtol=1e-3) +Test.@test isapprox(He[2,4] / He[2,3], 5; rtol=1e-3) +Test.@test isapprox(He[1,4] / He[1,3], 5; rtol=1e-3) @info("NMFk: nlopt: 2 sources, 5 sensors, 20 transients") Random.seed!(2015) @@ -157,39 +157,39 @@ W = [a b] H = [.1 1 0 0 .1; 0 0 .1 .5 .2] X = W * H X = [a a*10 b b*5 a+b*2] -@Suppressor.suppress global We, He, p, s = NMFk.execute(X, 2, 10; save=false, method=:nlopt, tolX=1e-6, tol=1e-19) -@Test.test isapprox(p, 0, atol=1e-3) -@Test.test isapprox(s, 1, rtol=1e-1) -@Test.test isapprox(He[1,2] / He[1,1], 10, rtol=1e-3) -@Test.test isapprox(He[2,2] / He[2,1], 10, rtol=1e-3) -@Test.test isapprox(He[2,4] / He[2,3], 5, rtol=1e-3) -@Test.test isapprox(He[1,4] / He[1,3], 5, rtol=1e-3) +Suppressor.@suppress global We, He, p, s = NMFk.execute(X, 2, 10; save=false, method=:nlopt, tolX=1e-6, tol=1e-19) +Test.@test isapprox(p, 0, atol=1e-3) +Test.@test isapprox(s, 1, rtol=1e-1) +Test.@test isapprox(He[1,2] / He[1,1], 10, rtol=1e-3) +Test.@test isapprox(He[2,2] / He[2,1], 10, rtol=1e-3) +Test.@test isapprox(He[2,4] / He[2,3], 5, rtol=1e-3) +Test.@test isapprox(He[1,4] / He[1,3], 5, rtol=1e-3) @info("NMFk: ipopt: 2 sources, 5 sensors, 21 transients") Random.seed!(2015) a = exp.(-(0:.5:10))*100 b = 100 .+ sin.(0:20)*10 X = [a a*10 b b*5 a+b*2] -@Suppressor.suppress global We, He, p, s = NMFk.execute(X, 2, 10; save=false, method=:ipopt, tolX=1e-3, tol=1e-7) -@Test.test isapprox(p, 0, atol=1e-3) -@Test.test isapprox(s, 1, rtol=1e-1) -@Test.test isapprox(He[1,2] / He[1,1], 10, rtol=1e-2) -@Test.test isapprox(He[2,3] / He[1,3], 0.38, rtol=0.1) -@Test.test isapprox(He[2,4] / He[2,3], 5, rtol=1e-3) -@Test.test isapprox(He[1,4] / He[1,3], 5, rtol=1e-3) +Suppressor.@suppress global We, He, p, s = NMFk.execute(X, 2, 10; save=false, method=:ipopt, tolX=1e-3, tol=1e-7) +Test.@test isapprox(p, 0, atol=1e-3) +Test.@test isapprox(s, 1, rtol=1e-1) +Test.@test isapprox(He[1,2] / He[1,1], 10, rtol=1e-2) +Test.@test isapprox(He[2,3] / He[1,3], 0.38, rtol=0.1) +Test.@test isapprox(He[2,4] / He[2,3], 5, rtol=1e-3) +Test.@test isapprox(He[1,4] / He[1,3], 5, rtol=1e-3) @info("NMFk: nlopt: 2 sources, 5 sensors, 100 transients") Random.seed!(2015) a = exp.(-(0:.5:10))*100 b = 100 .+ sin.(0:20)*10 X = [a a*10 b b*5 a+b*2] -@Suppressor.suppress global We, He, p, s = NMFk.execute(X, 2, 10; save=false, method=:nlopt, tolX=1e-12, tol=1e-19) -@Test.test isapprox(p, 0, atol=1e-3) -@Test.test isapprox(s, 1, rtol=1e-1) -@Test.test isapprox(He[1,2] / He[1,1], 10, rtol=1e-3) -@Test.test isapprox(He[1,3] / He[2,3], 3, rtol=0.9) -@Test.test isapprox(He[2,4] / He[2,3], 5, rtol=1e-3) -@Test.test isapprox(He[1,4] / He[1,3], 5, rtol=1e-3) +Suppressor.@suppress global We, He, p, s = NMFk.execute(X, 2, 10; save=false, method=:nlopt, tolX=1e-12, tol=1e-19) +Test.@test isapprox(p, 0, atol=1e-3) +Test.@test isapprox(s, 1, rtol=1e-1) +Test.@test isapprox(He[1,2] / He[1,1], 10, rtol=1e-3) +Test.@test isapprox(He[1,3] / He[2,3], 3, rtol=0.9) +Test.@test isapprox(He[2,4] / He[2,3], 5, rtol=1e-3) +Test.@test isapprox(He[1,4] / He[1,3], 5, rtol=1e-3) @info("NMFk: 3 sources, 5 sensors, 15 transients") Random.seed!(2015) @@ -198,42 +198,42 @@ b = rand(15) c = rand(15) X = [a+c*3 a*10 b b*5+c a+b*2+c*5] @info("NMFk: ipopt ...") -@Suppressor.suppress global We, He, p, s = NMFk.execute(X, 2:4, 10; maxiter=100, tol=1e-2, tolX=1e-2, save=false, method=:ipopt) +Suppressor.@suppress global We, He, p, s = NMFk.execute(X, 2:4, 10; maxiter=100, tol=1e-2, tolX=1e-2, save=false, method=:ipopt) @info("NMFk: nlopt ...") -@Suppressor.suppress global We, He, p, s = NMFk.execute(X, 2:4, 10; maxiter=100, tol=1e-2, tolX=1e-2, save=false, method=:nlopt) +Suppressor.@suppress global We, He, p, s = NMFk.execute(X, 2:4, 10; maxiter=100, tol=1e-2, tolX=1e-2, save=false, method=:nlopt) @info("NMFk: simple ...") -@Suppressor.suppress global We, He, p, s = NMFk.execute(X, 2:4, 10; maxiter=100, tol=1e-2, save=false, method=:simple) +Suppressor.@suppress global We, He, p, s = NMFk.execute(X, 2:4, 10; maxiter=100, tol=1e-2, save=false, method=:simple) @info("NMFk: nmf ...") -@Suppressor.suppress global We, He, p, s = NMFk.execute(X, 2:4, 10; maxiter=100, tol=1e-2, save=false, method=:nmf) +Suppressor.@suppress global We, He, p, s = NMFk.execute(X, 2:4, 10; maxiter=100, tol=1e-2, save=false, method=:nmf) @info("NMFk: with sparsity constraints ...") -@Suppressor.suppress global We, He, p, s = NMFk.execute(X, 2:4, 10; maxiter=100, tol=1e-2, save=false, method=:sparsity) +Suppressor.@suppress global We, He, p, s = NMFk.execute(X, 2:4, 10; maxiter=100, tol=1e-2, save=false, method=:sparsity) @info("NMFk: 3 sources, 7 sensors, 50 transients (no sparseness)") a = exp.(-(0.2:.2:10)) b = 1 .+ sin.(0:49) c = collect(0:49) / 49 X = [a b c a+c a+b b+c a+b+c] -@Suppressor.suppress global We, He, p, s = NMFk.execute(X, 2:4, 10; save=false, method=:simple) +Suppressor.@suppress global We, He, p, s = NMFk.execute(X, 2:4, 10; save=false, method=:simple) @info("NMFk: 3 sources, 7 sensors, 50 transients (with sparseness)") inan = rand(50, 7) .> .8 X[inan] .= NaN -@Suppressor.suppress global We, He, p, s = NMFk.execute(X, 2:4, 10; save=false, method=:simple) +Suppressor.@suppress global We, He, p, s = NMFk.execute(X, 2:4, 10; save=false, method=:simple) @info("NMFk: concentrations/delta tests ...") a0 = Float64[[20,10,1] [5,1,1]] b = NMFk.getisotopeconcentration(a0, [0.001,0.002], [[100,10,1] [500,50,5]]) a = NMFk.getisotopedelta(b, [0.001,0.002], [[100,10,1] [500,50,5]]) -@Test.test a0 ≈ a +Test.@test a0 ≈ a a0 = Float64[20,10,1] b = NMFk.getisotopeconcentration(a0, 0.001, [100,10,1]) a = NMFk.getisotopedelta(b, 0.001, [100,10,1]) -@Test.test a0 ≈ a +Test.@test a0 ≈ a a0 = 20.0 b = NMFk.getisotopeconcentration(a0, 0.001, 100) a = NMFk.getisotopedelta(b, 0.001, 100)[1] -@Test.test a0 ≈ a +Test.@test a0 ≈ a end