text:
using Noise
using Test
using Statistics
using Random
using ImageCore
using ColorVectorSpace

@testset "Util functions" begin
    x = 1f0
    @test 1f0 + 1f0im == Noise.complex_copy(x)

    x = 12f0 + 13f0im
    @test 12f0 + 13f0im == Noise.complex_copy(x)
end

Random.seed!(42)
tpl = (1000, 1000)
arr = rand(Float64, tpl)
arr_ones = ones(tpl)
arr_rand_int = convert(Array{Int64}, map(x -> round(1000*x), arr))
arr_zeros = zeros(tpl)
img_zeros = colorview(RGB, zeros(Normed{UInt8, 8}, 3, 1000, 1000))
img_zeros_gray = colorview(Gray, zeros(Normed{UInt8, 8}, 1000, 1000))

img_rand = colorview(RGB, rand(Normed{UInt8, 8}, 3, 1000, 1000))
img_rand_gray = colorview(Gray, rand(Normed{UInt8, 8}, 1000, 1000))

img_12 = colorview(RGB, 0.5 .* ones(Normed{UInt8, 8}, 3, 1000, 1000))
img_12_gray= colorview(Gray, 0.5 .* ones(Normed{UInt8, 8}, 1000, 1000))

img = colorview(RGB, ones(Normed{UInt8, 8}, 3, 1000, 1000))
img_gray= colorview(Gray,  ones(Normed{UInt8, 8}, 1000, 1000))

img_zeros_float = colorview(RGB, zeros(Float64, 3, 1000, 1000))
img_zeros_gray_float = colorview(Gray, zeros(Float64, 1000, 1000))

img_rand_float = colorview(RGB, rand(Float64, 3, 1000, 1000))
img_rand_gray_float = colorview(Gray, rand(Float64, 1000, 1000))

img_12_float = colorview(RGB, 0.5 .* ones(Float64, 3, 1000, 1000))
img_12_gray_float = colorview(Gray, 0.5 .* ones(Float64, 1000, 1000))

img_float = colorview(RGB, ones(Float64, 3, 1000, 1000))
img_gray_float = colorview(Gray,  ones(Float64, 1000, 1000))

@testset "Clip functions" begin

    @test mean(Noise.clip_v(1.1)) ≈ 1.0
    @test mean(Noise.clip_v(1.0)) ≈ 1.0
    @test mean(Noise.clip_v(10.0)) ≈ 1.0
    @test mean(Noise.clip_v(0.1)) ≈ 0.1
    @test mean(Noise.clip_v(-1.0)) ≈ 0.0
    @test mean(Noise.clip_v(0.0)) ≈ 0.0

end

@testset "Salt and Pepper noise" begin
    #check 100% salt
    @test salt_pepper(arr, 1.0, salt_prob=1.0) == ones(tpl)
    #check 100% pepper
    @test salt_pepper(arr, 1.0, salt_prob=0.0) == zeros(tpl)
    #check 100% salt with custom value
    @test salt_pepper(arr, 1.0, salt_prob=1.0, salt=0.23) == 0.23 .* ones(tpl)
    #check 100% pepper with custom value
    @test salt_pepper(arr, 1.0, salt_prob=0.0, pepper=0.63) == 0.63 .* ones(tpl)
    #check average salt and pepper roughly
    @test abs(mean(salt_pepper(arr, 1.0) .-0.5)) < 0.1
    #check that average salt and pepper are roughly 0.5
    @test abs(mean(salt_pepper(arr, 0.1) .-0.5)) < 0.05

    @test abs(mean(salt_pepper(arr, 0.1, salt_prob=1.0) .-0.5)) > 0.001
    @test abs(mean(salt_pepper(arr, 0.1, salt_prob=0.0) .-0.5)) > 0.001

    #check that average is different to 0.5 if salt or pepper are not 1.0 or 0.0
    @test abs(mean(salt_pepper(arr, 0.1, pepper=1.0) .-0.5)) > 0.001
    @test abs(mean(salt_pepper(arr, 0.1, salt=0.0) .-0.5)) > 0.001

    # check images chn
    @test abs(mean(channelview(salt_pepper_chn(img_zeros, 1.0))) - 0.5)< 0.001
    #check that average is different to 0.5 if salt or pepper are not 1.0 or 0.0
    @test abs(sum(channelview(salt_pepper_chn(img_zeros, 0.1, pepper=1.0)) .-0.5)/100/100) > 0.001
    @test abs(sum(channelview(salt_pepper_chn(img_rand, 0.1, salt=0.0)) .-0.5)/100/100) > 0.001

    # check RGB images
    @test abs(mean(channelview(salt_pepper(img_zeros, 1.0))) - 0.5)< 0.001
    #check that average is different to 0.5 if salt or pepper are not 1.0 or 0.0
    @test abs(mean(channelview(salt_pepper(img_zeros, 0.1, pepper=1.0)) .-0.5)) > 0.001
    @test abs(mean(channelview(salt_pepper(img_rand, 0.1, salt=0.0)) .-0.5)) > 0.001
    # check std
    @test abs(std(channelview(salt_pepper(img_rand, 1.0, salt=0.0, pepper=0.5))) -0.25) < 0.001
    @test abs(std(channelview(salt_pepper(img_rand, 1.0, salt=0.4, pepper=0.5))) -0.05) < 0.001

    # check gray images
    @test abs(mean(channelview(salt_pepper(img_zeros_gray, 1.0))) - 0.5)< 0.001
    #check that average is different to 0.5 if salt or pepper are not 1.0 or 0.0
    @test abs(mean(channelview(salt_pepper(img_zeros_gray, 0.1, pepper=1.0)) .-0.5)) > 0.001
    @test abs(mean(channelview(salt_pepper(img_rand_gray, 0.1, salt=0.0)) .-0.5)) > 0.001
    # check std
    @test abs(std(channelview(salt_pepper(img_rand_gray, 1.0, salt=0.0, pepper=0.5))) -0.25) < 0.001
    @test abs(std(channelview(salt_pepper(img_rand_gray, 1.0, salt=0.4, pepper=0.5))) -0.05) < 0.001

    # check images but with float
    @test abs(mean(channelview(salt_pepper_chn(img_zeros_float, 1.0))) - 0.5)< 0.001
    #check that average is different to 0.5 if salt or pepper are not 1.0 or 0.0
    @test abs(sum(channelview(salt_pepper_chn(img_zeros_float, 0.1, pepper=1.0)) .-0.5)/100/100) > 0.001
    @test abs(sum(channelview(salt_pepper_chn(img_rand_float, 0.1, salt=0.0)) .-0.5)/100/100) > 0.001

    # check RGB images
    @test abs(mean(channelview(salt_pepper(img_zeros_float, 1.0))) - 0.5)< 0.001
    #check that average is different to 0.5 if salt or pepper are not 1.0 or 0.0
    @test abs(mean(channelview(salt_pepper(img_zeros_float, 0.1, pepper=1.0)) .-0.5)) > 0.001
    @test abs(mean(channelview(salt_pepper(img_rand_float, 0.1, salt=0.0)) .-0.5)) > 0.001
    # check std
    @test abs(std(channelview(salt_pepper(img_rand_float, 1.0, salt=0.0, pepper=0.5))) -0.25) < 0.001
    @test abs(std(channelview(salt_pepper(img_rand_float, 1.0, salt=0.4, pepper=0.5))) -0.05) < 0.001

    # check gray images
    @test abs(mean(channelview(salt_pepper(img_zeros_gray_float, 1.0))) - 0.5)< 0.001
    #check that average is different to 0.5 if salt or pepper are not 1.0 or 0.0
    @test abs(mean(channelview(salt_pepper(img_zeros_gray_float, 0.1, pepper=1.0)) .-0.5)) > 0.001
    @test abs(mean(channelview(salt_pepper(img_rand_gray_float, 0.1, salt=0.0)) .-0.5)) > 0.001
    # check std
    @test abs(std(channelview(salt_pepper(img_rand_gray_float, 1.0, salt=0.0, pepper=0.5))) -0.25) < 0.001
    @test abs(std(channelview(salt_pepper(img_rand_gray_float, 1.0, salt=0.4, pepper=0.5))) -0.05) < 0.001

end

@testset "Additive white Gaussian" begin

    @test (std(imag(add_gauss(Complex.(arr_zeros), 13im))) -13) < 0.1
    @test (std(real(add_gauss(Complex.(arr_zeros), 13im)))) == 0
    @test (std(imag(add_gauss(Complex.(arr_zeros), 13))) -13) < 0.1
    @test (std(real(add_gauss(Complex.(arr_zeros), 13))) -13) < 0.1

    # check array with white gaussian noise
    @test abs(std(add_gauss(arr_zeros, 13.0)) - 13.0) < 0.1
    @test abs(mean(add_gauss(arr_zeros, 1.0))) < 0.05
    @test abs(mean(add_gauss(arr_zeros, 1.0, 10.0)) - 10) < 0.05
    @test mean(add_gauss(arr_zeros, 0.5, -1.0, clip=false)) < 0
    @test mean(add_gauss(arr_zeros, 0.5, -1.0, clip=true)) >= 0
    @test mean(add_gauss(arr_zeros, 0.5, 2.0, clip=false)) > 1
    @test mean(add_gauss(arr_zeros, 0.5, 2.0, clip=true)) <= 1

    # check images for gaussian white noise
    # check mean offset channelwise
    @test abs(mean(channelview(add_gauss_chn(img_zeros, 0.1, 0.5))) - 0.5) < 0.05
    @test abs(mean(channelview(add_gauss_chn(img_zeros, 0.2, 0.3))) - 0.3) < 0.05

    # check mean offset with clip
    @test abs(mean(channelview(add_gauss(img_zeros, 0.1, 0.5))) - 0.5) < 0.05
    @test abs(mean(channelview(add_gauss(img_zeros, 0.2, 0.3))) - 0.3) < 0.05
    @test abs(mean(channelview(add_gauss(img_zeros, 0.0, 10))) - 1.0) < 0.005

    # check mean offet with gray
    @test abs(mean(channelview(add_gauss(img_zeros_gray, 0.1, 0.5))) - 0.5) < 0.05
    @test abs(mean(channelview(add_gauss(img_zeros_gray, 0.2, 0.3))) - 0.3) < 0.05
    @test abs(mean(channelview(add_gauss(img_zeros_gray, 0.0, 10))) - 1.0) < 0.005

    # check the same but with float images
    # check images for gaussian white noise
    # check mean offset channelwise
    @test abs(mean(channelview(add_gauss_chn(img_zeros_float, 0.1, 0.5))) - 0.5) < 0.05
    @test abs(mean(channelview(add_gauss_chn(img_zeros_float, 0.2, 0.3))) - 0.3) < 0.05

    # check mean offset, now it's not clipped
    @test abs(mean(channelview(add_gauss(img_zeros_float, 0.1, 0.5))) - 0.5) < 0.05
    @test abs(mean(channelview(add_gauss(img_zeros_float, 0.2, 0.3))) - 0.3) < 0.05
    @test abs(mean(channelview(add_gauss(img_zeros_float, 0.0, 10))) - 10.0) < 0.005

    # check mean offet with gray
    @test abs(mean(channelview(add_gauss(img_zeros_gray_float, 0.1, 0.5))) - 0.5) < 0.05
    @test abs(mean(channelview(add_gauss(img_zeros_gray_float, 0.2, 0.3))) - 0.3) < 0.05
    @test abs(mean(channelview(add_gauss(img_zeros_gray_float, 0.0, 10))) - 10.0) < 0.005

end

@testset "Poisson Noise" begin
    #test std of poisson with clip and no clip
    @test abs(std(poisson(arr_ones .* 13.0)) - sqrt(13)) < 0.01
    @test abs(std(poisson(arr_ones .* 13.0, clip=true))) < 0.01
    #test mean with clip and no clip
    @test abs(mean(poisson(arr_ones .* 13.0)) - 13) < 0.05
    @test abs(mean(poisson(arr_ones .* 13.0, clip=true)) - 1) < 0.001
    #test mean with scaling
    @test abs(mean(poisson(arr, 1000.0))-0.5) < 0.005
    #test mean int
    @test abs(mean(poisson(arr_rand_int)))-500<1

    @test abs(std(channelview(poisson(img_12_gray, 100.0))) - 0.5 .* sqrt(100.0)./100.0) < 0.005
    @test abs(std(channelview(poisson(img_12, 100.0))) - 0.5 .* sqrt(100.0)./100.0) < 0.005
    @test abs(std(channelview(poisson(img_12_gray, 4200.0))) - 0.3 .* sqrt(4200.0)./4200.0) < 0.005
    @test abs(std(channelview(poisson(img_12, 4200.0))) - 0.3 .* sqrt(4200.0)./4200.0) < 0.005

    @test abs(std(channelview(poisson(img_gray, 10000.0, clip=true))))< 0.02
    @test abs(std(channelview(poisson(img, 10000.0, clip=true)))) < 0.02
    @test abs(mean(channelview(poisson(img_gray, 10000.0, clip=true)))-1.0)< 0.02
    @test abs(mean(channelview(poisson(img, 10000.0, clip=true))) -1.0) < 0.02

    @test abs(std(channelview(poisson(img_12_gray_float, 100.0))) - 0.5 .* sqrt(100.0)./100.0) < 0.005
    @test abs(std(channelview(poisson(img_12_float, 100.0))) - 0.5 .* sqrt(100.0)./100.0) < 0.005
    @test abs(std(channelview(poisson(img_12_gray_float, 4200.0))) - 0.3 .* sqrt(4200.0)./4200.0) < 0.005
    @test abs(std(channelview(poisson(img_12_float, 4200.0))) - 0.3 .* sqrt(4200.0)./4200.0) < 0.005

    @test abs(std(channelview(poisson(img_gray_float, 10000.0, clip=true))))< 0.02
    @test abs(std(channelview(poisson(img_float, 10000.0, clip=true)))) < 0.02
    @test abs(mean(channelview(poisson(img_gray_float, 10000.0, clip=true)))-1.0)< 0.02
    @test abs(mean(channelview(poisson(img_float, 10000.0, clip=true))) -1.0) < 0.02

end

@testset "Quantization" begin
    @test quantization(img_gray, 10000) ≈ img_gray
    @test quantization(img_gray, 255) ≈ img_gray
    @test abs(sum(quantization(img_gray, 255)) - sum(img_gray)) < 0.01 * sum(img_gray)

    @test quantization(img, 10000) ≈ img
    @test quantization(img, 255) ≈ img

    @test quantization(img_float, 1000) ≈ img_float
    @test quantization(img_float, 255) ≈ img_float

    @test abs((quantization(img_12_gray_float, 100) |> channelview |> mean) - 0.5) < 0.01
    @test abs((quantization(img_12_gray_float, 50) |> channelview |> mean) - 0.5) < 0.05
    @test abs(sum(quantization(img_12_gray_float, 255)) - sum(img_12_gray_float)) < 0.01 * sum(img_12_gray_float)

    @test abs((quantization(arr , 100) |> mean) - 0.5) < 0.01
    @test abs((quantization(arr, 50) |> mean) - 0.5) < 0.05
    @test abs(sum(quantization(arr, 255)) - sum(arr)) < 0.01 * sum(arr)

    @test minimum(quantization(img_rand_gray_float, 2, minv=0.32, maxv=0.98)) ≈ 0.32
    @test maximum(quantization(img_rand_gray_float, 2, minv=0.32, maxv=0.98)) ≈ 0.98

    @test minimum(quantization(img_rand_gray, 5, minv=0.1, maxv=0.96)) ≈ Gray{N0f8}(0.1)
    @test maximum(quantization(img_rand_gray, 2, minv=0.32, maxv=0.98)) ≈ Gray{N0f8}(0.98)
end

@testset "Multiplicative gauss Noise" begin

    @test (std(imag(mult_gauss(Complex.(arr_ones), 13im))) -13) < 0.1
    @test (std(real(mult_gauss(Complex.(arr_ones), 13im)))) == 0
    @test (std(imag(mult_gauss(Complex.(arr_ones), 13))) -13) < 0.1
    @test (std(real(mult_gauss(Complex.(arr_ones), 13))) -13) < 0.1

    # check array with white gauss noise
    @test abs(std(mult_gauss(arr_zeros, 13.0)) - 0) < 0.1
    @test mean(mult_gauss(arr_zeros, 0.5, -1.0, clip=false)) ≈ 0
    @test mean(mult_gauss(arr_zeros, 0.5, -1.0, clip=true)) ≈ 0

    # check images for gauss white noise
    # check mean offset channelwise
    @test abs(mean(channelview(mult_gauss_chn(img_zeros, 0.1, 0.5)))) ≈ 0
    @test abs(mean(channelview(mult_gauss_chn(img_zeros, 0.2, 0.3)))) ≈ 0

    # check mean offset with clip
    @test abs(mean(channelview(mult_gauss(img_zeros, 0.1, 0.5)))) ≈  0
    @test abs(mean(channelview(mult_gauss(img_zeros, 0.2, 0.3)))) ≈ 0
    @test abs(mean(channelview(mult_gauss(img_zeros, 0.0, 10)))) ≈ 0

    # check mean offset with clip
    @test abs(mean(channelview(mult_gauss(img, 0.1, 0.5))) - 0.5) < 0.05
    @test abs(std(channelview(mult_gauss(img, 0.2, 0.3))) - 0.2) < 0.05
    @test abs(mean(channelview(mult_gauss(img, 0.01, 0.2))) - 0.2) < 0.05
    @test abs(std(channelview(mult_gauss(img, 0.1, 0.5))) - 0.1) < 0.05

end

max_issues_repo_name: roflmaostc/Noise

max_forks_repo_name: roflmaostc/Noise


text:
using NLsolve

"Parameters for right-hand side function of projected Gauss-Legendre Runge-Kutta methods."
mutable struct ParametersVPRKpTableau{DT,TT,D,S,ET} <: Parameters{DT,TT}
    equ::ET
    tab::CoefficientsPGLRK{TT}
    Δt::TT

    t̅::TT
    t::TT

    q̅::Vector{DT}
    p̅::Vector{DT}

    λ::Vector{DT}
    A::Matrix{DT}
    W::Matrix{DT}

    function ParametersVPRKpTableau{DT,TT,D,S,ET}(equ, tab, Δt) where {DT,TT,D,S,ET}
        new(equ, tab, Δt, zero(TT), zero(TT), zeros(DT,D), zeros(DT,D), zeros(DT,D), zeros(DT,S,S), zeros(DT,S,S))
    end
end

"""
Projected Variational Gauss-Legendre Runge-Kutta integrator.
"""
struct IntegratorVPRKpTableau{DT, TT, PT <: ParametersVPRKpTableau{DT,TT},
                               ST <: NonlinearSolver{DT},
                               IT <: InitialGuessPODE{DT,TT}, N, D, S} <: IntegratorPRK{DT,TT}
    params::PT
    solver::ST
    iguess::IT
    cache::IntegratorCacheVPRK{DT,D,S}

    function IntegratorVPRKpTableau(N, params::ParametersVPRKpTableau{DT,TT,D,S,ET}, solver::ST, iguess::IT) where {DT, TT, D, S, ET, ST, IT}
        # create cache
        cache = IntegratorCacheVPRK{DT,D,S}()

        # create integrator
        new{DT, TT, typeof(params), ST, IT, N, D, S}(params, solver, iguess, cache)
    end
end

function IntegratorVPRKpTableau(equation::IODE{DT,TT,ΑT,FT,GT,HT,VT,N}, tableau::CoefficientsPGLRK{TT}, Δt::TT) where {DT,TT,ΑT,FT,GT,HT,VT,N}
    D = equation.d
    M = equation.n
    S = tableau.s

    # create params
    params = ParametersVPRKpTableau{DT,TT,D,S,typeof(equation)}(equation, tableau, Δt)

    # create solver
    solver  = create_nonlinear_solver(DT, D*S, params)

    # create initial guess
    iguess = InitialGuessPODE(get_config(:ig_interpolation), equation, Δt)

    # create integrator
    IntegratorVPRKpTableau(N, params, solver, iguess)
end

@inline nstages(integrator::IntegratorVPRKpTableau{DT,TT,PT,ST,IT,N,D,S}) where {DT,TT,PT,ST,IT,N,D,S} = S

function update_params!(params::ParametersVPRKpTableau, sol::AtomicSolutionPODE)
    # set time for nonlinear solver and copy previous solution
    params.t̅  = sol.t
    params.t  = sol.t + params.Δt
    params.q̅ .= sol.q
    params.p̅ .= sol.p
end

function update_tableau!(params::ParametersVPRKpTableau{DT,TT,D,S}, λ::Vector) where {DT,TT,D,S}
    # copy λ to parameters
    params.λ .= λ

    # compute tableaus
    for k in eachindex(λ)
        params.W[S-k+1, S-k] = +λ[k]
        params.W[S-k, S-k+1] = -λ[k]
    end

    # TODO make this more allocation-efficient
    params.A .= params.tab.P * params.W * params.tab.Q
    # simd_mult!(params.T, params.W, params.tab.Q)
    # simd_mult!(params.A, params.tab.P, params.T)
end

function compute_stages!(x::Vector{ST}, cache::IntegratorCacheVPRK{ST}, params::ParametersVPRKpTableau{DT,TT,D,S}) where {ST,DT,TT,D,S}
    compute_stages!(x, cache.Q, cache.V,
                       cache.P, cache.F,
                       cache.Y, cache.Z,
                       cache.q̃, cache.p̃, cache.θ̃,
                       cache.ṽ, cache.f̃,
                       params)
end

function compute_stages!(x::Vector{ST}, Q::Vector{Vector{ST}}, V::Vector{Vector{ST}},
                                        P::Vector{Vector{ST}}, F::Vector{Vector{ST}},
                                        Y::Vector{Vector{ST}}, Z::Vector{Vector{ST}},
                                        q::Vector{ST}, p::Vector{ST}, θ::Vector{ST},
                                        y::Vector{ST}, z::Vector{ST},
                                        params::ParametersVPRKpTableau{DT,TT,D,S}) where {ST,DT,TT,D,S}

    local tᵢ::TT
    local y₁::ST
    local y₂::ST
    local z₁::ST
    local z₂::ST

    # copy x to Y
    for i in 1:S
        for k in 1:D
            V[i][k] = x[D*(i-1)+k]
        end
    end

    # compute Y
    for i in 1:S
        for k in 1:D
            y₁ = 0
            y₂ = 0
            for j in 1:S
                y₁ += params.tab.a[i,j] * V[j][k]
                y₂ += params.A[i,j] * V[j][k]
            end
            Y[i][k] = y₁ + y₂
        end
    end

    # compute Q=q̅+Δt*Y
    for i in 1:S
        for k in 1:D
            Q[i][k] = params.q̅[k] + params.Δt * Y[i][k]
        end
    end

    # compute P=ϑ(Q,V) and F=f(Q,V)
    for i in 1:S
        tᵢ = params.t̅ + params.Δt * params.tab.c[i]
        params.equ.ϑ(tᵢ, Q[i], V[i], P[i])
        params.equ.f(tᵢ, Q[i], V[i], F[i])
    end

    # compute Z
    for i in 1:S
        for k in 1:D
            z₁ = 0
            z₂ = 0
            for j in 1:S
                z₁ += params.tab.a[i,j] * F[j][k]
                z₂ += params.A[i,j] * F[j][k]
            end
            Z[i][k] = z₁ + z₂
        end
    end

    # compute y=B*V and z=B*F
    y .= 0
    z .= 0
    for k in 1:D
        for j in 1:S
            y[k] += params.tab.b[j] * V[j][k]
            z[k] += params.tab.b[j] * F[j][k]
        end
    end

    # compute q=q̅+Δt*y and p=p̅+Δt*z
    q .= params.q̅ .+ params.Δt .* y
    p .= params.p̅ .+ params.Δt .* z

    # compute θ=ϑ(t,q)
    params.equ.ϑ(params.t, q, θ)
end

"Compute stages of projected Gauss-Legendre Runge-Kutta methods."
@generated function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersVPRKpTableau{DT,TT,D,S}) where {ST,DT,TT,D,S}

    cache = IntegratorCacheVPRK{ST,D,S}()

    quote
        compute_stages!(x, $cache, params)

        # compute b = [P-p-AF]
        for i in 1:S
            for k in 1:D
                b[D*(i-1)+k] = $cache.P[i][k] - params.p̅[k] - params.Δt * $cache.Z[i][k]
            end
        end
    end
end

function function_dirac_constraint!(λ::Vector, int::IntegratorVPRKpTableau{DT,TT}) where {DT,TT}
    # copy λ to integrator parameters
    update_tableau!(int.params, λ)

    # println(λ)

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector fields at internal stages
    compute_stages!(int.solver.x, int.cache, int.params)

    # compute and return ϑ(t,q)-p
    return int.cache.θ̃ .- int.cache.p̃
end

function initialize!(int::IntegratorVPRKpTableau, sol::AtomicSolutionPODE)
    sol.t̅ = sol.t - timestep(int)

    equation(int).v(sol.t, sol.q, sol.p, sol.v)
    equation(int).f(sol.t, sol.q, sol.p, sol.f)

    initialize!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f,
                            sol.t̅, sol.q̅, sol.p̅, sol.v̅, sol.f̅)
end

function initial_guess!(int::IntegratorVPRKpTableau, sol::AtomicSolutionPODE)
    for i in eachstage(int)
        evaluate!(int.iguess, sol.q, sol.p, sol.v, sol.f,
                              sol.q̅, sol.p̅, sol.v̅, sol.f̅,
                              int.cache.q̃, int.cache.ṽ,
                              tableau(int).c[i])

        for k in eachdim(int)
            int.solver.x[ndims(int)*(i-1)+k] = int.cache.ṽ[k]
        end
    end
end

"Integrate ODE with projected Gauss-Legendre Runge-Kutta integrator."
function integrate_step!(int::IntegratorVPRKpTableau{DT,TT}, sol::AtomicSolutionPODE{DT,TT}) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int.params, sol)

    # compute initial guess
    initial_guess!(int, sol)

    # reset solution
    reset!(sol, timestep(int))

    # determine parameter λ
    nlres = nlsolve(λ -> function_dirac_constraint!(λ, int), zero(int.params.λ);
                xtol=get_config(:nls_atol),
                ftol=maximum(int.params.p̅)*get_config(:nls_atol),
                iterations=100)
                #xtol=timestep(int)^nstages(int)*get_config(:nls_atol),

    int.params.λ .= nlres.zero

    # println("x converged = ", nlres.x_converged, ", ",
    #         "f converged = ", nlres.f_converged, ", ",
    #         "residual norm = ", nlres.residual_norm, ", ",
    #         "iterations = ", nlres.iterations)
    #
    # println("λ = ", int.params.λ)
    # println("Δ = ", int.cache.θ̃ .- int.cache.p̃)
    # println("p = ", int.cache.p̃)
    # println()

    # compute final update
    update_solution!(sol.q, int.cache.V, tableau(int).b, timestep(int))
    update_solution!(sol.p, int.cache.F, tableau(int).b, timestep(int))

    # copy solution to initial guess
    update!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f)
end

ext: jl


text:
using ModelingToolkit
using DifferentialEquations

@parameters t α
@variables x(t) y(t)
D = Differential(t)

eqs = [D(x) ~ y * x,
       D(y) ~ -α * x]

@named sys = ODESystem(eqs)

u0 = [ x => 1.0, y => 2.0]
tspan = (0.0, 1.0)
p = [α => 1.0]
prob = ODEProblem(sys, u0, tspan, p)
sol = solve(prob)

#=
using SymbolicUtils
r1 = @rule (~y)*(~x) => (~y)^2
r1(eqs[1].rhs)
=#

max_stars_repo_name: SteffenPL/Julia-for-mathematicians

max_issues_repo_name: SteffenPL/Julia-for-mathematicians

max_forks_repo_name: SteffenPL/Julia-for-mathematicians


text:
##------------------ Binary UNDIRECTED Networks

struct  ErgmBin1 <: ErgmBin #Bin stands for (Binary) Adjacency matrix
    "ERGM with degrees as statistics and possibility to have groups of nodes
    associated with the same parameter. "
    obs:: Array{<:Real,1} #[degrees]
    Par::Array{<:Real,1} # One parameter per Node  UNRestricted
    groupsInds::Array{<:Real,1} #The index of the groups each node belongs to
 end
fooErgmBin1 =  ErgmBin1(ones(3),zeros(3),ones(3))
ErgmBin1(deg:: Array{<:Real,1}) =  ErgmBin1(deg,zeros(length(deg)),Vector(1:length(deg)) )
function DegSeq2graphDegSeq(Model::ErgmBin1,Seq::Array{<:Real,1})
    "From a degree sequence find another one that is certanly graphicable.
     I choose a naive way to ensure graphicability.
      The matrix associated with the starting degree sequence is very hubs oriented
      the largest node catches as many links as it needs to fill all its strubs "

    Seq = sort(Seq,rev=true)#decreasing sort
    N = length(Seq)
    stubs = ceil.(Int,copy(Seq))
    mat = zeros(Int64,N,N)
    for i in N:-1:1
       j = 1
       while (stubs[i]>0)&(j<=N)
           if (mat[i,j] == 0)& (i!=j)& (stubs[j]>0)
               mat[i,j] = 1
               mat[j,i] = 1
               stubs[i] -= 1
               stubs[j] -= 1
               #println(j)
           else
               #println((i, j,mat[i,j],stubs[j]))
           end
           j += 1
       end
    end
    #display(stubs)
    #display(mat)
    ((sum(mat[Array{Bool}(eye(mat))])) != 0)  ?  error() : ()
    graphSeq= sumSq(mat,2)
    return graphSeq
 end
function linSpacedPar(Model::ErgmBin1,Nnodes::Int; Ngroups = Nnodes,deltaN::Int=3,graphConstr = true)
    "Generate an equally spaced degree sequence. Ensure that it is graphicable
    and then estimate the parameters "

          #generate syntetic degree distributions, equally spaced, and estimate the model
    if graphConstr
        tmpdegs  = DegSeq2graphDegSeq(fooErgmBin1, Array{Real,1}(range(1+deltaN,stop=Nnodes - round(0.1*Nnodes)-deltaN, length=Nnodes)))
    else
        #stima modello a Ngroups fitness costanti con degree dei nodi linearly spaced
        tmpdegs  =   Array{Real,1}(range(deltaN,stop=Nnodes-deltaN, length=Nnodes) )
        tmpdegs = sort(tmpdegs,rev =true)
    end
    #if Nnodes==Ngroups
    groupsInds = distributeAinVecN(Vector(1:Ngroups),Nnodes)
    out,it,estMod = estimate(ErgmBin1(tmpdegs, zeros(length(tmpdegs)), groupsInds))
    #else
    return out,tmpdegs
 end
function samplSingMatCan(Model::ErgmBin1,expMat::Array{<:Real,2})
    """
    Being p a matrix of entries in [0,1], samples one bernoulli r.v. for eache matrix element
    """
    N1,N2 = size(expMat)
    N1!=N2  ?  error() : N=N1
    out = zeros(Int8,N,N)
    for c = 1:N
        for r=1:c-1
            out[r,c] = rand(Bernoulli(expMat[r,c]))
        end
    end
    out = Symmetric(out)
    return out
 end

function expMatrix2(Model::ErgmBin1,par::Array{<:Real,1})
    "Given the vector of model parameters, return the product matrix(often useful
    in likelihood computation) and the expected matrix"
    parMat = putZeroDiag( exp.(Symmetric(par  .+ par') ))
    expMat = parMat ./ (1 .+ parMat)
    return parMat, expMat
 end
expMatrix(Model::ErgmBin1,par::Array{<:Real,1}) = expMatrix2(Model,par)[2] # return only the expected matrix given the parameters
function expValStatsFromMat(Model::ErgmBin1,expMat::Array{<:Real,2}  )
    "Expected value of the statistic for a model, given the expected matrix"
    expDeg = sumSq(expMat,2)
    expValStats = expDeg
    return expValStats
 end
function firstOrderCond(Model::ErgmBin1;deg::Array{<:Real,1} = Model.obs, par::Array{<:Real,1}=Model.Par, groupsInds::Array{<:Real,1} = Model.groupsInds )
    "Given the model, the degree sequence, the parameters and groups assignements,
    return the First order conditions. Gradient of the loglikelihood, or system
    of differences between expected and observed values of the statistics (either
    degrees of single nodes or degree of a group). "
    N = length(deg)
    groupsNames = unique(groupsInds)
    NG = length(groupsNames)

    if N==NG
        expMat = expMatrix(Model,par)
        expDeg = sumSq(expMat,2)
        foc = deg - expDeg
    elseif N<NG
        error()
    else
        tmpDic = countmap(groupsInds)
        # cumulative degree of each group
        degGroups = zeros(NG) ; for i =1:N degGroups[groupsInds[i]] += deg[i] end
        nodesPar = par[groupsInds]
        expMat = expMatrix(Model,nodesPar)
        #display(expMat)
        expDeg = sumSq(expMat,2)
        expDegGroups = zeros(NG) ; for i =1:N expDegGroups[groupsInds[i]] += expDeg[i] end
        foc = degGroups - expDegGroups
    end
    return foc
 end
function identify(Model::ErgmBin1,Par::Array{<:Real,1})
    "Given a vector of parameters, return the transformed vector that verifies
    an identification condition. Do not do anything if the model is identified."
    return Par
  end
function estimate(Model::ErgmBin1; deg::Array{<:Real,1} = Model.obs,
            groupsInds = Model.groupsInds, targetErr::Real=targetErrValStaticNets,bigConst = bigConstVal)
    "Given model type, observations of the statistics and groups assignments,
    estimate the parameters."
    # The parameters inside this function are Restricted to be in [0,1] hence
    # they need to be converted (log(par)) in order to be used outside
    N = length(deg)
    groupsNames = unique(groupsInds)
    NG = length(groupsNames)

    tmpDic = countmap(groupsInds)
    #number of nodes per group
    NgroupsMemb = sortrows([[i for i in keys(tmpDic)] [  tmpDic[i] for i in keys(tmpDic) ] ])
    # agrregate degree of each group
    degGroups = zeros(NG) ; for i =1:N degGroups[groupsInds[i]] += deg[i] end

    #sort the groups for degrees
    indsDegGroupsOrder = Int.(sortrows([degGroups Vector(1:NG)], rev = true)[:,2])

    Nzer = sum(degGroups.==0)
    L = sum(deg)
    LperLink  = L/(N*N-N)
    #println((LperLink,N,L))
    #uniform parameter as starting value
    unifstartval =  sqrt((LperLink/((1-LperLink))));
    #analytically solvable case
    if NG==1
        outParVec = unifstartval*ones(NG)
        outMod = ErgmBin1(Model.obs,outParVec,groupsInds)
        i = 0
    else
        par = unifstartval*ones(NG)
        eps = 1e-8
        i=0
        maxIt = 150
        relErr = 1
        tol_x = 10*Base.eps()
        function tmpFindZero(objFun::Function,brak::Array{<:Real,1},start::Real)
        #    wrapper find_zero options
            #println((brak[1],brak[2]))
            #println((objFun(brak[1]),objFun(brak[2])))
            out = find_zero(objFun,(brak[1],brak[2]),Bisection())
            #println((out,objFun(out)))
            return out
        end
        while (i<maxIt)&(relErr>targetErr)
            parMax = maximum(par)
            parMin = minimum(par[par.!=0])
            for nInd=1:NG
                #solve first the gradient of nodes with largest degrees
                ng = indsDegGroupsOrder[nInd]
                #put to zero each component of the gradient1
                function objfun_ng(par_ng::Real)
                     parVecOthGroups = par_ng.*par[1:end .!=ng]
                     expGroupDeg = NgroupsMemb[ng,2]*( sum( NgroupsMemb[1:end .!=ng ,2].*( parVecOthGroups ./ (1 + parVecOthGroups)) )    + (NgroupsMemb[ng,2]-1)* (par_ng^2)/(1+(par_ng^2)) )
                     return  degGroups[ng] - expGroupDeg
                end
                #find zero of given component

                if degGroups[ng]== (N-1-Nzer)*NgroupsMemb[ng,2]

                    par[ng] = (1/parMin)*(1-0.1*targetErr)/( 0.1*targetErr)
                elseif degGroups[ng]==0
                    par[ng] = 0
                else
                     tmpBnd = (degGroups[ng]/(NgroupsMemb[ng,2]*(N-1-Nzer)) )/(1-degGroups[ng]/(NgroupsMemb[ng,2]*(N-1-Nzer)) )

                    brakets = [ 0,2* tmpBnd/parMin]
                    par[ng] = tmpFindZero(objfun_ng,brakets,par[ng])
                end
                #println(par)
                #parUN = log.(par)
                #display(firstOrderCond(Model;deg = deg,par = parUN, groupsInds = groupsInds )')

            end
            # The parameters inside this function are Restricted to be in [0,1] hence
            # they need to be converted (log(par)) in order to be used outside
            parUN = log.(par)
            parUN[par.==0] = -bigConstVal
            err = firstOrderCond(Model;deg = deg,par = parUN, groupsInds = groupsInds )
            #    tmpMat = par.*par';expMat = putZeroDiag!(tmpMat ./ (1+tmpMat)) ;expDeg = sumSq(expMat,2)
            relErrVec = abs.(err ./ degGroups)
            relErr =  maximum(relErrVec[degGroups.!=0])
            # display(err')
             #println(relErr)
            i+=1
         end
         if i<maxIt
             # do nothing, optimization succeded
         else
             println(deg)
             if NG==N
                 error("with one node per group the optimization should succeed")
             else
                 error("Add a maximum likelihood refinement for the groups case if the map does not succeeed, as it might very well happen")
             end
         end
        parUN = log.(par)
        parUN[par.==0] = -bigConst
        outParVec = parUN
        outMod = ErgmBin1(Model.obs,outParVec,groupsInds)
    end
    return outParVec, i , outMod
 end
function sampl(Model::ErgmBin1,Nsample::Int;  par::Array{<:Real,1}=Model.Par ,groupsInds = Model.groupsInds)
    "Sample a given model."
    nodesPar = par[groupsInds]
    N = length(nodesPar)
    expMat = expMatrix(Model,nodesPar)
    SampleMats =zeros(Int8,N,N,Nsample)
    for s=1:Nsample
        SampleMats[:,:,s] = samplSingMatCan(Model,expMat)
    end
    return SampleMats
 end


text:
### Motion/Measurement Models ###

"""
`DynamicsModel` is an abstract type to encapsulate linear
and nonlinear dynamics models
"""
abstract type DynamicsModel end

"""
`ObservationModel` is an abstract type to encapsulate linear
and nonlinear observation models
"""
abstract type ObservationModel end

"""
    LinearDynamicsModel(A::AbstractMatrix,B::AbstractMatrix,W::Symmetric)
    LinearDynamicsModel(A::AbstractMatrix,B::AbstractMatrix,W::AbstractMatrix)

Construct linear dynamics model with transition matrix A,
control matrix B, and symmetric zero-mean process noise with
symmetric covariance matrix W
"""
struct LinearDynamicsModel{a<:AbstractMatrix, b<:AbstractMatrix,
                c<:Symmetric} <: DynamicsModel
    A::a
    B::b
    W::c
    #=
    function LinearDynamicsModel{a,b,c}(A::Matrix{a},B::Matrix{b},W::Symmetrix{c})
        where{a<:Number,b<:Number,c<:Number}

        (size(A,1) == size(B,1) == size(W,1)) ?
        new(A,B,W) : error("first dimensions not matching")
    end
    =#
end

function LinearDynamicsModel(A::AbstractMatrix, B::AbstractMatrix,
                             W::AbstractMatrix)

    return LinearDynamicsModel(A,B,Symmetric(W))
end

"""
    predict(m::LinearDynamicsModel, x::AbstractVector{<:Number}, u::AbstractVector{<:Number})
    predict(m::LinearDynamicsModel, x::AbstractVector{<:Number}, u::AbstractVector{<:Number}, rng::AbstractRNG)

Uses the linear dynamics model to propagate the state x one step forward in time with control input u.
If rng is given, it adds process noise.
"""
function predict(m::LinearDynamicsModel, x::AbstractVector{<:Number},
                 u::AbstractVector{<:Number})
    return m.A * x + m.B * u
end
function predict(m::LinearDynamicsModel, x::AbstractVector{T},
                 u::AbstractVector{T},
                 rng::AbstractRNG) where T<:Number
    return m.A * x + m.B * u + cholesky(m.W).L * randn(rng, T, size(m.W, 1))
end

function jacobian(m::LinearDynamicsModel, x, u)
    return m.A
end

"""
    LinearObservationModel(C::AbstractMatrix,D::AbstractMatrix,V::Symmetric)
    LinearObservationModel(C::AbstractMatrix,D::AbstractMatrix,V::AbstractMatrix)
    LinearObservationModel(C::AbstractMatrix,V::Symmetric)
    LinearObservationModel(C::AbstractMatrix,V::AbstractMatrix)

Construct linear observation dynamics model with output matrix C,
feedthrough matrix D, and symmetric zero-mean measurement noise with
symmetric covariance matrix V
"""
struct LinearObservationModel{a<:AbstractMatrix,b<:AbstractMatrix,c<:Symmetric} <: ObservationModel
    C::a
    D::b
    V::c
    #=
    LinearObservationModel(C,D,V) = (size(C,1) == size(D,1) == size(V,1)) ?
        new(C,D,V) : error("first dimensions not matching")
    =#
end

function LinearObservationModel(C::AbstractMatrix, D::AbstractMatrix,
                                V::AbstractMatrix)
    return LinearObservationModel(C,D,Symmetric(V))
end

function LinearObservationModel(C::AbstractMatrix, V::Symmetric)
    n = size(C,1)
    return LinearObservationModel(C,zeros(Bool,n,n),V)
end

function LinearObservationModel(C::AbstractMatrix, V::AbstractMatrix)
    n = size(C,1)
    return LinearObservationModel(C,zeros(Bool,n,n),Symmetric(V))
end

"""
    measure(m::LinearObservationModel, x::AbstractVector{<:Number}, u::AbstractVector{<:Number})
    measure(m::LinearObservationModel, x::AbstractVector{T}, u::AbstractVector{T}, rng::AbstractRNG) where T<:Number

Returns an observation of state x according to the linear observation model m, with control inputs u.
If rng is passed, adds additive Gaussian noise to the observation.
"""
function measure(m::LinearObservationModel, x::AbstractVector{<:Number},
                 u::AbstractVector{<:Number})
    return m.C * x + m.D * u
end
function measure(m::LinearObservationModel, x::AbstractVector{T},
                 u::AbstractVector{T},
                 rng::AbstractRNG) where T<:Number
    return m.C * x + m.D * u + cholesky(m.V).L * randn(rng, T, size(m.V, 1))
end

function jacobian(m::LinearObservationModel, x, u)
    m.C
end

"""
    NonlinearDynamicsModel(f::Function,W::Symmetric)
    NonlinearDynamicsModel(f::Function,W::AbstractMatrix)

Construct nonlinear dynamics model with transition function f
and symmetric zero-mean process noise with symmetric covariance matrix W
"""
struct NonlinearDynamicsModel{F<:Function, c<:Symmetric} <: DynamicsModel
    f::F
    W::c
end

function NonlinearDynamicsModel(f::Function, W::AbstractMatrix)
    return NonlinearDynamicsModel(f,Symmetric(W))
end

"""
    predict(m::NonLinearDynamicsModel, x::AbstractVector{<:Number}, u::AbstractVector{<:Number})
    predict(m::NonLinearDynamicsModel, x::AbstractVector{<:Number}, u::AbstractVector{<:Number}, rng::AbstractRNG)

Uses the non linear dynamics model to propagate the state x one step forward in time with control input u.
If rng is given, it adds process noise.
"""
function predict(m::NonlinearDynamicsModel, x::AbstractVector{<:Number},
                 u::AbstractVector{<:Number})
    return m.f(x, u)
end
function predict(m::NonlinearDynamicsModel, x::AbstractVector{T},
                 u::AbstractVector{T},
                 rng::AbstractRNG) where T<:Number
    return m.f(x, u) + cholesky(m.W).L * randn(rng, T, size(m.W, 1))
end

function jacobian(m::NonlinearDynamicsModel, x::AbstractVector, u::AbstractVector)
    return ForwardDiff.jacobian(μ -> m.f(μ, u), x)
end

"""
    NonlinearObservationModel(h::Function,V::Symmetric)
    NonlinearObservationModel(h::Function,V::AbstractMatrix)

Construct nonlinear observation dynamics model with measurement function h
and symmetric zero-mean measurement noise with symmetric covariance matrix V
"""
struct NonlinearObservationModel{F<:Function, c<:Symmetric} <: ObservationModel
    h::F
    V::c
end

function NonlinearObservationModel(h::Function, V::AbstractMatrix)
    return NonlinearObservationModel(h,Symmetric(V))
end

"""
    measure(m::NonlinearObservationModel, x::AbstractVector{<:Number}, u::AbstractVector{<:Number})
    measure(m::NonlinearObservationModel, x::AbstractVector{T}, u::AbstractVector{T}, rng::AbstractRNG) where T<:Number

Returns an observation of state x according to the non linear observation model m, with control inputs u.
If rng is passed, adds additive Gaussian noise to the observation.
"""
function measure(m::NonlinearObservationModel, x::AbstractVector{<:Number},
                 u::AbstractVector{<:Number})
    return m.h(x, u)
end
function measure(m::NonlinearObservationModel, x::AbstractVector{T},
                 u::AbstractVector{T},
                 rng::AbstractRNG) where T<:Number
    return m.h(x, u) + cholesky(m.V).L * randn(rng, T, size(m.V, 1))
end

function jacobian(m::NonlinearObservationModel, x::AbstractVector, u::AbstractVector)
    return ForwardDiff.jacobian(μ -> m.h(μ, u), x)
end


text:
module TestPseudoInverse
using LinearAlgebra
using LazyInverses
using Test

function scalar_tests(A, M, Inv)
    elty = eltype(A)
    x = randn(elty)
    @testset "scalar tests with eltype $elty" begin
        @test Matrix(Inv * x) ≈ M * x
        @test Matrix(x * Inv) ≈ x * M
        @test Matrix(Inv / x) ≈ M / x
        @test Matrix(x \ Inv) ≈ x \ M
        # @test Matrix(Inv \ x) ≈ A * x
        # @test Matrix(x / Inv) ≈ x * A
    end
end

@testset "pseudoinverse" begin
    n, m = 3, 2
    element_types = (Float32, Float64, ComplexF32, ComplexF64)
    for elty in element_types
        A = randn(elty, n, m)
        a = randn(elty, n)
        b = randn(elty, m)
        B = randn(elty, m, n)

        # pseudoinverse
        LInv = pseudoinverse(A)
        @test LInv isa PseudoInverse{<:elty}
        ML = Matrix(LInv)
        @test ML ≈ pinv(A)

        @test Matrix(adjoint(LInv)) ≈ adjoint(ML)
        @test Matrix(transpose(LInv)) ≈ transpose(ML)

        # scalar operations
        scalar_tests(A, ML, LInv)

        # vector operations
        @test LInv * a ≈ ML * a
        @test b' * LInv ≈ b' * ML

        # matrix operations
        @test LInv * A ≈ I(m)
        @test ML * A ≈ I(m)
        @test A * LInv ≈ A * ML
        @test B' * LInv ≈ B' * ML
        @test LInv * B' ≈ ML * B'
        @test B * LInv' ≈ B * ML'

        # factorize
        @test factorize(LInv) ≡ LInv

        # smart constructor
        A = randn(elty, n, n)
        @test pinverse(A) isa Inverse

        x = rand(elty, (1, 1))
        @test pinverse(x) isa elty
        @test pinverse(x) == inv(x[1]) # smart pseudo-constructor returns scalar on 1x1 matrix
        @test pinverse(Inverse(x)) == x[1]

        # scalar
        x = rand(elty)
        PInv = PseudoInverse(x)
        @test PInv * x ≈ 1
        @test x * PInv ≈ 1
        @test pinverse(x) ≈ 1/x
        @test AbstractMatrix(PInv) ≈ fill(1/x, (1, 1))
    end
end

end # module


text:
using Printf

struct NCMresults
    X::Symmetric{Float64,Array{Float64,2}}
    y::Vector{Float64}
    Λ::Symmetric{Float64,Array{Float64,2}}
    fgcountRef::Base.RefValue{Int32}
    fvals::Vector{Float64}
    resvals::Vector{Float64}
    distvals::Vector{Float64}
    rpRef::Base.RefValue{Float64}
    rdRef::Base.RefValue{Float64}
    εRef::Base.RefValue{Float64}

    function NCMresults(n, maxfgcalls)
        X = Symmetric(zeros(n, n))
        y = zeros(n)
        Λ = Symmetric(zeros(n, n))
        fgcountRef = Ref{Int32}(0)
        fvals = Vector{Float64}(undef, maxfgcalls)
        resvals = Vector{Float64}(undef, maxfgcalls)
        distvals = Vector{Float64}(undef, maxfgcalls)
        rpRef = Ref{Float64}(0.0)
        rdRef = Ref{Float64}(0.0)
        εRef = Ref{Float64}(0.0)

        new(X, y, Λ, fgcountRef, fvals, resvals, distvals, rpRef, rdRef, εRef)
    end
end

struct NCM
    n::Int32
    memlim::Int32
    maxfgcalls::Int32
    g::Vector{Float64}
    d::Vector{Float64}
    M::Symmetric{Float64,Array{Float64,2}}
    R::Symmetric{Float64,Array{Float64,2}}
    H2::Symmetric{Float64,Array{Float64,2}}
    Y::Symmetric{Float64,Array{Float64,2}}
    ∇fY::Symmetric{Float64,Array{Float64,2}}
    Γ::Symmetric{Float64,Array{Float64,2}}
    V::Symmetric{Float64,Array{Float64,2}}
    X::Symmetric{Float64,Array{Float64,2}}
    Z::Symmetric{Float64,Array{Float64,2}}
    Rd::Symmetric{Float64,Array{Float64,2}}
    Xold::Symmetric{Float64,Array{Float64,2}}
    wa::Vector{Float64}
    iwa::Vector{Int32}
    nbd::Vector{Int32}
    lower::Vector{Float64}
    upper::Vector{Float64}
    task::Vector{UInt8}
    task2::SubArray{UInt8,1,Vector{UInt8},Tuple{UnitRange{Int64}},true}
    csave::Vector{UInt8}
    lsave::Vector{Int32}
    isave::Vector{Int32}
    dsave::Vector{Float64}
    nRef::Base.RefValue{Int32}
    mRef::Base.RefValue{Int32}
    iprint::Base.RefValue{Int32}
    fRef::Base.RefValue{Float64}
    factr::Base.RefValue{Float64}
    pgtol::Base.RefValue{Float64}
    proj::ProjPSD
    res::NCMresults

    function NCM(n; memlim = 10, maxfgcalls = 100_000)
        g = zeros(n)
        d = zeros(n)

        M = Symmetric(zeros(n, n))
        R = copy(M)
        H2 = copy(M)
        Y = copy(M)
        ∇fY = copy(M)
        Γ = copy(M)
        V = copy(M)
        X = copy(M)
        Z = copy(M)
        Rd = copy(M)
        Xold = copy(M)

        nmax = n
        mmax = memlim

        wa = zeros(Cdouble, 2mmax * nmax + 5nmax + 11mmax * mmax + 8mmax)
        iwa = zeros(Cint, 3nmax)

        # provide nbd which defines the bounds on the variables:
        nbd = zeros(Cint, nmax)       # no bounds on the variables
        lower = zeros(Cdouble, nmax)    # the lower bounds
        upper = zeros(Cdouble, nmax)    # the upper bounds

        task = fill(Cuchar(' '), 60)   # fortran's blank padding
        task2 = view(task, 1:2)
        csave = fill(Cuchar(' '), 60)   # fortran's blank padding
        lsave = zeros(Cint, 4)
        isave = zeros(Cint, 44)
        dsave = zeros(Cdouble, 29)

        nRef = Ref{Cint}(n)
        mRef = Ref{Cint}(memlim)
        iprint = Ref{Cint}(-1)
        fRef = Ref{Cdouble}(0.0)
        factr = Ref{Cdouble}(0.0)
        pgtol = Ref{Cdouble}(0.0)

        proj = ProjPSD(n)

        res = NCMresults(n, maxfgcalls)

        new(
            n,
            memlim,
            maxfgcalls,
            g,
            d,
            M,
            R,
            H2,
            Y,
            ∇fY,
            Γ,
            V,
            X,
            Z,
            Rd,
            Xold,
            wa,
            iwa,
            nbd,
            lower,
            upper,
            task,
            task2,
            csave,
            lsave,
            isave,
            dsave,
            nRef,
            mRef,
            iprint,
            fRef,
            factr,
            pgtol,
            proj,
            res,
        )
    end
end

function (ncm::NCM)(
    G::Symmetric{Float64,Array{Float64,2}},
    H::Symmetric{Float64,Array{Float64,2}};
    method::Symbol = :IAPG,
    τ::Float64 = 1.0,
    α::Float64 = 0.0,
    σ::Float64 = 1.0,
    tol::Float64 = 1e-1,
    kmax::Int64 = 100_000,
    maxfgcalls::Int64 = 100_000,
    printlevel::Int64 = 1,
    lbfgsbprintlevel::Int64 = -1,
    exact::Bool = false,
    useXold::Bool = false,
    innerverbose::Bool = false,
    cleanvals::Bool = true,
    scaleX::Bool = true,
)

    # Loss function and gradient
    #f(X) = 0.5*norm(H.*(X .- G))^2
    #∇f(X) = Symmetric(H2.*(X .- G))

    n = size(G, 1)
    validmethod = (method == :IAPG || method == :IR || method == :IER)

    # Check for valid input
    if maxfgcalls > ncm.maxfgcalls
        error("require maxfgcalls ≤ ncm.maxfgcalls")
    end
    validmethod || error("method must be :IAPG or :IR or :IER")
    n == ncm.n || error("require n == ncm.n")
    size(G) == size(H) || error("G and H must be the same size")
    issymmetric(G) || error("G must be symmetric")
    issymmetric(H) || error("H must be symmetric")
    !iszero(H) || error("H must be nonzero")

    g = ncm.g
    H2 = ncm.H2
    Y = ncm.Y
    V = ncm.V
    R = ncm.R
    Xold = ncm.Xold

    proj = ncm.proj
    res = ncm.res

    Xnew = res.X
    y = res.y
    Λ = res.Λ

    Xnew .= Xold

    fgcountRef = res.fgcountRef
    rpRef = res.rpRef
    rdRef = res.rdRef
    εRef  = res.εRef
    fvals = res.fvals

    H2.data .= H .^ 2

    # Lipschitz constant of ∇f
    L = fronorm(H2, proj.work)

    if method == :IAPG
        τ == 1 || error("IAPG method requires τ = 1")
        t0 = 1.0
    end

    if method == :IR
        0 < τ ≤ 1 || error("IR method requires 0 < τ ≤ 1")
        0 ≤ α ≤ (1 - τ) * (L / τ) ||
        error("IR method requires 0 ≤ α ≤ $((1 - τ)*(L/τ))")
        t0 = 1.0
    end

    if method == :IER
        τ == 1 || error("IER method requires τ = 1")
        # α = 19/L
        α > 0 || error("IER method requires α > 0")
        0 ≤ σ ≤ 1 || error("IER method requires 0 ≤ σ ≤ 1")
        λ = α / (1 + α * L)
        t0 = 0.0
    end

    printlevel ≥ 1 && println("$method method, τ=$τ, α=$α, σ=$σ, tol=$tol")

    if !useXold
        fill!(Xold, 0.0)
    end

    fill!(y, 0.0)

    k = 0
    t = t0
    Y .= Xold
    innertol = NaN
    rp = rd = Inf
    innersuccess = true
    fgcountRef[] = 0
    fgcount = fgcountRef[]

    while (innersuccess &&
        max(rp, rd) > tol && k < kmax && fgcount < maxfgcalls
    )

        k += 1

        if method == :IER || method == :ER
            tnew = t + (λ + √(λ^2 + 4λ * t)) / 2
            Y.data .=
                (t / tnew) .* Xnew.data .+ ((tnew - t) / tnew) .* Xold.data
        else
            tnew = (1 + √(1 + 4 * t^2)) / 2
        end

        if method == :IAPG || method == :ER
            innertol = 1 / t^2
        end

        if exact
            innertol = 0.0
        end

        maxinnerfgcalls = maxfgcalls - fgcount

        # Solve the subproblem
        innersuccess = calllbfgsb!(
            ncm,
            G,
            H,
            tol,
            t,
            L,
            τ,
            α,
            σ;
            method = method,
            maxfgcalls = maxinnerfgcalls,
            innertol = innertol,
            exact = exact,
            verbose = innerverbose,
            lbfgsbprintlevel = lbfgsbprintlevel,
            cleanvals = cleanvals,
            scaleX = scaleX,
        )
        if !innersuccess
            printlevel ≥ 2 && println("Failed to solve subproblem.")
        end
        fgcount = fgcountRef[]
        innerfgcalls = fgcount - (maxfgcalls - maxinnerfgcalls)

        rp, rd, ε = rpRef[], rdRef[], εRef[]

        if printlevel ≥ 2
            mod(k, 20) == 1 && @printf(
                "%4s %8s %10s %10s %10s %10s %10s %10s\n",
                "k",
                "fgcalls",
                "||g||",
                "innertol",
                "f(X)",
                "rp",
                "rd",
                "ε"
            )
            @printf(
                "%4d %8d %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e\n",
                k,
                innerfgcalls,
                norm(g),
                innertol,
                fvals[fgcount],
                rp,
                rd,
                ε
            )
        end

        # Update
        if method == :IAPG
            Y.data .= Xnew.data .+ ((t - 1) / tnew) .* (Xnew.data .- Xold.data)
            Xold .= Xnew
        elseif method == :IR
            Y.data .=
                Xnew.data .- ((t / tnew) * (τ / L)) .* V.data .+
                ((t - 1) / tnew) .* (Xnew.data .- Xold.data)
            Xold .= Xnew
        elseif method == :IER
            Xold.data .-= (tnew - t) .* (V.data .+ (Y.data .- Xnew.data)./λ)
        end

        t = tnew
    end

    success = (max(rp, rd) ≤ tol)
    if printlevel ≥ 1
        resnorm = fronorm(R, proj.work)
        #pval = fvals[fgcount]
        #dval = 0.5*(dot(G, H2.*G) - dot(Xnew, H2.*Xnew)) + sum(y)
        println(success ? "Success." : "Failed.")
        @printf("%-20s: %-24d\n", "Outer iterations", k)
        @printf("%-20s: %-24d\n", "Function evals", fgcount)
        @printf("%-20s: %-24.16e\n", "Final ||H∘(X-G)||", resnorm)
        #@printf("%-20s: %-24.16e\n", "Primal objective",   pval)
        #@printf("%-20s: %-24.16e\n", "Dual objective",     dval)
        @printf("%-20s: %-24.6e\n", "Primal feasibility", rp)
        @printf("%-20s: %-24.6e\n", "Dual feasibility", rd)
        @printf("%-20s: %-24.6e\n", "Complementarity", symdot(Xnew, Λ))
        #@printf("%-20s: %-24d\n",    "rank(X)",            proj.m[])
    end

    return success, k
end

max_stars_repo_name: krislock/IR-FISTA

max_issues_repo_name: krislock/IR-FISTA

max_forks_repo_name: krislock/IR-FISTA


text:
using OrdinaryDiffEq, Test

f_ec(u,p,t) = exp(u)
u0 = 0.0 # explosion time is 1.0
tspan = (0.0, 10.0)
prob = ODEProblem(f_ec,u0,tspan)
options = [:reltol => 1e-8, :abstol => 1e-8, :verbose => false]
desired_code = :DtLessThanMin

# Test that sol.retcode is set to the correct value by various ways to
# invoke integrator.

sol = solve(prob,Tsit5(); options...)
@test sol.retcode == desired_code

integrator = init(prob,Tsit5(); options...)
solve!(integrator)
@test integrator.sol.retcode == desired_code

integrator = init(prob,Tsit5(); options...)
for _ in integrator end
@test integrator.sol.retcode == desired_code

integrator = init(prob,Tsit5(); options...)
step!(integrator, 10.0)
@test integrator.sol.retcode == desired_code

# Test check_error
integrator = init(prob,Tsit5(); options...)
step!(integrator)
@test check_error(integrator) == :Success
ok = false
for i in 1:integrator.opts.maxiters
  step!(integrator)
  if check_error(integrator) == desired_code
    global ok = true
    # @show i
    break
  end
end
@test ok


text:
struct MetricsScore
    name::String
    r::Array{Float64}
    p::Array{Float64}
    fpr::Array{Float64}
    tpr::Array{Float64}
    aupr::Array{Float64}
    auroc::Array{Float64}
end

show(io::IO,x::MetricsScore)=print(io,"The metrics of $(x.name)")

function compared_scores(data_dir::String)
    names=[
        "CBPred",
        "LRSSL ",
        "SCMFDD",
        "HGBI  ",
        "MBiRW "
    ]
    scores=MetricsScore[]
    for (i,name) in enumerate(names)
        P=readdlm(joinpath(data_dir,"$(i)mean5_P.txt"))[:]
        FPR=readdlm(joinpath(data_dir,"$(i)mean5_FPR.txt"))[:]
        TPR=readdlm(joinpath(data_dir,"$(i)mean5_TPR.txt"))[:]
        R=TPR
        AUROC=readdlm(joinpath(data_dir,"$(i)auc.txt"))[:]
        AUPR=readdlm(joinpath(data_dir,"$(i)aupr.txt"))[:]
        push!(scores,MetricsScore(name,R,P,FPR,TPR,AUPR,AUROC))
    end
    scores
end

function pr_roc_curve(scores::Array{MetricsScore})
    curve_roc=plot(title="ROC",xlabel="Recall",ylabel="Precision")
    curve_pr=plot(title="PR",xlabel="FPR",ylabel="TPR")
    for s in scores
        n=s.name
        au=auc(s.r,s.p)
        plot!(curve_pr,s.r,s.p,label="$n $(round(au,digits=3))")
        au=auc(s.fpr,s.tpr)
        plot!(curve_roc,s.fpr,s.tpr,label="$n $(round(au,digits=3))")
    end
    plot(curve_pr,curve_roc,layout=(1,2))
end

function topk_recall_curve(scores::Array{MetricsScore},ks::StepRange)
    666
end

export compared_scores,MetricsScore,pr_roc_curve,topk_recall_curve

max_stars_repo_name: yeruoforever/DrugRepositioning

max_issues_repo_name: yeruoforever/DrugRepositioning

max_forks_repo_name: yeruoforever/DrugRepositioning


text:
module COBS

export encode, decode

function _doblock(xs)
    n = length(xs)
    ys = Vector{UInt8}(undef, n + 1)
    ind = 0x01
    acc = 0x01
    for x in reverse(xs)
        if iszero(x)
            ys[ind] = acc
            acc = 0x00
        else
            ys[ind] = x
        end
        ind += 0x01
        acc += 0x01
    end
    ys[end] = acc
    return reverse(ys)
end

"""
    encode(sp, payload)

Encode each 254-byte block in the `payload` and directly send
the data through a specified serial port, `sp`. This is faster
than first encoding the whole payload and then seding it.
"""
function encode(sp, payload)
    for xs in Base.Iterators.partition(payload, 254)
        block = _doblock(xs)
        write(sp, block)
    end
    write(sp, 0x00)
end

"""
    encode(payload)

Return the encoded `payload`.
"""
function encode(payload)
    blocks = Vector{Vector{UInt8}}()
    for xs in Base.Iterators.partition(payload, 254)
        block = _doblock(xs)
        push!(blocks, block)
    end
    vcat(blocks..., 0x00)
end

#=function _f!(y, msg, i1, i2)
    for i in i1 + 1:i2 - 1
        y[i - 1] = msg[i]
    end
    return i2, i2 + msg[i2]
end=#

#=function decode(msg)
    y = Vector{UInt8}(undef, length(msg) - 2)
    i1 = 0x01
    i2 = i1 + msg[i1]
    i1, i2 = _f!(y, msg, i1, i2)
    while i1 < i2
        y[i1 - 1] = 0
        i1, i2 = _f!(y, msg, i1, i2)
    end
    return y
end=#

"""
    decode(msg)

Return the decoded `msg`.
"""
function decode(_msg::AbstractVector)
    msg = copy(_msg)
    pl = UInt8[]
    n = popfirst!(msg)
    c = 0
    b = popfirst!(msg)
    while b ≠ 0
        c += 1
        if c < n
            push!(pl, b)
        else
            push!(pl, 0)
            n = b
            c = 0
        end
        b = popfirst!(msg)
    end
    return pl
end

"""
    decode(sp)

Decode a message directly from a specified serial port, `sp`.
"""
function decode(sp)
    pl = UInt8[]
    n = read(sp, UInt8)
    c = 0
    b = read(sp, UInt8)
    while b ≠ 0
        c += 1
        if c < n
            push!(pl, b)
        else
            push!(pl, 0)
            n = b
            c = 0
        end
        b = read(sp, UInt8)
    end
    return pl
end

end # module


text:
############################################################
# matscipy neighbourlist functionality
############################################################

include("MatSciPy.jl")

function neighbourlist(at::ASEAtoms, cutoff::Float64)::MatSciPy.NeighbourList
   # TODO: also recompute if rcut is different !!!!!
   # if no previous neighbourlist is available, compute a new one
   # if !has_transient(at, (:nlist, cutoff)) || recompute
   #    # this nlist will be destroyed as soon as positions change
   #    set_transient!(at, (:nlist, cutoff), MatSciPy.NeighbourList(at, cutoff))
   # end
   # return get_transient(at, (:nlist, cutoff))
   return MatSciPy.NeighbourList(at, cutoff)
end

# TODO: no longer needed in the new ASE design
# """
# `static_neighbourlist(at::ASEAtoms, rcut::Float64)`
#
# This function first checks whether a static neighbourlist already exists
# with cutoff `rcut` and if it does then it returns the existing list.
# If it does not, then it computes a new neighbour list with the current
# configuration, stores it for later use and returns it.
# """
# function static_neighbourlist(at::ASEAtoms, rcut::Float64)
#    # if !has_transient(at, (:snlist, rcut))
#    #    set_transient!( at, (:snlist, rcut),
#    #                        MatSciPy.NeighbourList(at, rcut),
#    #                        Inf )    # Inf means this is never deleted!
#    # end
#    # return get_transient(at, (:snlist, rcut))
#    return MatSciPy.NeighbourList(at, rcut)
# end


text:
include("GridGraphs.jl")
@reexport using Gridap.GridGraphs

include("Grids.jl")
@reexport using Gridap.Grids

include("FaceLabels.jl")
@reexport using Gridap.FaceLabelsModule

include("DiscreteModels.jl")
@reexport using Gridap.DiscreteModels

include("UnstructuredGeometry.jl")
@reexport using Gridap.UnstructuredGeometry

include("GeometryWrappers.jl")
@reexport using Gridap.GeometryWrappers

include("CartesianGeometry.jl")
@reexport using Gridap.CartesianGeometry

include("GridPortions.jl")
@reexport using Gridap.GridPortions

include("BoundaryGrids.jl")
@reexport using Gridap.BoundaryGrids

include("SkeletonGrids.jl")
@reexport using Gridap.SkeletonGrids

include("Simplexify.jl")
@reexport using Gridap.Simplexify


text:
"""
x = fnnls(AtA, Atb; ...)

Returns x that solves A*x = b in the least-squares sense, subject to x >=0. The
inputs are the cross-products AtA = A'*A and Atb = A'*b. Uses the modified
active set method of Bro and De Jong (1997).

Optional arguments:
    tol: tolerance for nonnegativity constraints
    max_iter: maximum number of iterations (counts inner loop iterations)

References:
    Bro R, De Jong S. A fast non-negativitity-constrained least squares
    algorithm. Journal of Chemometrics. 11, 393–401 (1997)
"""
function fnnls(AtA, Atb::AbstractVector{T}; kwargs...) where {T}
    work = fnnls_work(AtA, Atb)
    fnnls!(work, AtA, Atb; kwargs...)
end

function fnnls_work(AtA, Atb::AbstractVector{T}) where {T}
    n    = size(AtA,1)
    x    = zeros(T, n)
    s    = zeros(T, n)
    w    = similar(x)
    P    = BitArray(undef, n)
    nP   = BitArray(undef, n)
    iP   = BitArray(undef, n)
    AtAP = similar(AtA)
    AtbP = similar(Atb)
    return (x=x, s=s, w=w, P=P, nP=nP, iP=iP, AtAP=AtAP, AtbP=AtbP)
end

function fnnls!(
        work,
        AtA,
        Atb::AbstractVector{T};
        tol::T = sqrt(eps(T)),
        max_iter::Int = 30 * size(AtA, 2)
    ) where {T}

    x, s, w, P, nP, iP, AtAP, AtbP = work
    n = size(AtA,1)
    x .= 0
    s .= 0

    # P is a bool array storing positive elements of x
    # i.e., x[P] > 0 and x[~P] == 0
    P  .= x .> tol
    nP .= .!(P)

    mul!(w, AtA, x)
    w .= Atb .- w; #w = Atb - AtA*x

    function anyP(f, x, P)
        @inbounds for i in eachindex(x,P)
            if P[i] && f(x[i])
                return true
            end
        end
        return false
    end

    function argmaxP(x, P)
        α = T(-Inf)
        out = 1
        @inbounds for i in eachindex(x, P)
            if P[i] && (x[i] > α)
                α = max(α, x[i])
                out = i # i = argmax(w .* nP)
            end
        end
        return out
    end

    function assigntoP!(dst, srcP, P)
        j = 0
        @inbounds for i in eachindex(dst, P)
            if P[i]
                j += 1
                dst[i] = srcP[j]
            end
        end
        return dst
    end

    function copytocorner!(dst, src, P)
        j1, j2 = 0, 0
        @inbounds for i2 in 1:size(src, 2)
            if P[i2]
                j2 += 1
                j1  = 0
                for i1 in 1:size(src, 1)
                    if P[i1]
                        j1 += 1
                        dst[j1,j2] = src[i1,i2]
                    end
                end
            end
        end
        return j1, j2
    end

    function copytostart!(dst, src, P)
        j1, j2, idx = 0, 0, 0
        @inbounds for i2 in 1:size(src, 2)
            if P[i2]
                j2 += 1
                j1  = 0
                for i1 in 1:size(src, 1)
                    if P[i1]
                        j1  += 1
                        idx += 1
                        dst[idx] = src[i1,i2]
                    end
                end
            end
        end
        return j1, j2
    end

    function copytotop!(dst, src, P)
        j1 = 0
        @inbounds for i1 in 1:size(src, 1)
            if P[i1]
                j1 += 1
                dst[j1] = src[i1]
            end
        end
        return j1
    end

    function alpha(x, s, P)
        α = oftype(Inf, one(eltype(x))/one(eltype(s)))
        @inbounds for i in eachindex(x, s, P)
            if P[i]
                α = min(α, x[i]/(x[i] - s[i]))
            end
        end
        return α
    end

    function solveP!(y, Atmp, ytmp, A, x, P) # s[P] = AtA[P,P] \ Atb[P]
        # (1)
        # y[P] .= A[P,P] \ x[P] # 644 - 721 us
        # (2)
        # assigntoP!(y, Symmetric(A[P,P]) \ x[P], P) # 723-836 us
        assigntoP!(y, A[P,P] \ x[P], P) # 611-709 us
        # (3)
        # invA = cholesky!(Symmetric(A[P,P]); check = false) # 650 us
        # invA = cholesky!(A[P,P]; check = false) # 650-850 ua
        # assigntoP!(y, ldiv!(invA, x[P]), P)
        # (4)
        # m_, n_ = copytocorner!(Atmp, A, P)
        # p_ = copytotop!(ytmp, x, P)
        # invA = cholesky!(Symmetric(view(Atmp, 1:m_, 1:n_)); check = false)
        # outP = ldiv!(invA, view(ytmp, 1:p_)) # 630-720 us
        # assigntoP!(y, outP, P)
        # (5)
        # m_, n_ = copytostart!(Atmp, A, P)
        # p_ = copytotop!(ytmp, x, P)
        # Atmp_view = reshape(view(Atmp, 1:m_*n_), m_, n_)
        # # invA = cholesky!(Atmp_view; check = false) # 660 - 730 us
        # invA = cholesky!(Symmetric(Atmp_view); check = false) # 650 - 715 us
        # outP = ldiv!(invA, view(ytmp, 1:p_))
        # assigntoP!(y, outP, P)
    end

    # We have reached an optimum when either:
    #   (a) all elements of x are positive (no nonneg constraints activated)
    #   (b) ∂f/∂x = A' * (b - A*x) > 0 for all nonpositive elements of x
    iter = 0
    while sum(P) < n && anyP(>(tol), w, nP) && iter < max_iter

        # find i that maximizes w, restricting i to indices not in P
        # Note: the while loop condition guarantees at least one w[~P]>0
        i = argmaxP(w, nP) # i = argmax(w .* nP)

        # Move i to P
        P[i]  = true
        nP[i] = false

        # Solve least-squares problem, with zeros for columns/elements not in P
        solveP!(s, AtAP, AtbP, AtA, Atb, P) # s[P] = AtA[P,P] \ Atb[P]
        @inbounds for i in eachindex(s, nP)
            nP[i] && (s[i] = 0) # s[nP] .= zero(eltype(s)) # zero out elements not in P
        end

        # Inner loop: deal with negative elements of s
        while anyP(<=(tol), s, P)
            iter += 1

            # find indices in P where s is negative
            @inbounds iP .= (s .<= tol) .& P

            # calculate step size, α, to prevent any xᵢ from going negative
            α = alpha(x, s, iP) # minimum(x[iP] ./ (x[iP] .- s[iP]))

            # update solution (pushes some xᵢ to zero)
            @inbounds x .+= α .* (s .- x)

            # Remove all i in P where x[i] == 0
            @inbounds for i = 1:n
                if P[i] && abs(x[i]) < tol
                    P[i]  = false # remove i from P
                    nP[i] = true
                end
            end

            # Solve least-squares problem again, zeroing nonpositive columns
            solveP!(s, AtAP, AtbP, AtA, Atb, P) # s[P] = AtA[P,P] \ Atb[P]
            @inbounds for i in eachindex(s, nP)
                nP[i] && (s[i] = 0) # s[nP] .= zero(eltype(s)) # zero out elements not in P
            end
        end

        # update solution
        x .= s
        mul!(w, AtA, x)
        w .= Atb .- w; #w = Atb - AtA*x
    end

    return x
end

function fnnls(A,
               B::AbstractMatrix;
               gram::Bool = false,
               use_parallel::Bool = true,
               kwargs...)

    n = size(A,2)
    k = size(B,2)

    if gram
        # A,B are actually Gram matrices
        AtA = A
        AtB = B
    else
        # cache matrix computations
        AtA = A'*A
        AtB = A'*B
    end

    if use_parallel && nprocs()>1
        X = @distributed (hcat) for i=1:k
            fnnls(AtA, AtB[:,i]; kwargs...)
        end
    else
        X = Array{eltype(B)}(undef,n,k)
        for i = 1:k
            X[:,i] = fnnls(AtA, AtB[:,i]; kwargs...)
        end
    end

    return X
end


text:
function  factorial_digit_sum(n::Int)
	ans=0

	while( n>=1)
		ans+=factorial( n%10);
		n=floor(Int,n/10)
	end

	return ans

end

function compute()

	ans =sum( i for i =3:10000000 if i==factorial_digit_sum(i))

	return ans
end

println(compute())


text:
abstract type CovarianceDesign{T} end

get_Σ(mat) = mat
nfeatures(mat) = size(mat, 1)

function spectrum(mat)
    eigs = eigvals(mat)
    probs = fill(1/length(eigs), length(eigs))
    DiscreteNonParametric(eigs, probs)
end

spectrum(cov::CovarianceDesign) = spectrum(get_Σ(cov))

nfeatures(cov::CovarianceDesign) = cov.p

function simulate_rotated_design(cov, n; rotated_measure = Normal())
    Σ = get_Σ(cov)
    Σ_chol = cholesky(Σ)
    p = nfeatures(cov)
    Z = rand(rotated_measure, n, p)
    X = randn(n, p) * Σ_chol.UL
    X
end

Base.@kwdef struct AR1Design{P<:Union{Missing,Int}} <: CovarianceDesign{P}
    p::P = missing
    ρ = 0.7
end

function get_Σ(cov::AR1Design{Int})
    p = nfeatures(cov)
    ρ = cov.ρ
    Σ = [ρ^(abs(i - j)) for i = 1:p, j = 1:p]
    Σ
end

abstract type DiagonalCovarianceDesign{T} <: CovarianceDesign{T} end

Base.@kwdef struct IdentityCovarianceDesign{P<:Union{Missing,Int}} <:
                   DiagonalCovarianceDesign{P}
    p::P = missing
end

spectrum(::IdentityCovarianceDesign) = DiscreteNonParametric([1.0],[1.0])

function get_Σ(cov::IdentityCovarianceDesign{Int})
    I(cov.p)
end

Base.@kwdef struct UniformScalingCovarianceDesign{P<:Union{Missing,Int}} <:
                   DiagonalCovarianceDesign{P}
    scaling::Float64 = 1.0
    p::P = missing
end

spectrum(unif::UniformScalingCovarianceDesign) = DiscreteNonParametric([unif.scaling],[1.0])

function get_Σ(cov::UniformScalingCovarianceDesign{Int})
    (cov.scaling * I)(cov.p)
end

Base.@kwdef struct ExponentialOrderStatsCovarianceDesign{P<:Union{Missing,Int}} <:
                   DiagonalCovarianceDesign{P}
    p::P = missing
    rate::Float64
end

function spectrum(cov::ExponentialOrderStatsCovarianceDesign)
    p = cov.p
    rate = cov.rate
    tmp = range(1 / (2p); stop = 1 - 1 / (2p), length = p)
    eigs = 1 / rate .* log.(1 ./ tmp)
    DiscreteNonParametric(eigs, fill(1/p, p))
end

get_Σ(cov::ExponentialOrderStatsCovarianceDesign) = Diagonal(support(spectrum(cov)))

struct BlockCovarianceDesign{T, S <: CovarianceDesign{T}, G} <: CovarianceDesign{T}
    blocks::Vector{S}
    groups::G
end

function BlockCovarianceDesign(blocks::Vector{S}) where S<:CovarianceDesign{Missing}
    BlockCovarianceDesign(blocks, missing)
end

nfeatures(cov::BlockCovarianceDesign) = sum(nfeatures.(cov.blocks))

function get_Σ(blockdesign::BlockCovarianceDesign)
    BlockDiagonal(get_Σ.(blockdesign.blocks))
end

function spectrum(blockdesign::BlockCovarianceDesign)
    @unpack blocks, groups = blockdesign
    spectra = spectrum.(blocks)
    mixing_prop = groups.ps ./ groups.p
    MixtureModel(spectra, mixing_prop)
end

function simulate_rotated_design(cov::BlockCovarianceDesign, n; rotated_measure = Normal())
    hcat(simulate_rotated_design.(cov.blocks, n; rotated_measure=rotated_measure)...)
end

# Set groups
function set_groups(design::CovarianceDesign, p::Integer)
    @set design.p = p
end

function set_groups(design::CovarianceDesign, groups::GroupedFeatures)
    set_groups(design, nfeatures(groups))
end

function set_groups(blockdesign::BlockCovarianceDesign, groups::GroupedFeatures)
    updated_blocks = set_groups.(blockdesign.blocks, groups.ps)
    BlockCovarianceDesign(updated_blocks, groups)
end