text:
"""
    funntion simpleFirstLoop_(platform, code)

Returns the sent code wrapped in parallelized loop over the particles of the system adapted for each platform.
"""
function simpleFirstLoop_(platform::String, code::Expr)

    if platform == "cpu"
        code =
            quote
                Threads.@threads for ic1_ in 1:N
                    $code
                end
            end
    elseif platform == "gpu"
        code =
            quote
                index_ = (threadIdx().x) + (blockIdx().x - 1) * blockDim().x
                stride_ = blockDim().x * gridDim().x
                for ic1_ in index_:stride_:N
                    $code
                end
            end
    else
        error("Platform should be or cpu or gpu. ", platform, " was given.")
    end

    return code
end

"""
    funntion simpleGridLoop_(platform, code)

Returns the sent code wrapped in parallelized loop over the grid dimensions of the system adapted for each platform.
"""
function simpleGridLoop_(platform::String, code::Expr, nloops::Int; indexes = [1,2,3])

    #Compute indexes
    vIc = [:ic1_,:ic2_,:ic3_][indexes[1:nloops]]
    vIndex = [:indexX_,:indexY_,:indexZ_][indexes[1:nloops]]
    vStride = [:strideX_,:strideY_,:strideZ_][indexes[1:nloops]]
    vNM = [:Nx_,:Ny_,:Nz_][indexes[1:nloops]]

    if platform == "cpu"
        if nloops == 0
            nothing
        elseif nloops == 1
            code =
                quote
                    Threads.@threads for $(vIc[1]) in 1:$(vNM[1])
                        $code
                    end
                end
        elseif nloops == 2
            code =
                quote
                    Threads.@threads for $(vIc[1]) in 1:$(vNM[1])
                        for $(vIc[2]) in 1:$(vNM[2])
                            $code
                        end
                    end
                end
        elseif nloops == 3
            code =
                quote
                    Threads.@threads for $(vIc[1]) in 1:$(vNM[1])
                        for $(vIc[2]) in 1:$(vNM[2])
                            for $(vIc[3]) in 1:$(vNM[3])
                                $code
                            end
                        end
                    end
                end
        end
    elseif platform == "gpu"
        code =
            if nloops == 0
                quote
                    indexX_ = (threadIdx().x) + (blockIdx().x - 1) * blockDim().x
                    strideX_ = blockDim().x * gridDim().x

                    indexY_ = (threadIdx().y) + (blockIdx().y - 1) * blockDim().y
                    strideY_ = blockDim().y * gridDim().y

                    indexZ_ = (threadIdx().z) + (blockIdx().z - 1) * blockDim().z
                    strideZ_ = blockDim().z * gridDim().z

                    if indexX_ == 1
                        $code
                    end
                end
            elseif nloops == 1
                code =
                    quote
                        indexX_ = (threadIdx().x) + (blockIdx().x - 1) * blockDim().x
                        strideX_ = blockDim().x * gridDim().x

                        indexY_ = (threadIdx().y) + (blockIdx().y - 1) * blockDim().y
                        strideY_ = blockDim().y * gridDim().y

                        indexZ_ = (threadIdx().z) + (blockIdx().z - 1) * blockDim().z
                        strideZ_ = blockDim().z * gridDim().z

                        for $(vIc[1]) in $(vIndex[1]):$(vStride[1]):$(vNM[1])
                            $code
                        end
                    end
            elseif nloops == 2
                code =
                    quote
                        indexX_ = (threadIdx().x) + (blockIdx().x - 1) * blockDim().x
                        strideX_ = blockDim().x * gridDim().x

                        indexY_ = (threadIdx().y) + (blockIdx().y - 1) * blockDim().y
                        strideY_ = blockDim().y * gridDim().y

                        indexZ_ = (threadIdx().z) + (blockIdx().z - 1) * blockDim().z
                        strideZ_ = blockDim().z * gridDim().z

                        for $(vIc[1]) in $(vIndex[1]):$(vStride[1]):$(vNM[1])
                            for $(vIc[2]) in $(vIndex[2]):$(vStride[2]):$(vNM[2])
                                $code
                            end
                        end
                    end
            elseif nloops == 3
                code =
                    quote
                        indexX_ = (threadIdx().x) + (blockIdx().x - 1) * blockDim().x
                        strideX_ = blockDim().x * gridDim().x

                        indexY_ = (threadIdx().y) + (blockIdx().y - 1) * blockDim().y
                        strideY_ = blockDim().y * gridDim().y

                        indexZ_ = (threadIdx().z) + (blockIdx().z - 1) * blockDim().z
                        strideZ_ = blockDim().z * gridDim().z

                        for $(vIc[1]) in $(vIndex[1]):$(vStride[1]):$(vNM[1])
                            for $(vIc[2]) in $(vIndex[2]):$(vStride[2]):$(vNM[2])
                                for ic3_ in $(vIndex[3]):$(vStride[3]):$(vNM[3])
                                    $code
                                end
                            end
                        end
                    end
            end
    else
        error("Platform should be or cpu or gpu. ", platform, " was given.")
    end

    return code
end

"""
    function wrapInFunction_(name, code)

Return **code** wrapped in a function called **name** and which returns nothing.
"""
function wrapInFunction_(name::Symbol, code::Expr)
    code =
        :(function $name(ARGS_)
            $code
            return nothing
        end
        )

    return code
end

"""
    function simpleLoopWrapInFunction_(platform, name, code)

Calls toguether functions simpleFirstLoop and wrapInFunction.
"""
function simpleFirstLoopWrapInFunction_(platform::String, name::Symbol, code::Expr)
    code = simpleFirstLoop_(platform,code)
    code =
        :(function $name(ARGS_)
            $code
            return nothing
        end
        )

    return code
end

"""
    function simpleGridLoopWrapInFunction_(platform, name, code)

Calls toguether functions simpleGridLoop and wrapInFunction.
"""
function simpleGridLoopWrapInFunction_(platform::String, name::Symbol, code::Expr, nLoop::Int; indexes = [1,2,3])
    code = simpleGridLoop_(platform,code,nLoop,indexes=indexes)
    code =
        :(function $name(ARGS_)
            $code
            return nothing
        end
        )

    return code
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
module RestartsTests

include("test_logparsing.jl")

end  # module RestartsTests

ext: jl


text:
using DigitSets
using Test
using Printf

const a = DigitSet(1, 2, 7)
const b = DigitSet(2, 5)

@test length(a) == 3
@test length(b) == 2
@test !isempty(a)
@test isempty(DigitSet())
@test collect(Int, a) == [1, 2, 7]
@test collect(Int, b) == [2, 5]
@test 1 in a
@test 2 in a
@test 7 in a
@test !(5 in a)
@test 2 in b
@test 5 in b
@test !(7 in b)
@test a == DigitSet(collect(a))
@test b == DigitSet(collect(b))
@test DigitSet(1:9) == DigitSet(1, 2, 3, 4, 5, 6, 7, 8, 9)

@test union(a, b) == DigitSet(1, 2, 5, 7)
@test union(a, b) == union(b, a)
@test intersect(a, b) == DigitSet(2)
@test intersect(a, b) == intersect(b, a)
@test setdiff(a, b) == DigitSet(1, 7)
@test setdiff(b, a) == DigitSet(5)
@test symdiff(a, b) == DigitSet(1, 5, 7)
@test symdiff(a, b) == symdiff(b, a)
@test union(DigitSet(), a) == a
@test union(a, DigitSet()) == a
@test intersect(a, DigitSet()) == DigitSet()
@test intersect(DigitSet(), a) == DigitSet()
@test setdiff(a, DigitSet()) == a
@test setdiff(DigitSet(), a) == DigitSet()
@test symdiff(a, DigitSet()) == a
@test symdiff(DigitSet(), a) == a

@test issubset(a, a)
@test a ⊆ a
@test issubset(b, b)
@test b ⊆ b
@test !issubset(a, b)
@test a ⊈ b
@test !issubset(b, a)
@test b ⊈ a
@test issubset(DigitSet(), a)
@test DigitSet() ⊆ a
@test issubset(DigitSet(), b)
@test DigitSet() ⊆ b
@test issubset(DigitSet(), DigitSet())
@test DigitSet() ⊆ DigitSet()
@test !issubset(a, DigitSet())
@test a ⊈ DigitSet()
@test !issubset(b, DigitSet())
@test b ⊈ DigitSet()
@test issubset(DigitSet(1), a)
@test DigitSet(1) ⊆ a
@test issubset(DigitSet(1, 2), a)
@test DigitSet(1, 2) ⊆ a
@test issubset(DigitSet(5), b)
@test DigitSet(5) ⊆ b
@test !(a < a)
@test !(a ⊊ a)
@test !(a < b)
@test !(a ⊊ b)
@test DigitSet(1, 2) < a
@test DigitSet(1, 2) ⊊ a

@test (@sprintf "%s" a) == "DigitSet(1,2,7)"
@test (@sprintf "%s" DigitSet()) == "DigitSet()"


text:
function test_diophantine_solve_diophantine_linear()
    @test solve(diophantine_equation_linear_xy()) == diophantine_anyx_anyy()
    @test solve(diophantine_equation_linear_xy(c0=1)) == diophantine_nonex_noney()

    @test solve(diophantine_equation_linear_xy(cy=2, c0=10)) == diophantine_anyx_oney(-5)
    @test solve(diophantine_equation_linear_xy(cy=2, c0=11)) == diophantine_nonex_noney()

    @test solve(diophantine_equation_linear_xy(cx=3, c0=-9)) == diophantine_onex_anyy(3)
    @test solve(diophantine_equation_linear_xy(cx=3, c0=2)) == diophantine_nonex_noney()

    eq = diophantine_equation_linear_xy(cx=10, cy=84, c0=16)
    sol = take(solve(eq), 1000) |> collect
    @test Int[evaluate(eq, s) for s in sol] == zeros(Int, 1000)

    @test solve(diophantine_equation_linear_xy(cx=10, cy=84, c0=-15)) == diophantine_nonex_noney()
end


text:
using ProbNumODE
using BSON
using DataFrames
using GaussianDistributions
using StructArrays

files = [
    "logistic.bson",
    "lotkavolterra.bson",
    "fitzhughnagumo.bson",
    "rigidbody.bson",
    "brusselator.bson",
    "logistic_fixedMLEMV.bson",
    "lotkavolterra_fixedMLEMV.bson",
    "fitzhughnagumo_fixedMLEMV.bson",
    "rigidbody_fixedMLEMV.bson",
    "brusselator_fixedMLEMV.bson",
]
files = "./data/" .* files

# r = results[collect(keys(results))[1]][1]
# config_keys = keys(r.config)
df = DataFrame(
    probname=String[],
    alg=String[],
    q=Int[],
    steprule=Symbol[],
    sigmarule=Symbol[],
    abstol=Real[],
    reltol=Real[],
    dt=Real[],
    smooth=Bool[],
    num_evals=Int[],
    nf=Int[],
    runtime=Real[],
    l∞=Real[],
    l2=Real[],
    L∞=Real[],
    L2=Real[],
    final=Real[],
    χ²=Real[],
    Χ²=Real[],
    final_χ²=Real[],
    final_est=Real[],
    l2_est=Real[],
    L2_est=Real[],
)

for file in files
    results = BSON.load(file)
    for (key, result_list) in results  # Iterate over the individual problems
        # probname = split(splitdir(file)[2], "_")[1]
        probname = splitext(basename(file))[1]
        if '_' in probname
            probname = split(probname, "_")[1]
        end
        for r in result_list
            config = r.config
            summaries = r.summaries
            # ssol = r.solution
            # final_vals = r.final_vals
            # errors = ProbNumODE.compute_errors(ssol)
            errors = summaries

            push!(df, (
                probname,
                config.alg,
                config.q,
                config.steprule,
                config.sigmarule,
                config.abstol,
                config.reltol,
                get(config, :dt, 0.0),
                config.smooth,
                summaries[:num_evals],
                summaries[:nf],
                summaries[:runtime],
                errors[:l∞],
                errors[:l2],
                errors[:L∞],
                errors[:L2],
                errors[:final],
                errors[:χ²],
                errors[:Χ²],
                errors[:final_χ²],
                errors[:final_est],
                errors[:l2_est],
                errors[:L2_est],
            ))
        end
    end
end

categorical!(df, :alg)
categorical!(df, :steprule)
categorical!(df, :sigmarule)

df.sigma_str = "sigma=" .* string.(Symbol.(df.sigmarule))
df.q_str = "IBM" .* string.(df.q)
df.step_str = "step=" .* string.(df.steprule)
df.smooth_str = "smooth=" .* string.(df.smooth)

@show df

max_stars_repo_name: nathanaelbosch/capos

max_issues_repo_name: nathanaelbosch/capos

max_forks_repo_name: nathanaelbosch/capos


text:
using FFTW
using Plots
using Roots
using Optim
using Distributions
using Random
using LinearAlgebra
using ProgressMeter
using Primes
using LaTeXStrings

"""
Nazarov's ``magic'' Phi function x*arctan(x)-(1/2)log(1+x*x)
"""
function nazarov_phi(x::Real)
    return x*atan(x) - 0.5*log(1+x*x)
end

"""
Derivative Phi'(x)
"""
function nazarov_phiprime(x::Real)
    return atan(x)
end

"""
Nazarov's I(f)
"""
function nazarov_potential(fv::Array)
    s = 0
    for x in fv
        s += nazarov_phi(x)
    end
    return s
end

"""
Real orthonormal basis for the DFT. We only return non-DC terms.
"""
function dft_mat(n::Integer)
    mat = zeros(n,n)
    q = 2*pi/n
    h_pt = ceil(Int, (n+1)/2)
    # ceil((n+1)/2) cosines, floor((n-1)/2) sines
    for i=1:h_pt
        for j=1:n
            mat[i,j] = sqrt(2)*cos(q*(i-1)*(j-1))
        end
    end
    for i=h_pt+1:n
        for j=1:n
            mat[i,j] = sqrt(2)*sin(q*(i-h_pt)*(j-1))
        end
    end
    mat[1,:] /= sqrt(2)
    if iseven(n)
        mat[h_pt,:] /= sqrt(2)
    end
    return mat
end

function nazarov_psi(n::Integer)
    return dft_mat(n)[2:end,:]
end

"""
Compute the Nazarov constant for the DFT. In the paper's notation, this is M(n).
The paper gives the asymptotics, here one can numerically study it.
"""
function nazarov_const(n::Integer)
    mat = dft_mat(n)
    nc = 3*pi/2
    max_rat = 0
    for i=1:n
        nrat = 1/norm(mat[i,:],1)
        if (nrat > max_rat)
            max_rat = nrat
        end
    end
    max_rat *= n
    return nc*max_rat*max_rat
end

function plot_nazarov_const(n1::Integer, n2::Integer)
    x = primes(n1,n2)
    plot(x, [nazarov_const(i) for i in x])
end

"""
Some functions for checking the l1 norm beta(n) for the DFT. Asymptotics are given
in the paper. Finite n behavior can be studied with these routines, basically
the above nazarov_const is governed by the minimal value of the l1 norm over the
divisors of n.
"""
function l1_cos(n::Integer)
    sum = 0
    omega = 2*pi/n
    for i=0:n-1
        sum += abs(cos(omega*i))
    end
    sum /= n
    return sum-2/pi
end

function plot_l1_cos(n1::Integer, n2::Integer)
    x = n1:n2
    y = [l1_cos(n) for n in x]
    p = plot(x,y)
    savefig(p, "/tmp/l1_dft.pdf")
end

"""
Seed the greedy local optimum search randomly, works well in practice.
"""
function nazarov_seed_cortege(n::Integer)
    return rand([1,-1], n)
end

"""
Perform the local optimum search, given a vector that describes the power control.
"""
function nazarov_search!(a::Array, psi::Array, cortege::Array{Int})
    n = length(cortege)
    vn = length(psi[1,:])
    iter_order = [i for i in 1:n]
    fv = sum([cortege[j]*a[j]*psi[j,:] for j=1:n])
    cur_pot = nazarov_potential(fv)
    iter = 0
    eval_ct = 1
    is_opt = false
    println("-----------------------------------------------------------------")
    println("Performing local optimum search...")
    println("(iter, pot)")
    while !is_opt
        iter += 1
        print("($iter, $(cur_pot/vn))")
        randperm!(iter_order)
        is_opt = true
        for j in iter_order
            fv -= 2*cortege[j]*a[j]*psi[j,:]
            new_pot = nazarov_potential(fv)
            eval_ct += 1
            if new_pot > cur_pot
                is_opt = false
                cur_pot = new_pot
                cortege[j] = -cortege[j]
                break
            end
            fv += 2*cortege[j]*a[j]*psi[j,:]
        end
    end
    println("\n---------------------------------------------------------------")
    println("Function evaluations (scaled by n^2): $(eval_ct/(n*n))")
    println("\n---------------------------------------------------------------")
end

"""
Perform the shift/scale to get a [0,1] sequence from the locally optimal cortege.
We do not do the scaling of the paper, as it is conservative and meant for theorems.
The scaling here works better in practice and is natural.
"""
function _nazarov_aperture(a::Array, psi::Array, cortege::Array{Int})
    n = length(cortege)
    vn = length(psi[1,:])
    fv = sum([cortege[j]*a[j]*psi[j,:] for j=1:n])
    aper = zeros(vn)
    for i=1:vn
        aper[i] = nazarov_phiprime(fv[i])
    end
    # flip signs if positive
    if (sum(aper) > 0)
        aper = -aper
    end

    min_val = minimum(aper)
    max_val = maximum(aper)
    # scaling is trivial if nonnegative (very unlikely!)
    if (min_val > 0)
        warn("Something really weird, everything is nonnegative a priori!")
        aper /= max_val
        return aper
    end
    # the standard case. Note that this scaling is not the same used in the proof,
    # as the one used in the proof is a conservative one for theoretical soundness,
    # and we can improve upon it in code. Note that this scaling does not guarantee
    # \rho <= 0.5, though it is very close to it.
    for i=1:vn
        aper[i] = (aper[i]-min_val)/(max_val-min_val)
    end
    return aper
end

"""
General function that hooks things together: just takes in an n, a power control
vector (not necessarily normalized), and calls the above routines to get an aperture.
"""
function nazarov_aperture(n::Integer, powers::Array)
    psi = nazarov_psi(n)
    a = powers/norm(powers)
    cortege = nazarov_seed_cortege(n-1)
    nazarov_search!(a, psi, cortege)
    #redirect_stdout(() -> nazarov_search!(a, psi, cortege), open("/dev/null", "w"))
    return _nazarov_aperture(a, psi, cortege)
end

"""
Sanity checks on the generated aperture, and computation of its spectrum
"""
function nazarov_ahat(n::Integer, powers::Array, ctx::FFTW.cFFTWPlan)
    aperture = nazarov_aperture(n, powers)
    @assert(minimum(aperture) >= 0)
    @assert(maximum(aperture) <= 1)
    ahat = ctx * aperture
    ahat = [abs(x) for x in ahat]
    rho = ahat[1]/n
    println("rho (Nazarov): $rho")
    return ahat
end

"""
Scene correlation functions, generated by sampling at n equally spaced points
in [0,1] of a nonnegative, bounded, continuous function with symmetry about 1/2.
"""

"""
We first define the underlying d(x), for 0 <= x <= 0.5. All are underscored to
signify their "internal" nature.
"""
function _exp_corr(x::Real, theta::Real, beta::Real)
    return theta*(beta^x)
end

function _iid_corr(x::Real, theta::Real)
    return theta
end

function _char_corr(x::Real, theta::Real, beta::Real, gamma::Real)
    if (x < beta - gamma)
        return theta
    elseif (x > beta + gamma)
        return 0
    else
        return theta*(beta + gamma - x)/(2*gamma)
    end
end

function _power_corr(x::Real, theta::Real, beta::Real, gamma::Real)
    return theta/(1+beta*(x^gamma))
end

"""
General d vector generator based on sampling and symmetry about 1/2.
"""
function gen_corr(corr_fn, n::Integer)
    f = x -> (x < 0.5) ? corr_fn(x) : corr_fn(1-x)
    d = (1/n)*ones(n)
    for i=1:n
        j = i-1
        d[i] *= f(j/n)
    end
    return d
end

"""
Specializations to various cases of interest.
"""
function exp_corr(theta::Real, beta::Real, n::Integer)
    return gen_corr(x -> _exp_corr(x, theta, beta), n)
end

function iid_corr(theta::Real, n::Integer)
    return gen_corr(x -> _iid_corr(x, theta), n)
end

function char_corr(theta::Real, beta::Real, gamma::Real, n::Integer)
    return gen_corr(x -> _char_corr(x, theta, beta, gamma), n)
end

function power_corr(theta::Real, beta::Real, gamma::Real, n::Integer)
    return gen_corr(x -> _power_corr(x, theta, beta, gamma), n)
end

"""
Computation of power loss factor, i.e using a suboptimal rho. Here, a
(corresponding to W/J) is fixed. rhos describes an array of wherever we have
or consider spectrally flat constructions.
"""
function _power_loss_factor!(a::Real, rhos::Array)
    num = 2*a - 2*sqrt(a*(a+1))+1
    max_denom = 0
    for x in rhos
        denom = x*(1-x)/(a+x)
        if denom > max_denom
            max_denom = denom
        end
    end
    return num/max_denom
end

"""
Here we optimize (i.e maximize) _power_loss_factor over a, with an array of rhos
given (e.g rhos=[0.5, 0.25]). One may somewhat painfully do this rigorously;
this is just a sanity check.
"""
function power_loss_factor(rhos)
    # The 10 is completely ad-hoc, but works fine (see the plot)
    x = range(0, stop=10, length=10^8)
    y = zeros(length(x))
    i = 1
    @showprogress for i=1:length(x)
        y[i] = _power_loss_factor(x[i], rhos)
    end
    println("Loss factor: $(maximum(y))")
    x = range(0, stop=10, length=10^3)
    y = [_power_loss_factor(a, rhos) for a in x]
    p = plot(x, y)
    savefig(p, "/tmp/power_loss.pdf")
end

"""
LMMSE expressions.
"""
function eval_mmse(d::Array, n::Integer, t::Real, w::Real, j::Real, ahat::Array)
    mmse = 0
    rho = ahat[1]/n
    gamma = t/(n*(w+rho*j))
    for i=1:n
        mmse += 1/((1/d[i])+gamma*ahat[i]*ahat[i])
    end
    return mmse
end

function eval_mmse_flat(d::Array, rho::Real, n::Integer, t::Real, w::Real, j::Real)
    mmse = 0
    r = n*rho
    gamma = t/(n*(w+rho*j))
    mmse += 1/((1/d[1])+gamma*r*r)
    r = sqrt((n*n*rho*(1-rho))/(n-1))
    for i=2:n
        mmse += 1/((1/d[i])+gamma*r*r)
    end
    return mmse
end

function eval_mmse_rand(d::Array, rho::Real, n::Integer, t::Real, w::Real, j::Real, ctx::FFTW.cFFTWPlan)
    a = zeros(n)
    pd = Bernoulli(rho)
    ntr = round(Int, 400) # 400 is arbitrary, larger will make plot smoother
    mmse_arr = zeros(ntr)
    gamma = t/(n*(w+rho*j))
    for i=1:ntr
        rand!(pd, a)
        ahat = ctx * a
        mmse = 0
        for i=1:n
            r = abs(ahat[i])
            mmse += 1/((1/d[i])+gamma*r*r)
        end
        mmse_arr[i] = mmse
    end
    mmse = sum(mmse_arr)/ntr
    return mmse
end

function eval_mmse_lens(d::Array, n::Integer, t::Real, w::Real, j::Real)
    mmse = 0
    for i=1:n
        mmse += 1/((1/d[i])+((t*n)/(w+j)))
    end
    return mmse
end

function _mmse_power_sum(c::Real, d::Array, gamma::Real)
    n = length(d)
    p = 0
    for i=2:n
        p += max(0, (c-(1/d[i]))/gamma)
    end
    return p
end

function frac(x::Real)
    return x-floor(x)
end

"""
Waterfilling computation for power control.
"""
function waterfill_powers!(powers::Array, d::Array, rho::Real, n::Integer, t::Real, w::Real, j::Real)
    gamma = t/(n*(w+rho*j))
    # A hack for t=0
    if (gamma == 0)
        for i=1:n-1
            j = i+1
            powers[i] = 1/sqrt(n-1)
        end
        return
    end
    r = n*rho
    p = n*(floor(r)+frac(r)*frac(r))-r*r
    ub = (1/d[1])+gamma*p
    c = find_zero(x -> _mmse_power_sum(x,d,gamma)-p, (0,1.5*ub))
    for i=1:n-1
        j = i+1
        powers[i] = max(0, (c-(1/d[j]))/gamma)
    end
end

"""
Lower bound on LMMSE based on waterfilling.
"""
function mmse_lb(d::Array, rho::Real, n::Integer, t::Real, w::Real, j::Real)
    gamma = t/(n*(w+rho*j))
    # A hack for t=0
    if (gamma == 0)
        return sum(d)
    end
    r = n*rho
    p = n*(floor(r)+frac(r)*frac(r))-r*r
    ub = (1/d[1])+gamma*p
    c = find_zero(x -> _mmse_power_sum(x,d,gamma)-p, (0,1.5*ub))
    mmse = 0
    p = n*n*rho*rho
    mmse += 1/((1/d[1])+gamma*p)
    for i=2:n
        p = max(0, (c-(1/d[i])))
        mmse += 1/((1/d[i])+p)
    end
    return mmse
end

"""
Read the values of n for which we have spectrally flat constructions from
data/ into a dictionary of n->rho_arr[n].
"""
function get_spec_flat()
    fdict = Dict("data/1_2_vals.txt" => 0.5, "data/1_4_vals.txt" => 0.25, "data/1_8_vals.txt" => 0.125)
    spec_flat_dict = Dict{Int64,Array{Float64,1}}()
    for (file, rho) in fdict
        for line in eachline(file)
            n = parse(Int64, line)
            if !haskey(spec_flat_dict, n)
                spec_flat_dict[n] = [rho]
            else
                push!(spec_flat_dict[n], rho)
            end
        end
    end
    return spec_flat_dict
end

"""
Finally, the function that generates the plots of the paper, saved in "/tmp/mmse.pdf".
Specifically, set globally:
julia> n=677; theta=1; w=0.001; j=0.001; start=-50; stop=50; res=256;

For plot (a), do:
julia> beta = 1; gamma=0; plot_lmmse(n,theta,beta,gamma,w,j,start,stop,res);
For plot (b), do:
julia> beta = 0.02, gamma=0.005; plot_lmmse(n,theta,beta,gamma,w,j,start,stop,res);
"""
function plot_lmmse(n::Integer, theta::Real, beta::Real, gamma::Real, w::Real, j::Real, start::Real=0, len::Real=10, res=256)
    d = char_corr(theta, beta, gamma, n)
    fft_ctx = plan_fft(zeros(n))
    waterfill = zeros(n-1)
    x = range(start, length=res, stop=len)
    y_lb = zeros(length(x))
    popt = zeros(length(x))
    y_prng = zeros(length(x))
    y_spec = zeros(length(x))
    y_nazarov = zeros(length(x))
    spec_flat_dict = get_spec_flat()
    if !haskey(spec_flat_dict, n)
        error("No spectrally flat construction available at n=$n")
    else
        rho_arr = spec_flat_dict[n]
    end
    i = 1
    @showprogress for t in x
        t /= 10
        t = 10^t
        t *= n
        # lower bound
        res = Optim.optimize(rho->mmse_lb(d,rho,n,t,w,j), 0, 1)
        y_lb[i] = 10*log10(Optim.minimum(res))
        popt[i] = Optim.minimizer(res)
        waterfill_powers!(waterfill,d,popt[i],n,t,w,j)
        # random sequence
        res = Optim.optimize(rho->eval_mmse_rand(d,rho,n,t,w,j,fft_ctx),0,1)
        y_prng[i] = 10*log10(Optim.minimum(res))
        # spectrally flat, rho in rho_arr
        y_spec[i] = 10*log10(minimum([eval_mmse_flat(d,rho,n,t,w,j) for rho in rho_arr]))
        # Nazarov construction
        ahat_nazarov = nazarov_ahat(n, waterfill, fft_ctx)
        y_nazarov[i] = 10*log10(eval_mmse(d,n,t,w,j,ahat_nazarov))
        i += 1
    end
    # force pyplot, GR has issues with LaTeX labels
    pyplot()
    xlabel_str = L"10\log_{10}(t/n)"
    p = plot(x, [y_lb y_prng y_spec y_nazarov],
             label=["lower bound" "optimal random on-off" "spectrally flat" "Nazarov"],
             xlabel=xlabel_str,
             ylabel="LMMSE (dB)")
    plot_str = "/tmp/mmse.pdf"
    savefig(p, plot_str)
    q = plot(x, popt,
             title="optimal p vs time for iid scene at\n (n,theta,w,j)=($n,$theta,$w,$j)",
             label="optimal p",
             xlabel=xlabel_str,
             ylabel="p")
    plot_str = "/tmp/optp_iid.pdf"
    savefig(q, plot_str)
end

max_stars_repo_name: gajjanag/apertures

max_issues_repo_name: gajjanag/apertures

max_forks_repo_name: gajjanag/apertures


text:
# Autogenerated wrapper script for AprilTags_jll for aarch64-linux-gnu
export libapriltag

JLLWrappers.@generate_wrapper_header("AprilTags")
JLLWrappers.@declare_library_product(libapriltag, "libapriltag.so.3")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libapriltag,
        "lib/libapriltag.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# %% imports

import Random
include("../../Modules/Gridcells/HexGrid.jl")

# %%

function grid_gen()
    Random.seed!(0)

    dilats = [2*sqrt(2), 2, sqrt(2), 1]

    grids = HexGrid[]

    for mod in 1:4
        n_neuron = 2^(mod+1)
        for neuron in 1:n_neuron
            res = 0.03*dilats[mod]
            node_size = 0.25*sqrt(3)*res*dilats[mod]
            theta = rand()*(2pi/6)-pi/6
            grid = HexGrid(res, theta = theta, node_size = node_size)
            push!(grids, grid)
        end
    end
    return grids
end

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
using Test
using ProxyInterfaces

# ProxyInterfaces.@dict
# ---------------------

struct MyDict{T}
  value::T
end

ProxyInterfaces.dict(::Type{MyDict{T}}) where T = T
ProxyInterfaces.dict(mydict::MyDict) = mydict.value
ProxyInterfaces.@dict MyDict

mydict = MyDict(Dict(:a => 4))

@test mydict[:a] == 4
@test length(mydict) == 1
@test (:a => 4) in mydict
@test haskey(mydict, :a)
@test !haskey(mydict, :b)
@test collect(values(mydict)) == [4]
@test collect(keys(mydict)) == [:a]

@test mydict.value == Dict(:a => 4)
@test merge(
  MyDict(Dict(:a => 4, :b => 5)),
  MyDict(Dict(:b => 32, :c => "hi"))
).value == Dict(:a => 4, :b => 32, :c => "hi")

@test_throws MethodError mydict[:b] = 42
@test_throws MethodError empty!(mydict)

# ProxyInterfaces.@dict_mutable
# -----------------------------

struct MyMutableDict{T}
  value::T
end

ProxyInterfaces.dict(::Type{MyMutableDict{T}}) where T = T
ProxyInterfaces.dict(mydict::MyMutableDict) = mydict.value
ProxyInterfaces.@dict_mutable MyMutableDict

mydict_mutable = MyMutableDict(Dict(:a => 4))
mydict_mutable[:b] = 42
@test mydict_mutable[:b] == 42

empty!(mydict_mutable)
@test mydict_mutable.value == Dict()


text:
"""
    restrict(refspace, element1, element2)

Computes the restriction of a set of local shape functions on `element1` as linear
combinations of the set of local shape functions on `element2`. More precisely `restrict`
returns an `NxM` matrix `P` such that the `i`-th local shape ``g_i`` function on element2
can be written as:

``g_i = sum_{j=1}^{M} P_{ij} f_j``
"""
function restrict end

function restrict(sp::Space, submesh::CompScienceMeshes.AbstractMesh)

    fns = similar(sp.fns)
    S = eltype(eltype(fns))
    for (i,fn) in pairs(sp.fns)
        fns[i] = Vector{S}()
        for sh in fn
            cell = CompScienceMeshes.restrict(submesh, sh.cellid)
            cell == 0 && continue
            push!(fns[i], S(cell, sh.refid, sh.coeff))
        end
    end

    similar(sp, submesh, fns, sp.pos)
end


text:
using CombinedParsers
using Test

@testset "CombinedParsers.jl" begin
    include("test-parser.jl")
    include("test-re.jl")
    include("pcretest-parser.jl")
end


text:
using Documenter
using LazyArrays

makedocs(
    sitename = "LazyArrays",
    modules = [LazyArrays],
    pages = [
        "Home" => "index.md",
        hide("internals.md"),
    ],
)

deploydocs(
    repo = "github.com/JuliaArrays/LazyArrays.jl"
)


text:
@testset "filter" begin
    inds = Indices([1,2,3,4,5])

    @test isequal(filter(iseven, inds)::Indices, Indices([2, 4]))
    @test isequal(filter(isodd, inds)::Indices, Indices([1, 3, 5]))

    @test isequal(filterview(iseven, inds)::Dictionaries.FilteredIndices, Indices([2, 4]))
    @test isequal(filterview(isodd, inds)::Dictionaries.FilteredIndices, Indices([1, 3, 5]))
    @test filterview(iseven, inds)[2] == 2
    @test length(filterview(iseven, inds)) == 2
    @test_throws IndexError filterview(isodd, inds)[2]
    @test in(2, filterview(iseven, inds))
    @test !in(2, filterview(isodd, inds))

    @test collect(reverse(filterview(iseven, inds))) == [4, 2]

    @test empty(filterview(iseven, inds))::Indices == Indices{Int}()
    @test keys(similar(filterview(iseven, inds))::Dictionary) == Indices{Int}([2, 4])

    filter!(iseven, inds)
    @test isequal(inds, Indices([2, 4]))

    d = Dictionary([1,2,3,4,5], [1,3,2,4,5])

    @test isequal(filter(iseven, d)::Dictionary, dictionary([3=>2, 4=>4]))
    @test isequal(filter(isodd, d)::Dictionary, dictionary([1=>1, 2=>3, 5=>5]))

    @test isequal(filterview(iseven, d)::Dictionaries.FilteredDictionary, dictionary([3=>2, 4=>4]))
    @test isequal(filterview(isodd, d)::Dictionaries.FilteredDictionary, dictionary([1=>1, 2=>3, 5=>5]))
    @test filterview(iseven, d)[3] == 2
    @test length(filterview(iseven, d)) == 2
    @test_throws IndexError filterview(isodd, d)[3]

    @test empty(filterview(iseven, d))::Dictionary{Int, Int} == Dictionary{Int, Int}()
    @test keys(similar(filterview(iseven, d))::Dictionary) == Indices{Int}([3, 4])

    filter!(iseven, d)
    @test isequal(d, Dictionary([3,4],[2,4]))

    for _ in 1:100
        @test rand(filterview(iseven, d)) in [2,4]
    end
    @test_throws ArgumentError rand(filterview(x -> false, d))
end


text:
"""
r = _create_grid_Greens_function(number_of_cells, resolution)\n
Computes the radius for the domain given by number_of_cells.\n
This radius is used to compute the weakened Greens functions.\n
"""
function _create_grid_Greens_function(number_of_cells, resolution)
    #get max dimensions
    dx,dy,dz = 2 .^(ceil.(Int, log2.(number_of_cells.+2).+1))
    lx = resolution[1].*collect(-dx/2f0:dx/2f0-1)
    ly = resolution[2].*collect(-dy/2f0:dy/2f0-1)
    lz = resolution[3].*collect(-dz/2f0:dz/2f0-1)

    # repeat the axes along the number corresponding dimensions
    X = repeat(lx, outer=(1, dy, dz))
    Y = repeat(reshape(ly,1,:), outer=(dx, 1, dz))
    Z = repeat(reshape(lz,1,1,:), outer=(dx, dy, 1))

    radius = @. sqrt(X^2 + Y^2 + Z^2)

    return radius
end

"""
Internal function that computes the greens function.
"""
function _create_Greens_function(radius, k_b, a)
    # sommerfeld radiation condition
    scalar_G = 3/(k_b*a)^2*(sin(k_b*a)/(k_b*a)-cos(k_b*a))
    # Green's function at r != 0
    G = @. exp.(-1im*k_b*radius)/(4*Float32(π)*radius)*scalar_G
    # Green's function at r = 0
    G₀ = 3/(4*Float32(π)*k_b^2*a^3)*((1+1im*k_b*a)*exp(-1im*k_b*a)-1)
    G = circshift(G, floor.(Int, size(G)./2))
    G[1,1,1] = G₀

    return fft(G)
end

"""
G = create_Greens_function(number_of_cells, resolution, k_b).\n
creates the green function on the larger grid.\n
"""
function create_Greens_function(number_of_cells, resolution, k_b)
    a = findmin(resolution)[1]/2
    radius = _create_grid_Greens_function(number_of_cells, resolution)
    G = _create_Greens_function(radius, k_b, a)

    return prod(resolution).*G
end

"""
_Greens_function_times_contrast_source!(vector_potential_1d, G, dielectric, E, w, vector_potential_3d, Ig, fft_tmp, planfft, planifft)\n
Using the Green's function and contrast source we compute the vector potential.\n
a = ∫GdielectricEdV is computed using an FFT.\n
or f = ∫GKdV is computed using an FFT.\n
"""
function _Greens_function_times_contrast_source!(
    vector_potential_1d,
    G,
    dielectric,
    electric_field,
    vector_potential_3d,
    Ig,
    fft_tmp,
    planfft,
    planifft
)
    # move E to larger domain for fft
    mul!(fft_tmp, Transpose(Ig), electric_field)
    # multiply with the contrast
    vector_potential_3d .= dielectric.*reshape(fft_tmp, size(G)..., 3)
    # do the fft of the contrast source (dielectric.*E)
    planfft*vector_potential_3d
    # multiply the green's function with the contrast source in freq domain
    @. vector_potential_3d = G*vector_potential_3d
    planifft*vector_potential_3d
    # map to smaller domain again
    mul!(vector_potential_1d, Ig, vec(vector_potential_3d))

    return nothing
end

"""
source_to_incident_field!(source_to_field, vector_potential_1d, vector_potential_3d, G, dielectric, Ig, operator, fft_tmp, planfft, planifft, scaling)\n
Compute the incident field in a vacuum given a source distribution\n
source_to_field will contain the initial source distribution and will contain the incident field when the function is finished\n
Use when electric source to electric field or magnetic source to magnetic field.\n
For an electric current density distribution to incident electric field:\n
    source_to_field: J\n
    vector_potential_1d: a\n
    vector_potential_3d: A\n
    operator: AtoE\n
    scaling: 1/(1im*ω*ϵ₀)\n
For a magnetic source distribution to incident magnetic field:\n
    source_to_field: K\n
    vector_potential_1d: f\n
    vector_potential_3d: F\n
    operator: FtoH\n
    scaling: 1/(1im*ω*μ₀)\n
source_to_incident_field!(incident_field, source, vector_potential_1d, vector_potential_3d, G, dielectric, Ig, operator, fft_tmp, planfft, planifft, scaling)\n
Use when electric source to magnetic field or magnetic source to electric field.\n
For a magnetic source distribution to incident electric field:\n
    source_to_field: K\n
    vector_potential_1d: f\n
    vector_potential_3d: F\n
    operator: FtoE\n
    scaling: 1/(1im*ω*μ₀)\n
For an electric current density distribution to incident magnetic field:\n
    source_to_field: J\n
    vector_potential_1d: a\n
    vector_potential_3d: A\n
    operator: AtoH\n
    scaling: 1/(1im*ω*ϵ₀)\n
"""
function source_to_incident_field!(
    source_to_field,
    vector_potential_1d,
    vector_potential_3d,
    G,
    dielectric,
    Ig,
    operator,
    fft_tmp,
    planfft,
    planifft,
    scaling
)
    _Greens_function_times_contrast_source!(vector_potential_1d, G, dielectric.vac, source_to_field, vector_potential_3d, Ig, fft_tmp, planfft, planifft)
    mul!(source_to_field, operator, vector_potential_1d)
    @. source_to_field *= scaling

    return nothing
end

function source_to_incident_field!(
    incident_field,
    source,
    vector_potential_1d,
    vector_potential_3d,
    G,
    dielectric,
    Ig,
    operator,
    fft_tmp,
    planfft,
    planifft,
    scaling
)
    _Greens_function_times_contrast_source!(vector_potential_1d, G, dielectric.vac, source, vector_potential_3d, Ig, fft_tmp, planfft, planifft)
    mul!(incident_field, operator, vector_potential_1d)
    @. incident_field *= scaling

    return nothing
end

"""
total_minus_scattered_electric_field!(v,total_electric_field,tmp,a,w,A,G,dielectric,Ig,AToE,fft_tmp,planfft,planifft)\n
calculate the incident electric field from the total electric field.\n
This can also be seen as Ax in the context of ||b-Ax||_2^2\n
"""
_total_minus_scattered_electric_field!(v, total_electric_field, vector_potential_1d, vector_potential_3d, G, dielectric, Ig, AToE, fft_tmp, planfft, planifft) = _total_minus_scattered_electric_field!(v, total_electric_field, vector_potential_1d, vector_potential_3d, G, dielectric, Ig, AToE, fft_tmp, planfft, planifft, 1)
function _total_minus_scattered_electric_field!(
    v,
    total_electric_field,
    vector_potential_1d,
    vector_potential_3d,
    G,
    dielectric,
    Ig,
    AToE,
    fft_tmp,
    planfft,
    planifft,
    scalar
)
    _Greens_function_times_contrast_source!(vector_potential_1d, G, dielectric.patient, total_electric_field, vector_potential_3d, Ig, fft_tmp, planfft, planifft)
    copyto!(v, total_electric_field)
    mul!(v, AToE, vector_potential_1d, -scalar, scalar)

    return nothing
end

"""
_total_minus_scattered_electric_field_for_residual!(v,total_electric_field,tmp,vector_potential_3d,w,vector_potential_3d,G,dielectric,Ig,AToE,fft_tmp,planfft,planifft,scalar)\n
similar to total_minus_scattered_electric_field!(), but now r-AΔx is computed\n
That is the change in the residual as a result of the update in x\n
"""
function _total_minus_scattered_electric_field_for_residual!(
    v,
    total_electric_field,
    vector_potential_1d,
    vector_potential_3d,
    G,
    dielectric,
    Ig,
    AToE,
    fft_tmp,
    planfft,
    planifft,
    scalar
)
    _Greens_function_times_contrast_source!(vector_potential_1d, G, dielectric.patient, total_electric_field, vector_potential_3d, Ig, fft_tmp, planfft, planifft)
    mul!(total_electric_field, AToE, vector_potential_1d, -scalar, scalar)
    axpy!(one(ComplexF32), total_electric_field, v)

    return nothing
end

"""
_set_source_values_in_domain!(sources, incident_electric_field, v, S)\n
will set the current density values to -source_values\n
Using the S matrix these are placed into the incident electric field.\n
"""
function _set_source_values_in_domain!(sources, incident_electric_field, source_values, S)
    # set source values to the vector we are multiplying A with
    sources .= -one(ComplexF32).*source_values
    # set the source values into the incident electric field (will be equal to Jinc before the function JIncToEinc!)
    mul!(incident_electric_field, S, sources)

    return nothing
end

"""
library_matrix_times_current_density!(current_density, nUpdates, sources, incident_electric_field, S, tmp, fft_tmp, G, vector_potential_3d, w, dielectric, vector_potential_3d, AToE, Ig, scaling, pfft, pifft, x, p, r, rt, u, v, q, uq).\n
calculate the scattered electric field given the current density.\n
"""
function library_matrix_times_current_density!(
    current_density,
    sources,
    incident_electric_field,
    S,
    fft_tmp,
    G,
    vector_potential_3d,
    dielectric,
    vector_potential_1d,
    AToE,
    Ig,
    scaling,
    pfft,
    pifft,
    x,
    p,
    r,
    rt,
    u,
    v,
    q,
    uq
)
    source_values_input = -current_density

    _set_source_values_in_domain!(sources, incident_electric_field, source_values_input, S)
    #from source in E incident to the actual incident electric field
    source_to_incident_field!(incident_electric_field, vector_potential_1d, vector_potential_3d, G, dielectric, Ig, AToE, fft_tmp, pfft, pifft, scaling)
    #compute total electric field
    cgs_efield!(incident_electric_field, fft_tmp, G, vector_potential_3d, dielectric, vector_potential_1d, AToE, Ig, pfft, pifft, x, p, r, rt, u, v, q, uq)

    return x
end

"""
library_matrix_times_current_density_gpu!(current_density, nUpdates, sources, incident_electric_field, S, tmp, fft_tmp, G, vector_potential_3d, w, dielectric, vector_potential_3d, AToE, Ig, scaling, pfft, pifft, x, p, r, rt, u, v, q, uq).\n
calculate the scattered electric field given the current density.\n
"""
function library_matrix_times_current_density_gpu!(
    current_density,
    sources,
    incident_electric_field,
    S,
    fft_tmp,
    G,
    vector_potential_3d,
    dielectric,
    vector_potential_1d,
    AToE,
    Ig,
    scaling,
    pfft,
    pifft,
    x,
    p,
    r,
    rt,
    u,
    v,
    q,
    uq
)
    source_values_input = -current_density |> cu

    _set_source_values_in_domain!(sources, incident_electric_field, source_values_input, S)
    #from source in E incident to the actual incident electric field
    source_to_incident_field!(incident_electric_field, vector_potential_1d, vector_potential_3d, G, dielectric, Ig, AToE, fft_tmp, pfft, pifft, scaling)
    #compute total electric field
    cgs_efield!(incident_electric_field, fft_tmp, G, vector_potential_3d, dielectric, vector_potential_1d, AToE, Ig, pfft, pifft, x, p, r, rt, u, v, q, uq)

    return x
end #TODO: check if cpu version/collect are required

"""
_volume_integral_equation_wrapper!(output, sources, incident_electric_field, C, input_vector, S, tmp, fft_tmp, G, vector_potential_3d, w, dielectric, vector_potential_3d, AToE, Ig, scaling, pfft, pifft, x, p, r, rt, u, v, q, uq)\n
In place version to calculate Av = (I-CSZ)v = v - CSZv.\n
-CSZv is calculated using the VIE method by setting the sources to -v and multiplying the edges with the corresponding C value.
"""
function _volume_integral_equation_wrapper!(
    output,
    sources,
    incident_electric_field,
    C,
    input_vector,
    S,
    fft_tmp,
    G,
    vector_potential_3d,
    dielectric,
    vector_potential_1d,
    AToE,
    Ig,
    scaling,
    pfft,
    pifft,
    x,
    p,
    r,
    rt,
    u,
    v,
    q,
    uq
)
    # set sources in the vacuum to the correct values
    _set_source_values_in_domain!(sources, incident_electric_field, input_vector, S)
    # calculate the incident electric field inside the vacuum
    source_to_incident_field!(incident_electric_field, vector_potential_1d, vector_potential_3d, G, dielectric, Ig, AToE, fft_tmp, pfft, pifft, scaling)
    # solve for the field inside the patient anatomy
    cgs_efield!(incident_electric_field, fft_tmp, G, vector_potential_3d, dielectric, vector_potential_1d, AToE, Ig, pfft, pifft, x, p, r, rt, u, v, q, uq)
    # calculate the (v - CSZv) vector for the minimization process
    mul!(output, Transpose(S), x)
    @. output *= C
    # output = input_vector - CSZinput_vector
    axpy!(one(ComplexF32), input_vector, output)

    return nothing
end

"""
total_minus_scattered_electric_field_map!(v, total_electric_field, vector_potential_3d, vector_potential_3d, G, dielectric, Ig, AToE, fft_tmp, planfft, planifft, scalar).\n
used as a map to solve for the total electric field.\n
"""
function total_minus_scattered_electric_field_map!(
    v,
    total_electric_field,
    vector_potential_1d,
    vector_potential_3d,
    G,
    dielectric,
    Ig,
    AToE,
    fft_tmp,
    planfft,
    planifft,
    scalar
)
    copyto!(v, total_electric_field)
    _Greens_function_times_contrast_source!(vector_potential_1d, G, dielectric.patient, v, vector_potential_3d, Ig, fft_tmp, planfft, planifft)
    mul!(v, AToE, vector_potential_1d, -scalar, scalar)

    return v
end


text:
###############
# Status type #
###############

@ghdef mutable struct Status
    id::Union{Int, Nothing}
    total_count::Union{Int, Nothing}
    state::Union{String, Nothing}
    description::Union{String, Nothing}
    context::Union{String, Nothing}
    sha::Union{String, Nothing}
    url::Union{HTTP.URI, Nothing}
    target_url::Union{HTTP.URI, Nothing}
    created_at::Union{Dates.DateTime, Nothing}
    updated_at::Union{Dates.DateTime, Nothing}
    creator::Union{Owner, Nothing}
    repository::Union{Repo, Nothing}
    statuses::Union{Vector{Status}, Nothing}
end

Status(id::Real) = Status(Dict("id" => id))

namefield(status::Status) = status.id

###############
# API Methods #
###############

@api_default function create_status(api::GitHubAPI, repo, sha; options...)
    result = gh_post_json(api, "/repos/$(name(repo))/statuses/$(name(sha))"; options...)
    return Status(result)
end

@api_default function statuses(api::GitHubAPI, repo, ref; options...)
    results, page_data = gh_get_paged_json(api, "/repos/$(name(repo))/commits/$(name(ref))/statuses"; options...)
    return map(Status, results), page_data
end

@api_default function status(api::GitHubAPI, repo, ref; options...)
    result = gh_get_json(api, "/repos/$(name(repo))/commits/$(name(ref))/status"; options...)
    return Status(result)
end


text:
# --------------------------------------------------------------------
struct FitResult
    timestamp::DateTime
    elapsed::Float64
    ndata::Int
    nfree::Int
    dof::Int
    fitstat::Float64
    gofstat::Float64
    log10testprob::Float64
    resid::Vector{Float64}
    mzer::Union{Nothing, AbstractMinimizerStatus}
end

function fit!(model::Model, data::Measures{N};
              minimizer=lsqfit(), dry=false) where N
    timestamp = now()
    evaluate!(model)

    data1d = flatten(data)
    resid1d = fill(NaN, length(data1d))
    params = model.meval.params[model.meval.ifree]
    @assert length(params) > 0 "No free parameter in the model"

    function private_func(pvalues::Vector{Float64})
        model.meval.pvalues[model.meval.ifree] .= pvalues
        eval2!(model)
        eval3!(model)
        resid1d .= (model() .- data1d.val) ./ data1d.unc

        evaluate_showvalues(x) = () -> begin
            dof = (length(resid1d) - length(pvalues))
            [(:fit_stat, sum(abs2, x) / dof)]
        end
        ProgressMeter.next!(prog; showvalues=evaluate_showvalues(resid1d))
        return resid1d
    end

    dof = length(resid1d) - length(params)
    prog = ProgressUnknown("Model (dof=$dof) evaluations:", dt=0.5, showspeed=true)
    if !dry
        result = minimize(minimizer, private_func, model.meval.params[model.meval.ifree])
        private_func(result.best)
        eval4!(model, result.unc)
    else
        private_func(getfield.(params, :val))
        eval4!(model, fill(NaN, length(model.meval.ifree)))
    end
    ProgressMeter.finish!(prog)

    # Prepare output
    ndata = length(resid1d)
    nfree = length(params)
    dof = ndata - nfree
    fitstat = sum(abs2, resid1d)
    gofstat = sum(abs2, resid1d)
    elapsed = now() - timestamp
    @assert isa(elapsed, Millisecond)
    return FitResult(timestamp, elapsed.value / 1.e3,
                     ndata, nfree, dof, fitstat, gofstat,
                     logccdf(Chisq(dof), gofstat) * log10(exp(1)),
                     resid1d, (dry  ?  nothing  :  result))
end

function fit!(multi::MultiModel, data::Vector{Measures{N}};
              minimizer=lsqfit(), dry=false) where N
    timestamp = now()
    evaluate!(multi)

    data1d = [flatten(data[id]) for id in 1:length(multi)]
    resid1d = fill(NaN, sum(length.(data1d)))
    params = Vector{Parameter}()
    for id in 1:length(multi)
        append!(params, multi[id].meval.params[multi[id].meval.ifree])
    end
    @assert length(params) > 0 "No free parameter in the model"

    function private_func(pvalues::Vector{Float64})
        i1 = 1
        for id in 1:length(multi)
            nn = length(multi[id].meval.ifree)
            i2 = i1 + nn - 1
            multi[id].meval.pvalues[multi[id].meval.ifree] = pvalues[i1:i2]
            i1 += nn
        end
        eval2!(multi)
        eval3!(multi)

        i1 = 1
        for id in 1:length(multi)
            nn = length(data1d[id])
            i2 = i1 + nn - 1
            resid1d[i1:i2] .= (multi[id]() .- data1d[id].val) ./ data1d[id].unc
            i1 += nn
        end
        evaluate_showvalues(x) = () -> begin
            dof = (length(resid1d) - length(pvalues))
            [(:fit_stat, sum(abs2, x) / dof)]
        end
        ProgressMeter.next!(prog; showvalues=evaluate_showvalues(resid1d))
        return resid1d
    end

    dof = length(resid1d) - length(params)
    prog = ProgressUnknown("Model (dof=$dof) evaluations:", dt=0.5, showspeed=true)
    if !dry
        result = minimize(minimizer, private_func, params)
        private_func(result.best)
        i1 = 1
        for id in 1:length(multi)
            nn = length(multi[id].meval.ifree)
            i2 = i1 + nn - 1
            eval4!(multi[id], result.unc[i1:i2])
            i1 += nn
        end
    else
        private_func(getfield.(params, :val))
        i1 = 1
        for id in 1:length(multi)
            nn = length(multi[id].meval.ifree)
            i2 = i1 + nn - 1
            eval4!(multi[id], fill(NaN, i2-i1+1))
            i1 += nn
        end
    end
    ProgressMeter.finish!(prog)

    # Prepare output
    ndata = length(resid1d)
    nfree = length(params)
    dof = ndata - nfree
    fitstat = sum(abs2, resid1d)
    gofstat = sum(abs2, resid1d)
    elapsed = now() - timestamp
    @assert isa(elapsed, Millisecond)
    return FitResult(timestamp, elapsed.value / 1.e3,
                     ndata, nfree, dof, fitstat, gofstat,
                     logccdf(Chisq(dof), gofstat) * log10(exp(1)),
                     resid1d, (dry  ?  nothing  :  result))
end

function free_param_names(model::Model)
    out = Vector{String}()
    for (cname, ceval) in model.cevals
        for (pid, par) in ceval.params
            if (!par.fixed)  &&  (model.cevals[cname].cfixed == 0)
                if pid.index == 0
                    push!(out, "[:$(cname)].$(pid.name)")
                else
                    push!(out, "[:$(cname)].$(pid.name)[$(pid.index)]")
                end
            end
        end
    end
    return out
end

function free_param_names(multi::MultiModel)
    out = Vector{String}()
    for id in 1:length(multi)
        append!(out, "[$(string(id))]" .* free_param_names(multi[id]))
    end
    return out
end

function print_param_covariance(model::Union{Model, MultiModel}, fitres::FitResult;
                                param::Union{Nothing, String}=nothing, sort=false, thresh=0.)
    @assert isa(fitres.mzer.specific, CMPFit.Result)
    if isa(model, Model)
        @assert isnothing(model.parent)
    end
    names = free_param_names(model)
    @assert length(names)^2 == length(fitres.mzer.specific.covar)
    ii = Vector{Int}()
    jj = Vector{Int}()
    covar = Vector{Float64}()
    for i in 1:length(names)
        for j in i+1:length(names)
            push!(covar, fitres.mzer.specific.covar[i, j])
            push!(ii, i)
            push!(jj, j)
        end
    end
    if sort
        ii    = ii[   sortperm(abs.(covar))]
        jj    = jj[   sortperm(abs.(covar))]
        covar = covar[sortperm(abs.(covar))]
    end
    for i in 1:length(ii)
        if !isnothing(param)
            (names[ii[i]] != param)  &&  continue
        end
        (abs(covar[i]) < thresh)  &&  continue
        @printf "%-30s  %-30s  %10.4f\n" names[ii[i]] names[jj[i]] covar[i]
    end
end


text:
using Primes

function totient_func(n)
    prime_fact = factor(n)
    for (p, k) in prime_fact
        n *= (1-1/p)
    end
    return n
end

function find_increment()
    p = primes(30)
    for i in 1:length(p)
	x = prod((1-1/prime) for prime in p[1:i])
	n = prod(p[1:i])
	if x*n/(n-1) < 15499/94744
	    return prod(p[1:i-1])
	end
    end
end

function find_denominator(inc)
    n = 0
    while true
        n += inc
        x = totient_func(n)
        if x/(n-1) < 15499/94744
            return n
        end
    end
end

function main()
    inc = find_increment()
    denom = find_denominator(inc)
    println(denom)
end

main()

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function genbo2grad(foldername, filename, u, xij, qi, qj, ti, tj, ai, aj, rho, mu, eta, kappa, natomtype)

# if (u[1]==0)
#     m = 0;
# else
#     n,m = size(u);
#     if (n==0) || (n == 1)
#         m = 0;
#     elseif (n == natomtype+1)
#         if u[n,1] == 0
#             m = 0;
#         end
#     else
#         error("Number of functions for the two-body bond order potential must be equal to 1 plus the number of atom types");
#     end
# end

if length(u)>0
    if (u[1]==0)
        m = 0;
        n = 0;
    else
        n,m = size(u);
        # if (n==0) || (n == 1)
        #     m = 0;
        # elseif (n == natomtype+1)
        #     if u[end,1] == 0
        #         m = 0;
        #     end
        # else
        #     error("Number of functions for the two-body bond order potential must be equal to 1 plus the number of atom types");
        # end
    end
else
    n = 0;
    m = 0;
end

gen = 1;
ifile = 0;
stropu = "";
strcpu = "";
strgpu = "";
potnum = 0;
fn = "PaircDensityGradient";

for i = 1:m
    ui = u[:,i];
    if (ui[n] != 0)
        potnum = potnum + 1;
        stropu1, strcpu1, strgpu1 = genpairgrad(foldername, filename*string(potnum), ui[1:(n-1)], xij, qi, qj, ti, tj, ai, aj, mu, eta, kappa, gen, ifile);
        stropu1 = modifystr(stropu1, natomtype);
        strcpu1 = modifystr(strcpu1, natomtype);
        strgpu1 = modifystr(strgpu1, natomtype);
        strgpu1 = replace(strgpu1, "\t\ti += blockDim.x * gridDim.x;\n\t\t}" => "\t\t}\n\t\ti += blockDim.x * gridDim.x;");
        stropu2, strcpu2, strgpu2 = gendensitygrad(foldername, fn*string(potnum), ui[n], rho, mu, eta, kappa, gen, ifile);
        stropu = stropu * stropu1 * "\n" * stropu2 * "\n";
        strcpu = strcpu * strcpu1 * "\n" * strcpu2 * "\n";
        strgpu = strgpu * strgpu1 * "\n" * strgpu2 * "\n";
    end
end

if potnum==0
    ifile = 0;
    gen = 0;
    stropu, strcpu, strgpu = genpairgrad(foldername, filename, u, xij, qi, qj, ti, tj, ai, aj, mu, eta, kappa, gen, ifile);
    tmpopu, tmpcpu, tmpgpu = gendensitygrad(foldername, fn, u, mu, mu, eta, kappa, gen, ifile);
    stropu = stropu  * "\n" * tmpopu * "\n";
    strcpu = strcpu  * "\n" * tmpcpu * "\n";
    strgpu = strgpu  * "\n" * tmpgpu * "\n";
    cppfiles(foldername, filename, stropu, strcpu, strgpu, 1);
else
    sp0, sp1, sp2, sp3 = getpotstr(2);
    sp0 = replace(sp0, "int ng)" => "int ng, int potnum)");
    sp1 = replace(sp1, "int);" => "int, int);");
    sp2 = replace(sp2, "int);" => "int, int);");
    tmpopu, tmpcpu, tmpgpu = genpotnum(filename, potnum, sp0, sp1, sp2, sp3);
    tmpopu = modstr11(tmpopu);
    tmpcpu = modstr11(tmpcpu);
    tmpgpu = modstr11(tmpgpu);
    stropu = stropu * tmpopu * "\n";
    strcpu = strcpu * tmpcpu * "\n";
    strgpu = strgpu * tmpgpu * "\n";

    sp0, sp1, sp2, sp3 = getpotstr(0);
    sp0 = replace(sp0, "int ng)" => "int ng, int potnum)");
    sp1 = replace(sp1, "int);" => "int, int);");
    sp2 = replace(sp2, "int);" => "int, int);");
    tmpopu, tmpcpu, tmpgpu = genpotnum(fn, potnum, sp0, sp1, sp2, sp3);
    tmpopu = modstr12(tmpopu);
    tmpcpu = modstr12(tmpcpu);
    tmpgpu = modstr12(tmpgpu);
    stropu = stropu * tmpopu;
    strcpu = strcpu * tmpcpu;
    strgpu = strgpu * tmpgpu;

    cppfiles(foldername, filename, stropu, strcpu, strgpu, 0);
end

end

function modifystr(oldstr, natomtype)

newstr = oldstr;
newstr = replace(newstr, "u_xij[" => "\tu_xij[");
for i = 1:natomtype
    s = "u[" * string(i-1) * " * i*" * string(natomtype) * "]";
    if contains(newstr, s)
        news = "if (tj[i] == " * string(i) * ") \n\t\t\t" * "u[i]";
        if (i>1)
            news = "}\n \t\t" * news;
        end
        newstr = replace(newstr, s => news);
    end
end
newstr = replace(newstr, "\t}\n}" => "\t\t}\n\t}\n}");

return newstr

end

function modstr11(tmpgpu)

tmpgpu = replace(tmpgpu, "(T *u," => "(T *u => T *u_xij,");
tmpgpu = replace(tmpgpu, "(u," => "(u, u_xij,");
tmpgpu = replace(tmpgpu, "(double *" => "(double *, double *");
tmpgpu = replace(tmpgpu, "(float *" => "(float *, float *");

return tmpgpu

end

function modstr12(tmpgpu)

tmpgpu = replace(tmpgpu, "(T *u," => "(T *u, T *u_rho,");
tmpgpu = replace(tmpgpu, "(u," => "(u, u_rho,");
tmpgpu = replace(tmpgpu, "(double *" => "(double *, double *");
tmpgpu = replace(tmpgpu, "(float *" => "(float *, float *");

return tmpgpu

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
import DataFrames
import ModelSanitizer
import Test

a = DataFrames.DataFrame()

a[:x] = [1, 2, 3, 4, 5, 6, 7, 8, 9]
a[:y] = [10, 20, 30, 40, 50, 60, 70, 80, 90]
a[:z] = [100, 200, 300, 400, 500, 600, 700, 800, 900]

Test.@test a[1, :x] == 1
Test.@test a[2, :x] == 2
Test.@test a[3, :x] == 3
Test.@test a[4, :x] == 4
Test.@test a[5, :x] == 5
Test.@test a[6, :x] == 6
Test.@test a[7, :x] == 7
Test.@test a[8, :x] == 8
Test.@test a[9, :x] == 9

Test.@test a[1, :y] == 10
Test.@test a[2, :y] == 20
Test.@test a[3, :y] == 30
Test.@test a[4, :y] == 40
Test.@test a[5, :y] == 50
Test.@test a[6, :y] == 60
Test.@test a[7, :y] == 70
Test.@test a[8, :y] == 80
Test.@test a[9, :y] == 90

Test.@test a[1, :z] == 100
Test.@test a[2, :z] == 200
Test.@test a[3, :z] == 300
Test.@test a[4, :z] == 400
Test.@test a[5, :z] == 500
Test.@test a[6, :z] == 600
Test.@test a[7, :z] == 700
Test.@test a[8, :z] == 800
Test.@test a[9, :z] == 900

ModelSanitizer.sanitize!(ModelSanitizer.Model(a), ModelSanitizer.Data(a))

for column in names(a)
    Test.@test all(a[:, column] .== 0)
end

Test.@test a[1, :x] == 0
Test.@test a[2, :x] == 0
Test.@test a[3, :x] == 0
Test.@test a[4, :x] == 0
Test.@test a[5, :x] == 0
Test.@test a[6, :x] == 0
Test.@test a[7, :x] == 0
Test.@test a[8, :x] == 0
Test.@test a[9, :x] == 0

Test.@test a[1, :y] == 0
Test.@test a[2, :y] == 0
Test.@test a[3, :y] == 0
Test.@test a[4, :y] == 0
Test.@test a[5, :y] == 0
Test.@test a[6, :y] == 0
Test.@test a[7, :y] == 0
Test.@test a[8, :y] == 0
Test.@test a[9, :y] == 0

Test.@test a[1, :z] == 0
Test.@test a[2, :z] == 0
Test.@test a[3, :z] == 0
Test.@test a[4, :z] == 0
Test.@test a[5, :z] == 0
Test.@test a[6, :z] == 0
Test.@test a[7, :z] == 0
Test.@test a[8, :z] == 0
Test.@test a[9, :z] == 0


text:
"""
Type: StoppingMeta

Methods: no methods.

Attributes:
- atol: absolute tolerance.
- rtol: relative tolerance.
- optimality0: optimality score at the initial guess.
- tol_check: Function of `atol`, `rtol` and `optimality0` testing a score to zero.
- tol_check_neg: Function of `atol`, `rtol` and `optimality0` testing a score to zero.
- check_pos: pre-allocation for positive tolerance
- check_neg: pre-allocation for negative tolerance
- recomp_tol: true if tolerances are updated
- optimality_check: a stopping criterion via an admissibility function
- unbounded_threshold: threshold for unboundedness of the problem.
- unbounded_x: threshold for unboundedness of the iterate.
- max_f: maximum number of function (and derivatives) evaluations.
- max_cntrs: Dict contains the maximum number of evaluations
- max_eval:  maximum number of function (and derivatives) evaluations.
- max_iter: threshold on the number of stop! call/number of iteration.
- max_time: time limit to let the algorithm run.
- nb\\_of\\_stop: keep track of the number of stop! call/iteration.
- start_time: keep track of the time at the beginning.
- fail\\_sub\\_pb: status.
- unbounded: status.
- unbounded_pb: status.
- tired: status.
- stalled: status.
- iteration_limit: status.
- resources: status.
- optimal: status.
- infeasible: status.
- main_pb: status.
- domainerror: status.
- suboptimal: status.
- stopbyuser: status.
- exception: status.

- meta_user_struct:  Any
- user_check_func!: Function (AbstractStopping, Bool) -> callback.

`StoppingMeta(;atol :: Number = 1.0e-6, rtol :: Number = 1.0e-15, optimality0 :: Number = 1.0, tol_check :: Function = (atol,rtol,opt0) -> max(atol,rtol*opt0), tol_check_neg :: Function = (atol,rtol,opt0) -> -max(atol,rtol*opt0), unbounded_threshold :: Number = 1.0e50, unbounded_x :: Number = 1.0e50, max_f :: Int = typemax(Int), max_eval :: Int = 20000, max_iter :: Int = 5000, max_time :: Number = 300.0, start_time :: Float64 = NaN, meta_user_struct :: Any = nothing, kwargs...)`

an alternative with constant tolerances:

`StoppingMeta(tol_check :: T, tol_check_neg :: T;atol :: Number = 1.0e-6, rtol :: Number = 1.0e-15, optimality0 :: Number = 1.0, unbounded_threshold :: Number = 1.0e50, unbounded_x :: Number = 1.0e50, max_f :: Int = typemax(Int), max_eval :: Int = 20000, max_iter :: Int = 5000, max_time :: Number = 300.0, start_time :: Float64 = NaN, meta_user_struct :: Any = nothing, kwargs...)`

Note:
- It is a mutable struct, therefore we can modify elements of a `StoppingMeta`.
- The `nb_of_stop` is incremented everytime `stop!` or `update_and_stop!` is called
- The `optimality0` is modified once at the beginning of the algorithm (`start!`)
- The `start_time` is modified once at the beginning of the algorithm (`start!`)
      if not precised before.
- The different status: `fail_sub_pb`, `unbounded`, `unbounded_pb`, `tired`, `stalled`,
      `iteration_limit`, `resources`, `optimal`, `main_pb`, `domainerror`, `suboptimal`, `infeasible`
- `fail_sub_pb`, `suboptimal`, and `infeasible` are modified by the algorithm.
- `optimality_check` takes two inputs (`AbstractNLPModel`, `NLPAtX`)
 and returns a `Number` or an `AbstractVector` to be compared to `0`.
- `optimality_check` does not necessarily fill in the State.

Examples: `StoppingMeta()`, `StoppingMeta(1., -1.)`
"""
mutable struct StoppingMeta{TolType <: Number,
                            CheckType, #Type of the tol_check output
                            MUS, #Meta User Struct
                            IntType <: Int
                            } <: AbstractStoppingMeta

  # problem tolerances
  atol                :: TolType # absolute tolerance
  rtol                :: TolType # relative tolerance
  optimality0         :: TolType # value of the optimality residual at starting point
  tol_check           :: Union{Function, CheckType} #function of atol, rtol and optimality0
                                  #by default: tol_check = max(atol, rtol * optimality0)
                                  #other example: atol + rtol * optimality0
  tol_check_neg       :: Union{Function, CheckType} # function of atol, rtol and optimality0
  check_pos           :: CheckType #pre-allocation for positive tolerance
  check_neg           :: CheckType #pre-allocation for negative tolerance
  optimality_check    :: Function # stopping criterion
                                  # Function of (pb, state; kwargs...)
                                  #return type  :: Union{Number, eltype(stp.meta)}
  recomp_tol          :: Bool #true if tolerances are updated

  unbounded_threshold :: TolType # beyond this value, the problem is declared unbounded
  unbounded_x         :: TolType # beyond this value, ||x||_\infty is unbounded

  # fine grain control on ressources
  max_f               :: IntType    # max function evaluations allowed TODO: used?
  max_cntrs           :: Dict{Symbol, Int} #contains the detailed max number of evaluations

  # global control on ressources
  max_eval            :: IntType    # max evaluations (f+g+H+Hv) allowed TODO: used?
  max_iter            :: IntType    # max iterations allowed
  max_time            :: Float64 # max elapsed time allowed

  #intern Counters
  nb_of_stop          :: IntType
  #intern start_time
  start_time          :: Float64

  # stopping properties status of the problem)
  fail_sub_pb         :: Bool
  unbounded           :: Bool
  unbounded_pb        :: Bool
  tired               :: Bool
  stalled             :: Bool
  iteration_limit     :: Bool
  resources           :: Bool
  optimal             :: Bool
  infeasible          :: Bool
  main_pb             :: Bool
  domainerror         :: Bool
  suboptimal          :: Bool
  stopbyuser          :: Bool
  exception           :: Bool

  meta_user_struct    :: MUS
  user_check_func!    :: Function #called dans Stopping._user_check!(stp, x)

end

function StoppingMeta(tol_check           :: CheckType,
                      tol_check_neg       :: CheckType;
                      atol                :: Number   = 1.0e-6,
                      rtol                :: Number   = 1.0e-15,
                      optimality0         :: Number   = 1.0,
                      optimality_check    :: Function = (a,b) -> Inf,
                      recomp_tol          :: Bool     = false,
                      unbounded_threshold :: Number   = 1.0e50, #typemax(Float64)
                      unbounded_x         :: Number   = 1.0e50,
                      max_f               :: Int      = typemax(Int),
                      max_cntrs           :: Dict{Symbol, Int} = Dict{Symbol, Int}(),
                      max_eval            :: Int      = 20000,
                      max_iter            :: Int      = 5000,
                      max_time            :: Float64  = 300.0,
                      start_time          :: Float64  = NaN,
                      meta_user_struct    :: Any      = nothing,
                      user_check_func!    :: Function = (stp :: AbstractStopping, start :: Bool) -> nothing,
                      kwargs...) where CheckType

  check_pos = tol_check
  check_neg = tol_check_neg

  # This might be an expansive step.
  # if (true in (check_pos .< check_neg)) #any(x -> x, check_pos .< check_neg)
  #     throw(ErrorException("StoppingMeta: tol_check should be greater than tol_check_neg."))
  # end

  fail_sub_pb     = false
  unbounded       = false
  unbounded_pb    = false
  tired           = false
  stalled         = false
  iteration_limit = false
  resources       = false
  optimal         = false
  infeasible      = false
  main_pb         = false
  domainerror     = false
  suboptimal      = false
  stopbyuser      = false
  exception       = false

  nb_of_stop = 0

  #new{TolType, typeof(check_pos), typeof(meta_user_struct)}
  return StoppingMeta(atol, rtol, optimality0,
                      tol_check, tol_check_neg,
                      check_pos, check_neg, optimality_check, recomp_tol,
                      unbounded_threshold, unbounded_x,
                      max_f, max_cntrs, max_eval, max_iter,
                      max_time, nb_of_stop, start_time,
                      fail_sub_pb, unbounded, unbounded_pb, tired, stalled,
                      iteration_limit, resources, optimal, infeasible, main_pb,
                      domainerror, suboptimal, stopbyuser, exception,
                      meta_user_struct, user_check_func!)
end

function StoppingMeta(;atol               :: Number   = 1.0e-6,
                      rtol                :: Number   = 1.0e-15,
                      optimality0         :: Number   = 1.0,
                      tol_check           :: Function = (atol :: Number, rtol :: Number, opt0 :: Number) -> max(atol,rtol*opt0),
                      tol_check_neg       :: Function = (atol :: Number, rtol :: Number, opt0 :: Number) -> - tol_check(atol,rtol,opt0),
                      optimality_check    :: Function = (a,b) -> Inf,
                      recomp_tol          :: Bool     = true,
                      unbounded_threshold :: Number   = 1.0e50, #typemax(Float64)
                      unbounded_x         :: Number   = 1.0e50,
                      max_f               :: Int      = typemax(Int),
                      max_cntrs           :: Dict{Symbol, Int} = Dict{Symbol, Int}(),
                      max_eval            :: Int      = 20000,
                      max_iter            :: Int      = 5000,
                      max_time            :: Float64  = 300.0,
                      start_time          :: Float64  = NaN,
                      meta_user_struct    :: Any      = nothing,
                      user_check_func!    :: Function = (stp :: AbstractStopping, start :: Bool) -> nothing,
                      kwargs...)

  check_pos = tol_check(atol, rtol, optimality0)
  check_neg = tol_check_neg(atol, rtol, optimality0)

  # This might be an expansive step.
  # if (true in (check_pos .< check_neg)) #any(x -> x, check_pos .< check_neg)
  #     throw(ErrorException("StoppingMeta: tol_check should be greater than tol_check_neg."))
  # end

  fail_sub_pb     = false
  unbounded       = false
  unbounded_pb    = false
  tired           = false
  stalled         = false
  iteration_limit = false
  resources       = false
  optimal         = false
  infeasible      = false
  main_pb         = false
  domainerror     = false
  suboptimal      = false
  stopbyuser      = false
  exception       = false

  nb_of_stop = 0

  #new{TolType, typeof(check_pos), typeof(meta_user_struct)}
  return StoppingMeta(atol, rtol, optimality0,
                      tol_check, tol_check_neg,
                      check_pos, check_neg, optimality_check, recomp_tol,
                      unbounded_threshold, unbounded_x,
                      max_f, max_cntrs, max_eval, max_iter,
                      max_time, nb_of_stop, start_time,
                      fail_sub_pb, unbounded, unbounded_pb, tired, stalled,
                      iteration_limit, resources, optimal, infeasible, main_pb,
                      domainerror, suboptimal, stopbyuser, exception,
                      meta_user_struct, user_check_func!)
end

const meta_statuses = [:fail_sub_pb,
                       :unbounded,
                       :unbounded_pb,
                       :tired,
                       :stalled,
                       :iteration_limit,
                       :resources,
                       :optimal,
                       :suboptimal,
                       :main_pb,
                       :domainerror,
                       :infeasible,
                       :stopbyuser,
                       :exception]

"""
`OK_check(meta :: StoppingMeta)`

Return true if one of the decision boolean is true.
"""
function OK_check(meta :: StoppingMeta{TolType, CheckType, MUS, IntType}
                  ) where {TolType, CheckType, MUS, IntType}
  #13 checks
  OK = meta.optimal || meta.tired || meta.iteration_limit || meta.resources || meta.unbounded || meta.unbounded_pb || meta.main_pb || meta.domainerror || meta.suboptimal || meta.fail_sub_pb || meta.stalled || meta.infeasible || meta.stopbyuser
  return OK
end

"""
`tol_check(meta :: StoppingMeta)`

Return the pair of tolerances, recomputed if `meta.recomp_tol` is `true`.
"""
function tol_check(meta :: StoppingMeta{TolType, CheckType, MUS, IntType}
                  ) where {TolType, CheckType, MUS, IntType}

  if meta.recomp_tol
    atol, rtol, opt0 = meta.atol, meta.rtol, meta.optimality0
    setfield!(meta, :check_pos, meta.tol_check(atol, rtol, opt0))
    setfield!(meta, :check_neg, meta.tol_check_neg(atol, rtol, opt0))
  end

  return (meta.check_pos, meta.check_neg)
end

"""
`update_tol!(meta :: StoppingMeta; atol = meta.atol, rtol = meta.rtol, optimality0 = meta.optimality0)`

Update the tolerances parameters. Set `meta.recomp_tol` as `true`.
"""
function update_tol!(meta     :: StoppingMeta{TolType, CheckType, MUS, IntType};
                     atol        :: TolType = meta.atol,
                     rtol        :: TolType = meta.rtol,
                     optimality0 :: TolType = meta.optimality0
                     ) where {TolType, CheckType, MUS, IntType}

  setfield!(meta, :recomp_tol, true)
  setfield!(meta, :atol, atol)
  setfield!(meta, :rtol, rtol)
  setfield!(meta, :optimality0, optimality0)

  return meta
end

function reinit!(meta :: StoppingMeta{TolType, CheckType, MUS, IntType}
                ) where {TolType, CheckType, MUS, IntType}

  for k in meta_statuses
    setfield!(meta, k, false)
  end

  return meta
end

function checktype(meta :: StoppingMeta{TolType, CheckType, MUS, IntType}
                   ) where {TolType, CheckType, MUS, IntType}
  return CheckType
end

function toltype(meta :: StoppingMeta{TolType, CheckType, MUS, IntType}
                   ) where {TolType, CheckType, MUS, IntType}
  return TolType
end

function metausertype(meta :: StoppingMeta{TolType, CheckType, MUS, IntType}
                   ) where {TolType, CheckType, MUS, IntType}
  return MUS
end

function inttype(meta :: StoppingMeta{TolType, CheckType, MUS, IntType}
                   ) where {TolType, CheckType, MUS, IntType}
  return IntType
end

import Base.show
function show(io :: IO, meta :: AbstractStoppingMeta)
  varlines = "$(typeof(meta)) has"
  if OK_check(meta)
    ntrue = 0
    for f in meta_statuses
      if getfield(meta, f) && ntrue == 0
        ntrue += 1
        varlines=string(varlines, @sprintf(" %s", f))
      elseif getfield(meta, f)
        ntrue += 1
        varlines=string(varlines, @sprintf(",\n %s", f))
      end
    end
    varlines= ntrue == 1 ? string(varlines, " as only true status.\n") : string(varlines, " as true statuses.\n")
  else
    varlines=string(varlines, " now no true statuses.\n")
  end
  varlines=string(varlines,  "The return type of tol check functions is $(checktype(meta))")
  if meta.recomp_tol
    varlines=string(varlines, ", and these functions are reevaluated at each stop!.\n")
  else
    varlines=string(varlines, ".\n")
  end
  varlines=string(varlines,  "Current tolerances are: \n")
  for k in [:atol, :rtol, :optimality0, :unbounded_threshold, :unbounded_x, :max_f, :max_eval, :max_iter, :max_time]
    varlines=string(varlines,  @sprintf("%19s: %s (%s) \n", k, getfield(meta, k), typeof(getfield(meta, k))))
  end
  if metausertype(meta) != Nothing
    varlines=string(varlines,  "The user defined structure in the meta is a $(metausertype(meta)).\n")
  else
    varlines=string(varlines,  "There is no user defined structure in the meta.\n")
  end
  println(io, varlines)
end


text:
function make_gui_html(cur_id, cur_widgets)
  cur_script = """

  """

  display(HTML(cur_html * cur_script))

  cur_id
end


text:
# This script will build a `sysimage` for Windows.
# Please submit a pull-request if you modify it to support additional operating systems.

# Load the libraries
using NeXLSpectrum, Gadfly, DataFrames

# Load the package compiler
using PackageCompiler

# This specifies which libraries to include and also specifies a "precompile.jl" file which forces
# the just-in-time compiler to compile frequently used functions.
PackageCompiler.create_sysimage(
    [ "Gadfly", "DataFrames", "BoteSalvatICX", "FFAST", "NeXLUncertainties", "NeXLCore", "NeXLMatrixCorrection", "NeXLSpectrum" ];
    sysimage_path=joinpath(homedir(), ".julia", "NeXLSysimage.dll"),
    precompile_execution_file=joinpath(pkgdir(NeXLSpectrum),"scripts", "precompile.jl"))

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
# import Pkg
# Pkg.add("Gadfly")
#
# ^ you have to fight with system at the moment because of Compact conflict.
# I edited Manifest to remove specific reference to Compat version then added
# Gadfly, which downgrades Compat to 2.2

using TreeParzen

import Gadfly

######### PCHOICE #############
example_space = Dict(
    :example => HP.PChoice(
        :example,
        [
            Prob(0.1, 0),
            Prob(0.2, 1),
            Prob(0.7, 2),
        ]
    )
)

trials = [ask(example_space) for i in 1:1000]
samples = getindex.(getproperty.(trials, :hyperparams), :example)
vals = sort(unique(samples))
counts = sum(samples .== vals'; dims=1)
probs = dropdims(counts'/sum(counts), dims=2)
p = Gadfly.plot(x=vals, y=probs, Gadfly.Geom.hair, Gadfly.Geom.point, Gadfly.Scale.y_continuous(minvalue=0.0), Gadfly.Guide.xticks(ticks=vals), Gadfly.Guide.yticks(ticks=Float64.(0:0.1:1)));

p |> Gadfly.SVGJS("$(@__DIR__)/../hp_images/pchoice.svg")
######### UNIFORM #############
example_space = Dict(
    :example => HP.Uniform(:example, 0.0, 1.0),
)

trials = [ask(example_space) for i in 1:1000]
samples = getindex.(getproperty.(trials, :hyperparams), :example)
p = Gadfly.plot(x=samples, Gadfly.Stat.density(bandwidth=0.05), Gadfly.Geom.polygon(fill=true, preserve_order=true));

p |> Gadfly.SVGJS("$(@__DIR__)/../hp_images/uniform.svg")
######### QUNIFORM #############
example_space = Dict(
    :example => HP.QuantUniform(:example, 0.0, 10.0, 2.0),
)

trials = [ask(example_space) for i in 1:1000]
samples = getindex.(getproperty.(trials, :hyperparams), :example)
vals = sort(unique(samples))
counts = sum(samples .== vals'; dims=1)
probs = dropdims(counts'/sum(counts), dims=2)
p = Gadfly.plot(x=vals, y=probs, Gadfly.Geom.hair, Gadfly.Geom.point, Gadfly.Scale.y_continuous(minvalue=0.0), Gadfly.Guide.xticks(ticks=vals));

p |> Gadfly.SVGJS("$(@__DIR__)/../hp_images/quniform.svg")
######### NORMAL #############
example_space = Dict(
    :example => HP.Normal(:example, 4.0, 5.0),
)

trials = [ask(example_space) for i in 1:1000]
samples = getindex.(getproperty.(trials, :hyperparams), :example)
p = Gadfly.plot(x=samples, Gadfly.Stat.density(bandwidth=1), Gadfly.Geom.polygon(fill=true, preserve_order=true));

p |> Gadfly.SVGJS("$(@__DIR__)/../hp_images/normal.svg")
######### QNORMAL #############
example_space = Dict(
    :example => HP.QuantNormal(:example, 2., 0.5, 1.0),
)

trials = [ask(example_space) for i in 1:1000]
samples = getindex.(getproperty.(trials, :hyperparams), :example)
vals = sort(unique(samples))
counts = sum(samples .== vals'; dims=1)
probs = dropdims(counts'/sum(counts), dims=2)
p = Gadfly.plot(x=vals, y=probs, Gadfly.Geom.hair, Gadfly.Geom.point, Gadfly.Scale.y_continuous(minvalue=0.0), Gadfly.Guide.xticks(ticks=vals));

p |> Gadfly.SVGJS("$(@__DIR__)/../hp_images/qnormal.svg")
######### LOGNORMAL #############
example_space = Dict(
    :example => HP.LogNormal(:example, log(3.0), 1.0),
)

trials = [ask(example_space) for i in 1:1000]
samples = getindex.(getproperty.(trials, :hyperparams), :example)
p = Gadfly.plot(x=samples, Gadfly.Stat.density(bandwidth=1.0), Gadfly.Geom.polygon(fill=true, preserve_order=true), Gadfly.Scale.x_log10, Gadfly.Guide.xlabel("x (log)"));

p |> Gadfly.SVGJS("$(@__DIR__)/../hp_images/lognormallog.svg")

p = Gadfly.plot(x=samples, Gadfly.Stat.density(bandwidth=0.5), Gadfly.Geom.polygon(fill=true, preserve_order=true));

p |> Gadfly.SVGJS("$(@__DIR__)/../hp_images/lognormal.svg")
######### QLOGNORMAL #############
example_space = Dict(
    :example => HP.LogQuantNormal(:example, log(3.0), 0.5, 2.0),
)

trials = [ask(example_space) for i in 1:1000]
samples = getindex.(getproperty.(trials, :hyperparams), :example)
vals = sort(unique(samples))
counts = sum(samples .== vals'; dims=1)
probs = dropdims(counts'/sum(counts), dims=2)
p = Gadfly.plot(x=vals, y=probs, Gadfly.Geom.hair, Gadfly.Geom.point, Gadfly.Scale.y_continuous(minvalue=0.0), Gadfly.Guide.xticks(ticks=vals));

p |> Gadfly.SVGJS("$(@__DIR__)/../hp_images/qlognormal.svg")
######### LOGUNIFORM #############
example_space = Dict(
    :example => HP.LogUniform(:example, log(1.0), log(5.0)),
)

trials = [ask(example_space) for i in 1:1000]
samples = getindex.(getproperty.(trials, :hyperparams), :example)
p = Gadfly.plot(x=samples, Gadfly.Stat.density(bandwidth=0.25), Gadfly.Geom.polygon(fill=true, preserve_order=true));

p |> Gadfly.SVGJS("$(@__DIR__)/../hp_images/loguniform.svg")
######### QLOGUNIFORM #############
example_space = Dict(
    :example => HP.QuantLogUniform(:example, log(1.0), log(5.0), 1.0),
)

trials = [ask(example_space) for i in 1:1000]
samples = getindex.(getproperty.(trials, :hyperparams), :example)
vals = sort(unique(samples))
counts = sum(samples .== vals'; dims=1)
probs = dropdims(counts'/sum(counts), dims=2)
p = Gadfly.plot(x=vals, y=probs, Gadfly.Geom.hair, Gadfly.Geom.point, Gadfly.Scale.y_continuous(minvalue=0.0), Gadfly.Guide.xticks(ticks=vals));

p |> Gadfly.SVGJS("$(@__DIR__)/../hp_images/qloguniform.svg")
######### THEEND #############