text:
"""Integrates ShallowWaters forward in time."""
function time_integration(  Prog::PrognosticVars{Tprog},
                            Diag::DiagnosticVars{T,Tprog},
                            S::ModelSetup{T,Tprog}) where {T<:AbstractFloat,Tprog<:AbstractFloat}

    @unpack u,v,η,sst = Prog
    @unpack u0,v0,η0 = Diag.RungeKutta
    @unpack u1,v1,η1 = Diag.RungeKutta
    @unpack du,dv,dη = Diag.Tendencies
    @unpack um,vm = Diag.SemiLagrange

    @unpack dynamics,RKo,tracer_advection = S.parameters
    @unpack RKaΔt,RKbΔt = S.constants

    @unpack nt,dtint = S.grid
    @unpack nstep_advcor,nstep_diff,nadvstep,nadvstep_half = S.grid

    if dynamics == "linear"
        Ix!(Diag.VolumeFluxes.h_u,S.forcing.H)
        Iy!(Diag.VolumeFluxes.h_v,S.forcing.H)
    end

    # propagate initial conditions
    copyto!(u0,u)
    copyto!(v0,v)
    copyto!(η0,η)

    # feedback, output initialisation and storing initial conditions
    feedback = feedback_init(S)
    netCDFfiles = NcFiles(feedback,S)
    output_nc!(0,netCDFfiles,Prog,Diag,S)

    nans_detected = false
    t = 0           # model time
    for i = 1:nt

        # ghost point copy for boundary conditions
        ghost_points!(u,v,η,S)
        copyto!(u1,u)
        copyto!(v1,v)
        copyto!(η1,η)

        # Runge-Kutta 4th order / 3rd order
        for rki = 1:RKo
            if rki > 1
                ghost_points!(u1,v1,η1,S)
            end

            # type conversion for mixed precision
            u1rhs = convert(Diag.PrognosticVarsRHS.u,u1)
            v1rhs = convert(Diag.PrognosticVarsRHS.v,v1)
            η1rhs = convert(Diag.PrognosticVarsRHS.η,η1)

            rhs!(u1rhs,v1rhs,η1rhs,Diag,S,t)

            if rki < RKo
                caxb!(u1,u,RKbΔt[rki],du)   #u1 .= u .+ RKb[rki]*Δt*du
                caxb!(v1,v,RKbΔt[rki],dv)   #v1 .= v .+ RKb[rki]*Δt*dv
                caxb!(η1,η,RKbΔt[rki],dη)   #η1 .= η .+ RKb[rki]*Δt*dη
            end

            # sum RK-substeps on the go
            axb!(u0,RKaΔt[rki],du)          #u0 .+= RKa[rki]*Δt*du
            axb!(v0,RKaΔt[rki],dv)          #v0 .+= RKa[rki]*Δt*dv
            axb!(η0,RKaΔt[rki],dη)          #η0 .+= RKa[rki]*Δt*dη
        end

        ghost_points!(u0,v0,η0,S)

        # type conversion for mixed precision
        u0rhs = convert(Diag.PrognosticVarsRHS.u,u0)
        v0rhs = convert(Diag.PrognosticVarsRHS.v,v0)
        η0rhs = convert(Diag.PrognosticVarsRHS.η,η0)

        # ADVECTION and CORIOLIS TERMS
        # although included in the tendency of every RK substep,
        # only update every nstep_advcor steps if nstep_advcor > 0
        if dynamics == "nonlinear" && nstep_advcor > 0 && (i % nstep_advcor) == 0
            advection_coriolis!(u0rhs,v0rhs,η0rhs,Diag,S)
        end

        # DIFFUSIVE TERMS - SEMI-IMPLICIT EULER
        # use u0 = u^(n+1) to evaluate tendencies, add to u0 = u^n + rhs
        # evaluate only every nstep_diff time steps
        if (i % nstep_diff) == 0
            bottom_drag!(u0rhs,v0rhs,η0rhs,Diag,S)
            diffusion!(u0rhs,v0rhs,Diag,S)
            add_drag_diff_tendencies!(u0,v0,Diag,S)
        end

        t += dtint

        # TRACER ADVECTION
        tracer!(i,u0rhs,v0rhs,Prog,Diag,S)

        # feedback and output
        feedback.i = i
        feedback!(Prog,feedback,S)
        output_nc!(i,netCDFfiles,Prog,Diag,S)

        if feedback.nans_detected
            break
        end

        # RK3/4 copy back from substeps
        copyto!(u,u0)
        copyto!(v,v0)
        copyto!(η,η0)
    end

    # finalise feedback and output
    feedback_end!(feedback)
    output_close!(netCDFfiles,feedback,S)

    return PrognosticVars{Tprog}(remove_halo(u,v,η,sst,S)...)
end

"""Add to a x multiplied with b. a += x*b """
function axb!(a::Array{T,2},x::T,b::Array{T,2}) where {T<:AbstractFloat}
    m,n = size(a)
    @boundscheck (m,n) == size(b) || throw(BoundsError())

    #TODO @simd?
    @inbounds for j ∈ 1:n
        for i ∈ 1:m
           a[i,j] += x*b[i,j]
        end
    end
end

"""c equals add a to x multiplied with b. c = a + x*b """
function caxb!(c::Array{T,2},a::Array{T,2},x::T,b::Array{T,2}) where {T<:AbstractFloat}
    m,n = size(a)
    @boundscheck (m,n) == size(b) || throw(BoundsError())
    @boundscheck (m,n) == size(c) || throw(BoundsError())

    #TODO @simd?
    @inbounds for j ∈ 1:n
        for i ∈ 1:m
           c[i,j] = a[i,j] + x*b[i,j]
        end
    end
end

"""Convert function for two arrays, X1, X2, in case their eltypes differ.
Convert every element from X1 and store it in X2."""
function Base.convert(X2::Array{T2,N},X1::Array{T1,N}) where {T1,T2,N}

    @boundscheck size(X2) == size(X1) || throw(BoundsError())

    @inbounds for i in eachindex(X1)
            X2[i] = T2(X1[i])
    end

    return X2
end

"""Convert function for two arrays, X1, X2, in case their eltypes are identical.
Just pass X1, such that X2 is pointed to the same place in memory."""
function Base.convert(X2::Array{T,N},X1::Array{T,N}) where {T,N}
    @boundscheck size(X2) == size(X1) || throw(BoundsError())
    return X1
end


text:
# This file includes several general functions which helps the solver.
# This includes:
# - threaded matrix multiplication (matmul_threaded)
# - symmetrizing an ArbRefMatrix (symmetric!)
# - inner product for ArbRefMatrices & BlockDiagonals (LinearAlgebra.dot)
# - Arblib Cholesky without error bounds (approx_cholesky!)

# extending LinearAlgebra.dot for the BlockDiagonal matrices and Arb(Ref)Matrices.
# In principle dot(A,B) is also (mathematically) defined if they have different block sizes, and can be done more efficient than materalizing the matrices.
# But in that case it is far more difficult to implement, so we just give the normal matrix results then
function LinearAlgebra.dot(A::BlockDiagonal, B::BlockDiagonal)
    # tot = Arb(prec=precision(A))
    # for (a,b) in zip(blocks(A),blocks(B))
    #     Arblib.add!(tot,tot,dot(a,b))
    # end
    # return tot
    if length(blocks(A)) == length(blocks(B))
        return sum(dot(a, b) for (a, b) in zip(blocks(A), blocks(B)))
    else
        return dot(Matrix(A),Matrix(B))
    end
end

#Extending the LinearAlgebra.dot for ArbMatrices. Arblib does not have an inner product for matrices
function LinearAlgebra.dot(A::T, B::T) where T <: Union{ArbMatrix,ArbRefMatrix}
    #This is also faster for ArbMatrices due to the addmul instead of += A[i,j]*B[i,j], which allocates
    @assert size(A) == size(B)
    res = Arb(0, prec = precision(A))
    for i = 1:size(A, 1)
        for j = 1:size(A, 2) #Arb is row major?
            Arblib.addmul!(res, A[i,j], B[i,j])
        end
    end
    return res
end

function LinearAlgebra.transpose(A::T;prec=precision(A)) where T<:Union{ArbRefMatrix,ArbMatrix}
    B = T(size(A,2),size(A,1),prec=prec)
    Arblib.transpose!(B,A)
    return B
end

function symmetric!(A::ArbRefMatrix,uplo=:U)
    #make the matrix symmetric, using the upper triangular part (uplo=:U) or the lower triangular part otherwise
    size(A,1) == size(A,2) || error("A must be square")

    for j=1:size(A,2)
        for i=j+1:size(A,1)
            if uplo == :U
                A[i,j] = A[j,i]
            else
                A[j,i] = A[i,j]
            end
        end
    end
    A
end

function approx_submul!(a::T,b::T,c::T, prec=precision(a)) where T <: Union{Arb,ArbRef}
    Arblib.submul!(Arblib.midref(a), Arblib.midref(b), Arblib.midref(c), prec=prec)
end

function approx_div!(a::T,b::T,c::T, prec=precision(a)) where T <: Union{Arb,ArbRef}
    Arblib.div!(Arblib.midref(a), Arblib.midref(b), Arblib.midref(c),prec= prec)
end
function approx_sqrt!(a::T,b::T, prec=precision(a)) where T <: Union{Arb,ArbRef}
    Arblib.sqrt!(Arblib.midref(a), Arblib.midref(b),prec= prec)
end
function approx_cholesky!(A::ArbRefMatrix,B::ArbRefMatrix;prec=precision(A))
    size(A) == size(B) || error("A and B must have the same size")
    Arblib.set!(A,B)
    approx_cholesky!(A,prec=prec)
end

function approx_cholesky!(A::ArbRefMatrix;prec=precision(A))
    #The cholesky of Arblib, but with zeroed error bounds during the calculations
    size(A,1) == size(A,2) || error("The matrix must be square")

    n = size(A,1)
    Arblib.get_mid!(A,A)
    for i=1:n
        for j=1:i-1
            for k=1:j-1
                approx_submul!(A[i,j], A[i,k], A[j,k], prec)
            end
            approx_div!(A[i,j], A[i,j], A[j,j], prec)
            # Arblib.get_mid!(A[i,j],A[i,j])
        end
        for k=1:i-1
            approx_submul!(A[i,i], A[i,k], A[i,k], prec)
        end
        # Arblib.get_mid!(A[i,i], A[i,i])
        if Arblib.is_positive(A[i,i]) == 0
            return 0
        end
        approx_sqrt!(A[i,i], A[i,i], prec)
        # Arblib.get_mid!(A[i,i], A[i,i])
    end

    #zero the (strict) upper triangular part:
    for i=1:n
        for j=i+1:n
            A[i,j] = 0
        end
    end
    return 1
end

function unique(x::Vector{T}) where T<:Union{ArbMatrix,ArbRefMatrix}
    unique_idx = Int[]
    for i=1:length(x)
        #We compare this one with the elements we did already
        duplicate = false
        for j=1:i-1
            if x[i] == x[j]
                # this is a duplicate
                duplicate = true
                break
            end
        end
        if !duplicate
            push!(unique_idx,i)
        end
    end
    return x[unique_idx]
end

function matmul_threaded!(C::T, A::T, B::T; n = Threads.nthreads(), prec=precision(C), opt::Int=0) where T<:Union{ArbMatrix, ArbRefMatrix}
    # matrix multiplication C = A*B, using multithreading.
    # Note that this is not the same as Arblib.mul_threaded!(), because that uses the classical matrix multiplication with flint threads
    # We will have to tune the parameters, because Arblib does classical mat mul when some dimension of the matrix has size <= 40 (i.e. A has <= 40 rows or columns, or B has <=40 rows or columns)
    # We might have very rectangular matrices (i.e., the number of rows much larger than the number of columns or vice versa), in which case it is easy do distribute over threads. When the number of remaining rows/columns can be below their cutoff due to our threading, we might want to consider using less threads? or just call the block version directly?

    # n = Threads.nthreads()
    #check dimensions
    @assert size(C,1) == size(A,1) && size(C,2) == size(B,2) && size(A,2) == size(B,1)
    #check for the special case of 1 thread:
    if n==1 || size(A,1)*size(A,2)*size(B,2) * div(prec,256) < 1000
        #TODO: determine up to what sizes (& what number of threads) running multithreaded is better than threaded.
        return Arblib.approx_mul!(C,A,B,prec=prec)
    end
    # determine how to do threading
    if !(opt in [1,2,3,4])
        opt = argmax([size(B,2),size(A,1),size(A,2)])
    end

    # opt == 1 => distribute the columns of B over threads (preferred above rows when equal due to julia/arb being column-major)
    # opt == 2 => distribute the rows of A over threads
    # opt == 3 => distribute the columns of A/rows of B over threads and add the results
    # Option 1 and 2 are preferred above 3, because in those cases we don't have to add matrices. argmax gives the first maximal element
    # for now, we'll just do one of the options. Maybe in some cases it is better to have for example a 3x3 grid instead of 1x9, so to use multiple options instead of 1
    # That will just be generating different idx_rows, idx_cols and inner stuff

    if opt == 1
        thread_func_right!(C,A,B,Arblib.approx_mul!,prec=prec)
        return C
    else
        idx_cols = [1:size(B,2) for i=1:n]
    end
    if opt == 2
        thread_func_left!(C,A,B,Arblib.approx_mul!,prec=prec)
        return C
    else
        idx_rows = [1:size(A,1) for i=1:n]
    end

    #If we are threading over the outer part, we only need one result. So n_inner gives the number of results, idx_inner gives the indices thread i uses, res_idx gives which result is used
    if opt == 3
        n_inner = n
        nr = size(A,2)
        # every thread gets at least floor(nr/n) rows, and we add the remaining parts until they are all distributed
        sizes = [div(nr,n_inner) for i=1:n_inner] .+ [div(nr,n_inner)*n_inner+i <= nr ? 1 : 0 for i=1:n_inner]
        idxs = [0, cumsum(sizes)...]
        idx_inner = [idxs[i]+1:idxs[i+1] for i=1:n_inner]
        res_idx = [i for i=1:n_inner]
    else
        n_inner = 1
        idx_inner = [1:size(A,2) for i=1:n]
        res_idx = [1 for i=1:n]
    end

    if opt == 4 #both over rows and columns, no inner. This option can only be chosen manually. It seems to be far less efficient
        # we want the blocks to be about square. So we need to find numbers a, b such that rows/a ≈ cols/b and a*b ≈ n
        if size(C,1) == size(C,2) #square, so we take squareroot
            sqrtn = isqrt(n) #integer squareroot, so sqrtn^2 < n probably. We
            if sqrtn^2 != n
                a = sqrtn
                b = sqrtn+1
            else
                a = b = sqrtn
            end
        else
            r = size(C,1)
            c = size(C,2)
            a = isqrt(div(n*r,c))
            b = div(n,a)
            #it is possible that we can increase a by one while keeping a*b <= n
            a = max(a, div(n,b))
        end
        n = a*b #new n based on the partition

        rsizes =  [div(size(C,1),a) for i=1:a] .+ [div(size(C,1),a)*a+i <= size(C,1) ? 1 : 0 for i=1:a]
        ridxs = [0, cumsum(rsizes)...]
        idx_rows_part = [ridxs[i]+1:ridxs[i+1] for i=1:a]

        csizes =  [div(size(C,2),b) for i=1:b] .+ [div(size(C,2),b)*b+i <= size(C,2) ? 1 : 0 for i=1:b]
        cidxs = [0, cumsum(csizes)...]
        idx_cols_part = [cidxs[i]+1:cidxs[i+1] for i=1:b]

        #final partition: iteration i does block C[idx_rows[i],idx_cols[i]]
        idx_rows = [idx_rows_part[div(i,b)+1] for i=0:n-1]
        idx_cols = [idx_cols_part[i%b+1] for i=0:n-1]

        n_inner = 1
        idx_inner = [1:size(A,2) for i=1:n]
        res_idx = [1 for i=1:n]
    end

    # Do the threaded matrix multiplication
    # We could use 1 matrix less if we use C with n_inner > 1 too.
    if n_inner > 1
        res_matrices = [T(size(C,1),size(C,2),prec=prec) for i=1:n_inner]
    end
    part_res = [T(length(idx_rows[i]),length(idx_cols[i]),prec=prec) for i=1:n] # when n_inner =1, this is max one copy of C
    As = [A[idx_rows[i],idx_inner[i]] for i=1:n] #in total max n copies of A, but in total 1 copy if n_inner = 1
    Bs = [B[idx_inner[i],idx_cols[i]] for i=1:n] # in total max n copies of B
    Threads.@threads for i=1:n
        # cur_res = T(length(idx_rows[i]),length(idx_cols[i]),prec=prec)
        # I am not sure whether allocating inside a threaded loop + finalizers can give memory leakage.
        # for ArbMatrices we would probably want to use the ref versions of the slice & set operations
        # Arblib.approx_mul!(part_res[i],A[idx_rows[i],idx_inner[i]],B[idx_inner[i],idx_cols[i]]) #do we want to use a different mul here?
        Arblib.approx_mul!(part_res[i],As[i],Bs[i])
        if n_inner > 1 # race conditions if we set/add to C directly, so we set the parts
            res_matrices[res_idx[i]][idx_rows[i],idx_cols[i]] = part_res[i]
        else
            C[idx_rows[i],idx_cols[i]] = part_res[i] # we actually would want to use a view/window of C[...,...] for the matrix product. But the Arblib window does not work nicely due to init/clear and normal view does not work with ArbMatrices.
        end
    end

    #add the results if needed
    if n_inner > 1
        # this way we don't have to zero C, we do that implicitely.
        # if n_inner>1, then we at least have 2 matrices to add
        Arblib.add!(C,res_matrices[1],res_matrices[2],prec=prec)
        for i=3:n_inner
            Arblib.add!(C,C,res_matrices[i],prec=prec)
        end
    end
    return C
end
function thread_func_right!(C::T,A::T,B::T,func!;n=Threads.nthreads(),prec=precision(C)) where T <: AbstractMatrix
    # distribute B over several cores and apply func!(C[part],A,B[part],prec=prec)
    nc = size(B,2)
    sizes = [div(nc,n) for i=1:n] .+ [div(nc,n)*n+i <= nc ? 1 : 0 for i=1:n]
    idxs = [0, cumsum(sizes)...]
    idx_cols = [idxs[i]+1:idxs[i+1] for i=1:n]

    parts_B = [B[:,idx_cols[i]] for i=1:n]
    part_res = [T(size(C,1),length(idx_cols[i]),prec=prec) for i=1:n] # when n_inner =1, this is max one copy of C
    Threads.@threads for i=1:n
        func!(part_res[i],A,parts_B[i],prec=prec)
        C[:,idx_cols[i]] = part_res[i]
    end
end

function thread_func_left!(C::T,A::T,B::T,func!;n=Threads.nthreads(),prec=precision(C)) where T <: AbstractMatrix
    # distribute B over several cores and apply func!(C[part],A,B[part],prec=prec)
    nr = size(A,1)
    sizes = [div(nr,n) for i=1:n] .+ [div(nr,n)*n+i <= nr ? 1 : 0 for i=1:n]
    idxs = [0, cumsum(sizes)...]
    idx_rows = [idxs[i]+1:idxs[i+1] for i=1:n]

    parts_A = [A[idx_rows[i],:] for i=1:n]
    part_res = [T(length(idx_rows[i]),size(C,2),prec=prec) for i=1:n] # when n_inner =1, this is max one copy of C
    Threads.@threads for i=1:n
        func!(part_res[i],parts_A[i],B,prec=prec)
        C[idx_rows[i],:] = part_res[i]
    end
end

function malloc_trim(k::Int)
    if Sys.islinux()
        #max because the parameter has to be at least 0
        ccall(:malloc_trim,Cvoid,(Cint,),max(k,0))
    end
end

# For BlockDiagonals it is not clear what the precision is, because different blocks have potentially a different precision
# So we make sure not to extend the Base.precision function for BlockDiagonals.
# However, then we need to define precision for ArbMatrices and BigFloats
precision(P::BlockDiagonal) = precision(P.blocks[1]) #we assume that all blocks have the same precision
precision(x) = Base.precision(x) #for other things, just use the Base version

ext: jl


text:
my_f(x, y) = 2 * x + y
derivative_of_my_f(x, y) = ForwardDiff.derivative(x -> my_f(x, y), x)


text:
using SparseArrays, Test
LarA = Lar.Arrangement
using LinearAlgebraicRepresentation
LAR = LinearAlgebraicRepresentation

@testset "Planar Arrangement - Big Random Data" begin

    @testset "Closure of 2-cells" begin
        for i = 1 : size(model2, 2, 1)
            edges  = findnz(model2.T[2][i, :])[1]
            points = abs.(model2.T[1][edges, :])
            @test sum(points) == 2 * length(edges)
        end
    end

    @testset "Well formed cells" begin
        @test sum(sum(abs.(model2.T[1]), dims = 1) .<= 1) == 0
        @test sum(sum(abs.(model2.T[1]), dims = 2) .!= 2) == 0
        @test sum(sum(abs.(model2.T[2]), dims = 1) .>= 3) == 0
        @test sum(sum(abs.(model2.T[2]), dims = 2) .<= 2) == 0
    end

    @testset "Euler" begin
        @test size(model2, 0, 2) - size(model2, 1, 1) + size(model2, 2, 1) == 1
    end

    @testset "1-Cell Compactness" begin
        #@test sum(abs.(model2.T[1])) == 2 * size(model, 1, 2)
    end

    @testset "2-Cell Compactness" begin
        #@test sum(abs.(model2.T[2])) == 2 * size(model2, 2, 2)
    end

    @testset "Biconnected Components" begin
        # for comp in bicon_comps
        #     sum(
        #         (sum(abs.(model.T[1][comp, :]), dims=1) .!= 0)
        #         + (sum(abs.(model.T[1][setdiff(1:size(model, 1, 1), comp), :]), dims=1) .!= 0)
        #         .!= 1
        #     ) == 0
        # end
    end
end


text:
# ideally this should be in the GenerativeAD pkg
function compute_ranks(rt)
	mask_nan_max = (x) -> (isnan(x) ? -Inf : x)
	rs = []
	for row in eachrow(rt)
		push!(rs, StatsBase.competerank(mask_nan_max.(Vector(row)), rev = true))
	end
	# each row represents ranks of a method
	reduce(hcat, rs)
end

# top 10 mean-max
function aggregate_stats_mean_max_top_10(df::DataFrame,
										criterion_col=:val_auc;
										min_samples=("anomaly_class" in names(df)) ? 10 : 3,
										add_col=nothing, verbose=true)
	agg_cols = vcat(_prefix_symbol.("val", BASE_METRICS), _prefix_symbol.("tst", BASE_METRICS))
	agg_cols = vcat(agg_cols, _prefix_symbol.("val", PAT_METRICS), _prefix_symbol.("tst", PAT_METRICS))
	agg_cols = vcat(agg_cols, _prefix_symbol.("val", PATN_METRICS), _prefix_symbol.("tst", PATN_METRICS))
	agg_cols = vcat(agg_cols, _prefix_symbol.("val", PAC_METRICS), _prefix_symbol.("tst", PAC_METRICS))
	agg_cols = vcat(agg_cols, Symbol.(TRAIN_EVAL_TIMES))
	agg_cols = (add_col !== nothing) ? vcat(agg_cols, add_col) : agg_cols
	top10_std_cols = _prefix_symbol.(agg_cols, "top_10_std")
	std_cols = _prefix_symbol.(agg_cols, "std")

	# agregate by seed over given hyperparameter and then choose best
	results = []
	for (dkey, dg) in pairs(groupby(df, :dataset))
		for (mkey, mg) in pairs(groupby(dg, :modelname))
			pg = groupby(mg, :phash)
			mg_suff = reduce(vcat, [g for g in pg if nrow(g) >= min_samples])
			pg_agg = combine(groupby(mg_suff, :phash),
						nrow => :psamples,
						agg_cols .=> mean .=> agg_cols,
						agg_cols .=> std,
						:parameters => unique => :parameters)

			# sort by criterion_col
			sort!(pg_agg, order(criterion_col, rev=true))

			# take best 10 models
			best_10 = first(pg_agg, 10)
			best = combine(best_10,
							agg_cols .=> mean .=> agg_cols,
							std_cols .=> mean .=> std_cols,
							agg_cols .=> std  .=> top10_std_cols)

			# add grouping keys
			best[:dataset] = dkey.dataset
			best[:modelname] = mkey.modelname

			push!(results, best)
		end
	end
	vcat(results...)
end

# top 10 max-mean
function aggregate_stats_max_mean_top_10(df::DataFrame,
										criterion_col=:val_auc; add_col=nothing, verbose=true)
	agg_cols = vcat(_prefix_symbol.("val", BASE_METRICS), _prefix_symbol.("tst", BASE_METRICS))
	agg_cols = vcat(agg_cols, _prefix_symbol.("val", PAT_METRICS), _prefix_symbol.("tst", PAT_METRICS))
	agg_cols = vcat(agg_cols, _prefix_symbol.("val", PATN_METRICS), _prefix_symbol.("tst", PATN_METRICS))
	agg_cols = vcat(agg_cols, _prefix_symbol.("val", PAC_METRICS), _prefix_symbol.("tst", PAC_METRICS))
	agg_cols = vcat(agg_cols, Symbol.(TRAIN_EVAL_TIMES))
	agg_cols = (add_col !== nothing) ? vcat(agg_cols, add_col) : agg_cols
	top10_std_cols = _prefix_symbol.(agg_cols, "top_10_std")
	std_cols = _prefix_symbol.(agg_cols, "std")

	agg_keys = ("anomaly_class" in names(df)) ? [:seed, :anomaly_class] : [:seed]

	# choose best for each seed/anomaly_class cobination and average over them
	results = []
	for (dkey, dg) in pairs(groupby(df, :dataset))
		for (mkey, mg) in pairs(groupby(dg, :modelname))
			partial_results = []

			if ("anomaly_class" in names(df))
				classes = unique(mg.anomaly_class)
				dif = setdiff(collect(1:10), classes)
				if (length(classes) < 10) && verbose
					@warn "$(mkey.modelname) - $(dkey.dataset): missing runs on anomaly_class $(dif)."
				end
			else
				seeds = unique(mg.seed)
				dif = setdiff(collect(1:5), seeds)
				if (length(seeds) < 3) && verbose
					@warn "$(mkey.modelname) - $(dkey.dataset): missing runs on seed $(dif)."
				end
			end

			# iterate over seed-anomaly_class groups
			for (skey, sg) in pairs(groupby(mg, agg_keys))
				ssg = sort(sg, order(criterion_col, rev=true))
				# best hyperparameter after sorting by criterion_col
				best_10 = first(ssg, 10)
				best = combine(best_10,
							agg_cols .=> mean .=> agg_cols,
							agg_cols .=> std  .=> top10_std_cols)

				push!(partial_results, best)
			end

			best_per_seed = reduce(vcat, partial_results)
			# average over seed-anomaly_class groups
			best = combine(best_per_seed,
						agg_cols .=> mean .=> agg_cols,
						top10_std_cols .=> mean .=> top10_std_cols,
						agg_cols .=> std)

			# add grouping keys
			best[:dataset] = dkey.dataset
			best[:modelname] = mkey.modelname

			push!(results, best)
		end
	end
	vcat(results...)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiBatchV1beta1CronJob <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    spec::Any # spec type: Union{ Nothing, IoK8sApiBatchV1beta1CronJobSpec } # spec name: spec
    status::Any # spec type: Union{ Nothing, IoK8sApiBatchV1beta1CronJobStatus } # spec name: status

    function IoK8sApiBatchV1beta1CronJob(;apiVersion=nothing, kind=nothing, metadata=nothing, spec=nothing, status=nothing)
        o = new()
        validate_property(IoK8sApiBatchV1beta1CronJob, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiBatchV1beta1CronJob, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiBatchV1beta1CronJob, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiBatchV1beta1CronJob, Symbol("spec"), spec)
        setfield!(o, Symbol("spec"), spec)
        validate_property(IoK8sApiBatchV1beta1CronJob, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        o
    end
end # type IoK8sApiBatchV1beta1CronJob

const _property_map_IoK8sApiBatchV1beta1CronJob = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("spec")=>Symbol("spec"), Symbol("status")=>Symbol("status"))
const _property_types_IoK8sApiBatchV1beta1CronJob = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("spec")=>"IoK8sApiBatchV1beta1CronJobSpec", Symbol("status")=>"IoK8sApiBatchV1beta1CronJobStatus")
Base.propertynames(::Type{ IoK8sApiBatchV1beta1CronJob }) = collect(keys(_property_map_IoK8sApiBatchV1beta1CronJob))
Swagger.property_type(::Type{ IoK8sApiBatchV1beta1CronJob }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiBatchV1beta1CronJob[name]))}
Swagger.field_name(::Type{ IoK8sApiBatchV1beta1CronJob }, property_name::Symbol) =  _property_map_IoK8sApiBatchV1beta1CronJob[property_name]

function check_required(o::IoK8sApiBatchV1beta1CronJob)
    true
end

function validate_property(::Type{ IoK8sApiBatchV1beta1CronJob }, name::Symbol, val)
end


text:
# fallback methods

function project(a::Sequence, space_dest::VectorSpace, ::Type{T}=eltype(a)) where {T}
    _iscompatible(space(a), space_dest) || return throw(DimensionMismatch)
    c = Sequence(space_dest, Vector{T}(undef, dimension(space_dest)))
    _project!(c, a)
    return c
end

function project!(c::Sequence, a::Sequence)
    _iscompatible(space(a), space(c)) || return throw(DimensionMismatch)
    _project!(c, a)
    return c
end

function _project!(c::Sequence, a::Sequence)
    space_a = space(a)
    space_c = space(c)
    if space_a == space_c
        coefficients(c) .= coefficients(a)
    elseif space_c ⊆ space_a
        @inbounds for α ∈ indices(space_c)
            c[α] = a[α]
        end
    else
        coefficients(c) .= zero(eltype(c))
        @inbounds for α ∈ indices(space_a ∩ space_c)
            c[α] = a[α]
        end
    end
    return c
end

#

function project(A::LinearOperator, domain_dest::VectorSpace, codomain_dest::VectorSpace, ::Type{T}=eltype(A)) where {T}
    _iscompatible(domain(A), domain_dest) & _iscompatible(codomain(A), codomain_dest) || return throw(DimensionMismatch)
    C = LinearOperator(domain_dest, codomain_dest, Matrix{T}(undef, dimension(codomain_dest), dimension(domain_dest)))
    _project!(C, A)
    return C
end

function project!(C::LinearOperator, A::LinearOperator)
    _iscompatible(domain(A), domain(C)) & _iscompatible(codomain(A), codomain(C)) || return throw(DimensionMismatch)
    _project!(C, A)
    return C
end

function _project!(C::LinearOperator, A::LinearOperator)
    domain_A, codomain_A = domain(A), codomain(A)
    domain_C, codomain_C = domain(C), codomain(C)
    if domain_A == domain_C && codomain_A == codomain_C
        coefficients(C) .= coefficients(A)
    elseif domain_C ⊆ domain_A && codomain_C ⊆ codomain_A
        @inbounds for β ∈ indices(domain_C), α ∈ indices(codomain_C)
            C[α,β] = A[α,β]
        end
    else
        coefficients(C) .= zero(eltype(A))
        @inbounds for β ∈ indices(domain_A ∩ domain_C), α ∈ indices(codomain_A ∩ codomain_C)
            C[α,β] = A[α,β]
        end
    end
    return C
end

# Cartesian spaces

function _project!(c::Sequence{<:CartesianSpace}, a::Sequence{<:CartesianSpace})
    space_c = space(c)
    if space(a) == space_c
        coefficients(c) .= coefficients(a)
    else
        @inbounds for i ∈ 1:nb_cartesian_product(space_c)
            _project!(component(c, i), component(a, i))
        end
    end
    return c
end

#

function _project!(C::LinearOperator{<:CartesianSpace,<:CartesianSpace}, A::LinearOperator{<:CartesianSpace,<:CartesianSpace})
    domain_C = domain(C)
    codomain_C = codomain(C)
    if domain(A) == domain_C && codomain(A) == codomain_C
        coefficients(C) .= coefficients(A)
    else
        @inbounds for j ∈ 1:nb_cartesian_product(domain_C), i ∈ 1:nb_cartesian_product(codomain_C)
            _project!(component(C, i, j), component(A, i, j))
        end
    end
    return C
end

function _project!(C::LinearOperator{<:CartesianSpace,<:VectorSpace}, A::LinearOperator{<:CartesianSpace,<:VectorSpace})
    domain_C = domain(C)
    if domain(A) == domain_C && codomain(A) == codomain(C)
        coefficients(C) .= coefficients(A)
    else
        @inbounds for j ∈ 1:nb_cartesian_product(domain_C)
            _project!(component(C, j), component(A, j))
        end
    end
    return C
end

function _project!(C::LinearOperator{<:VectorSpace,<:CartesianSpace}, A::LinearOperator{<:VectorSpace,<:CartesianSpace})
    codomain_C = codomain(C)
    if domain(A) == domain(C) && codomain(A) == codomain_C
        coefficients(C) .= coefficients(A)
    else
        @inbounds for i ∈ 1:nb_cartesian_product(codomain_C)
            _project!(component(C, i), component(A, i))
        end
    end
    return C
end


text:
using DifferentialEquations
using Plots; gr()
H = <HAMILTONIAN> :: Array{Complex{Float64},2}
f(u,p,t) = im*(u*H - H*u)
u0 = <START> :: Array{Complex{Float64},2}
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob)
t = sol.t
u = sol.u
#pl1 = plot(sol)
#gui()

max_stars_repo_name: vleschuk/odesolve

max_issues_repo_name: vleschuk/odesolve

max_forks_repo_name: vleschuk/odesolve


text:
#------------------------------------------------------------------------------
function SolveModel(A, obj, b, n)
    # Step 1: Create the model]
        m1 = Model(solver=GurobiSolver(OutputFlag = 0))
        @variable(m1, x[1:n])
        @objective(m1, Max, dot(obj, x))
        for i = 1:length(A[:, 1])
            @constraint(m1, dot(A[i, :], x) <= b[i])
        end

    # Step 2: Solve the model
        status = solve(m1; suppress_warnings = true)
    # Step 3: output optimal value;
        val = getvalue(x);
    end
    #------------------------------------------------------------------------------

    function SolveFP(cone, iteration)
        epsilon = 1
        # M2 = 10
        m5 = Model(solver=GurobiSolver(OutputFlag = 0))
        @variable(m5, y[1:n]<=-1e-4)
        @variable(m5, gamma[1:iteration, 1:c] >= 0)
        @variable(m5, w, Bin)

        for i = 1:iteration
            for j = 1:n
                @constraint(m5, y[j] == sum(gamma[i,k]*cone[i][k,j] for k = 1:length(cone[i][:,1])))
            end
        end
        # @constraint(m5, epsilon <= M2*w + sum(y))
        # @constraint(m5, epsilon <= M2*(1-w) - sum(y))
        @constraint(m5, sum(y) == -epsilon)

        @objective(m5, Min, 0)
        status = solve(m5; suppress_warnings = true)
        return status
    end

    #------------------------------------------------------------------------------

    function SampleC(cone)
        iteration = length(cone[:])
        m3 = Model(solver=GurobiSolver(OutputFlag = 0, Seed = rand(rng, 1:1000)))
        @variable(m3, y[1:n] <= 0)
        @variable(m3, 1 >= gamma[1:iteration, 1:c] >= 0)
        T = rand(rng, iteration, c)
        for i = 1:iteration
            for j = 1:n
                @constraint(m3, y[j] == sum(gamma[i,k]*cone[i][k,j] for k = 1:length(cone[i][:,1])))
            end
        end
        @objective(m3, Max, sum(T[i,k]*gamma[i,k] for i = 1:iteration for k = 1:length(cone[i][:,1])))
        status = solve(m3; suppress_warnings = true)
        return getvalue(y)
    end

    #------------------------------------------------------------------------------
    function tightind(A, b, x)
        # Outputs an array with the indices of active constraints
        epsilon = 5
        istight = []
        z = zeros(Int16, c)
        for i = 1:c
            if (abs(-dot(A[i,:], x) + b[i]) <= epsilon)
                push!(istight, i)
                z[i] = 1
            end
        end
        return istight, z
    end

    #------------------------------------------------------------------------------
    function SolveP2(Cprime, cone, iteration)
        m3 = Model(solver=GurobiSolver(OutputFlag = 0))
        @variable(m3, y[1:n])
        @variable(m3, gamma[1:iteration, 1:c] >= 0)

        for i = 1:iteration
            for j = 1:n
                @constraint(m3, y[j] == sum(gamma[i,k]*cone[i][k,j] for k = 1:length(cone[i][:,1])))
            end
        end

        @objective(m3, Min, sum((y[j]-Cprime[j])^2 for j = 1:n))
        status = solve(m3; suppress_warnings = true)
        return getvalue(y), getobjectivevalue(m3)
    end

    #------------------------------------------------------------------------------

    function SolveP1(A, b, x, NoS, J)
        n = length(A[1, :, 1])
        m2 = Model(solver = GurobiSolver(OutputFlag = 0, IntFeasTol = 1e-9, MIPGapAbs = 1e-5, MIPGap = 0))
        M = 2000
        epsilon = 1

        @variable(m2, xbar[1:NoS, 1:n])
        @variable(m2, lambda[1:NoS, 1:c] >= 0)
        @variable(m2, s[1:NoS, 1:c] >= 0)
        @variable(m2, z[1:NoS, 1:c], Bin)
        @variable(m2, C[1:n]<=-1e-3)
        @variable(m2, t[1:J, 1:n, 1:NoS] >= 0)

        @constraint(m2, DF[k = 1:NoS, j = 1:n], C[j] - sum(lambda[k, i]*A[i, j, k] for i = 1:c) == 0)
        @constraint(m2, PF[k = 1:NoS, i = 1:c], sum(A[i, j, k]*xbar[k, j] for j = 1:n) + s[k, i] == b[i])
        @constraint(m2, CS1[k = 1:NoS, i = 1:c], s[k, i] <= M*(1-z[k, i]))
        @constraint(m2, CS2[k = 1:NoS, i = 1:c], lambda[k, i] <= M*z[k, i])
        @constraint(m2, VC[k = 1:NoS] , sum(z[k,:]) >= n)
        @constraint(m2, linearization1[j = 1:J, i = 1:n, k = 1:NoS], t[j,i,k] >= xbar[k,i] - x[j,i,k])
        @constraint(m2, linearization2[j = 1:J, i = 1:n, k = 1:NoS], t[j,i,k] >= -(xbar[k,i] - x[j,i,k]))
        @constraint(m2, sum(C) == -epsilon)
        @objective(m2, Min, sum(t))
        status = solve(m2; suppress_warnings = true)
        loss1 = getobjectivevalue(m2)
        return getvalue(xbar), loss1, getvalue(z)
    end

    #-------------------------------------------------------------------------------
    function SolveP1m(A, b, x, NoS, J, trial, initx, initz)
        n = length(A[1, :, 1])
        initz = round.(initz)
        m2 = Model(solver = GurobiSolver(OutputFlag = 1, TimeLimit = 7200, ConcurrentMIP = 3, IntFeasTol = 1e-9, MIPGapAbs = 1e-5, MIPGap = 0))
        M = 2000
        epsilon = 1

        @variable(m2, xbar[j = 1:NoS, k = 1:n])
        @variable(m2, lambda[1:NoS, 1:c] >= 0)
        @variable(m2, s[1:NoS, 1:c] >= 0)
        @variable(m2, z[j = 1:NoS, k = 1:c], Bin)
        for j = 1:NoS-1
            for k = 1:c
                setvalue(z[j,k], initz[j,k])
            end
        end
        @variable(m2, C[1:n]<=-1e-3)
        @variable(m2, t[1:J, 1:n, 1:NoS] >= 0)

        @constraint(m2, DF[k = 1:NoS, j = 1:n], C[j] - sum(lambda[k, i]*A[i, j, k] for i = 1:c) == 0)
        @constraint(m2, PF[k = 1:NoS, i = 1:c], sum(A[i, j, k]*xbar[k, j] for j = 1:n) + s[k, i] == b[k, i])
        @constraint(m2, CS1[k = 1:NoS, i = 1:c], s[k, i] <= M*(1-z[k, i]))
        @constraint(m2, CS2[k = 1:NoS, i = 1:c], lambda[k, i] <= M*z[k, i])
        @constraint(m2, VC[k = 1:NoS] , sum(z[k,:]) >= n)
        @constraint(m2, -epsilon == sum(C))
        @constraint(m2, linearization1[j = 1:J, i = 1:n, k = 1:NoS], t[j,i,k] >= xbar[k,i] - x[j,i,k])
        @constraint(m2, linearization2[j = 1:J, i = 1:n, k = 1:NoS], t[j,i,k] >= -(xbar[k,i] - x[j,i,k]))
        # JuMP.build(m2)
        # grb = JuMP.internalmodel(m2).inner
        # Gurobi.write_model(grb, "test.mps")

        @objective(m2, Min, sum(t))
        status = solve(m2; suppress_warnings = true)
        loss1 = getobjectivevalue(m2)
        return getvalue(xbar), loss1, getvalue(z), status
    end

    #-------------------------------------------------------------------------------
    function infnorm(xtest, xhat)
        NoS = length(xtest[:, 1])
        n = length(xtest[1,:])
        error = zeros(Float64, NoS)
        for i = 1:NoS
            error[i] = norm(xtest[i,:] - xhat[i,:], Inf)
        end
        return sum(error)/NoS
    end
    #-------------------------------------------------------------------------------

max_stars_repo_name: ddolab/InvLP

max_issues_repo_name: ddolab/InvLP

max_forks_repo_name: ddolab/InvLP


text:
###### Operations for JuMPExpressionMatrix ######

function remove_undef!(expression_array::T) where {T <: JuMPExpressionMatrix}
    for j in 1:size(expression_array)[2]
        for i in 1:size(expression_array)[1]
            if !isassigned(expression_array, i, j)
                expression_array[i, j] = zero(eltype(expression_array))
            else
                continue
            end
        end
    end

    return
end

###### Operations for JuMPDenseAxisExpressionMatrix ######

function remove_undef!(axis_expression_array::JuMP.Containers.DenseAxisArray)
    remove_undef!(axis_expression_array.data)
    return
end

function add_to_expression!(
    expression_array::T,
    ix::Int,
    jx::Int,
    var::JV,
    multiplier::Float64,
) where {T, JV <: JuMP.AbstractVariableRef}
    if isassigned(expression_array, ix, jx)
        JuMP.add_to_expression!(expression_array[ix, jx], multiplier, var)
    else
        expression_array[ix, jx] = multiplier * var
    end

    return
end

function add_to_expression!(
    expression_array::T,
    ix::Int,
    jx::Int,
    var::JV,
    multiplier::Float64,
    constant::Float64,
) where {T, JV <: JuMP.AbstractVariableRef}
    if isassigned(expression_array, ix, jx)
        JuMP.add_to_expression!(expression_array[ix, jx], multiplier, var)
        JuMP.add_to_expression!(expression_array[ix, jx], constant)
    else
        expression_array[ix, jx] = multiplier * var + constant
    end

    return
end

function add_to_expression!(expression_array::T, ix::Int, jx::Int, value::Float64) where {T}
    if isassigned(expression_array, ix, jx)
        expression_array[ix, jx].constant += value
    else
        expression_array[ix, jx] = zero(eltype(expression_array)) + value
    end

    return
end

function add_to_expression!(
    expression_array::T,
    ix::Int,
    jx::Int,
    parameter::PJ.ParameterRef,
) where {T}
    if isassigned(expression_array, ix, jx)
        JuMP.add_to_expression!(expression_array[ix, jx], 1.0, parameter)
    else
        expression_array[ix, jx] = zero(eltype(expression_array)) + parameter
    end

    return
end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using UnitCommitment, JSON, GZip, DataStructures

basedir = @__DIR__

function parse_case14()
    return JSON.parse(
        GZip.gzopen("$basedir/../../instances/test/case14.json.gz"),
        dicttype = () -> DefaultOrderedDict(nothing),
    )
end

@testset "repair!" begin
    @testset "Cost curve should be convex" begin
        json = parse_case14()
        json["Generators"]["g1"]["Production cost curve (MW)"] = [100, 150, 200]
        json["Generators"]["g1"]["Production cost curve (\$)"] = [10, 25, 30]
        instance = UnitCommitment._from_json(json, repair = false)
        @test UnitCommitment.repair!(instance) == 4
    end

    @testset "Startup limit must be greater than Pmin" begin
        json = parse_case14()
        json["Generators"]["g1"]["Production cost curve (MW)"] = [100, 150]
        json["Generators"]["g1"]["Production cost curve (\$)"] = [100, 150]
        json["Generators"]["g1"]["Startup limit (MW)"] = 80
        instance = UnitCommitment._from_json(json, repair = false)
        @test UnitCommitment.repair!(instance) == 1
    end

    @testset "Startup costs and delays must be increasing" begin
        json = parse_case14()
        json["Generators"]["g1"]["Startup costs (\$)"] = [300, 200, 100]
        json["Generators"]["g1"]["Startup delays (h)"] = [8, 4, 2]
        instance = UnitCommitment._from_json(json, repair = false)
        @test UnitCommitment.repair!(instance) == 4
    end
end


text:
"""
  smm_slices(sMMProblem::SMMProblem)

Function to plot slices of the objective function.
"""
function smm_slices(sMMProblem::SMMProblem, paramValues::Vector, nbPoints::Int64; showPlots::Bool = true)

    listPlots = []
    p = Plots.plot()

    # Loop over parameter values
    #---------------------------
    for (keyIndex, keyValue) in enumerate(keys(sMMProblem.priors))

        vXGrid = SharedArray(collect(linspace(sMMProblem.priors[keyValue][2], sMMProblem.priors[keyValue][3], nbPoints)))
        vYGrid = SharedArray(zeros(nbPoints))

        info("slicing along $(keyValue)")
        # Options A.
        #-----------
        @sync @parallel for xIndex = 1:1:length(vXGrid)

          # Move along one dimension, keeping other values constant
          #--------------------------------------------------------
          localParamValues = copy(paramValues)
          localParamValues[keyIndex] = vXGrid[xIndex]

          vYGrid[xIndex] = sMMProblem.objective_function(localParamValues)

        end

        p = Plots.plot(vXGrid, vYGrid, title = "$(keyValue)")
        push!(listPlots, p)

        if showPlots == true
            display(p)
        end

    end

    return listPlots

end


text:
"""
    WeightedSum(T = Float64)

Simple weighted sum, tracked as type `T`.

# Example:
    o = fit!(WeightedSum(), rand(100), rand(100))
    sum(o)
"""
mutable struct WeightedSum{T} <: WeightedOnlineStat{T}
    ∑::T
    W::T
    n::Int
    function WeightedSum{T}(∑ = T(0), W = T(0), n = 0) where T
        new{T}(T(∑), T(W), Int(n))
    end
end

WeightedSum(∑::T, W::T, n::Int) where T = WeightedSum{T}(∑, W, n)
WeightedSum(::Type{T}) where T = WeightedSum(T(0), T(0), 0)
WeightedSum() = WeightedSum(Float64)

function WeightedOnlineStats._fit!(o::WeightedSum{T}, x, w) where T
    ww = convert(T, w)
    xx = convert(T, x)

    o.n += 1
    o.W += smooth(o.W, ww, T(1) / o.n)
    o.∑ += xx * ww
    o
end

function WeightedOnlineStats._merge!(o::WeightedSum{T}, o2::WeightedSum) where T
    o.n += o2.n
    o.W = smooth(o.W, convert(T, o2.W), convert(T, o2.n / o.n))
    o.∑ += convert(T, o2.∑)
    o
end

OnlineStatsBase.value(o::WeightedSum) = o.∑
Base.sum(o::WeightedSum) = value(o)
Base.copy(o::WeightedSum) = WeightedSum(o.∑, o.W, o.n)


text:
module ExtendableSparse
using DocStringExtensions
using SparseArrays
using LinearAlgebra

include("sparsematrixcsc.jl")
include("sparsematrixlnk.jl")
include("extendable.jl")
include("preconditioners.jl")
include("sprand.jl")

export SparseMatrixLNK,ExtendableSparseMatrix,flush!,nnz, sprand!,sprand_sdd!, fdrand,fdrand!
export JacobiPreconditioner, ILU0Preconditioner, updateindex!

export colptrs
end # module


text:
"""
Covariance Matrix Adaptation Evolution Strategy Implementation: (μ/μ_{I,W},λ)-CMA-ES

The constructor takes following keyword arguments:

- `μ`/`mu` is the number of parents
- `λ`/`lambda` is the number of offspring
- `c_1` is a learning rate for the rank-one update of the covariance matrix update
- `c_c` is a learning rate for cumulation for the rank-one update of the covariance matrix
- `c_mu` is a learning rate for the rank-``\\mu`` update of the covariance matrix update
- `c_sigma` is a learning rate for the cumulation for the step-size control
- `c_m` is the learning rate for the mean update, ``c_m \\leq 1``
- `σ0`/`sigma0` is the initial step size `σ`
- `weights` are recombination weights, if the weights are set to ``1/\\mu`` then the *intermediate* recombination is activated.
- `metrics` is a collection of convergence metrics.
"""
struct CMAES{T} <: AbstractOptimizer
    μ::Int
    λ::Int
    c_1::T
    c_c::T
    c_μ::T
    c_σ::T
    σ₀::T
    cₘ::T
    wᵢ::Vector{T}
    metrics::ConvergenceMetrics

    function CMAES(; μ::Int=10, λ::Int=2μ, mu::Int=μ, lambda::Int=0, weights::Vector{T}=zeros(lambda),
                     c_1::Real=NaN, c_c::Real=NaN, c_mu::Real=NaN, c_sigma::Real=NaN,
                     sigma0::Real=0.5, c_m::Real=1, metrics=ConvergenceMetric[]) where {T}
        @assert c_m ≤ 1 "cₘ > 1"
        if lambda == 0
            lambda = μ == mu ? λ : 2*mu
            weights = zeros(lambda)
        end
        if mu != lambda/2
            mu = lambda >> 1
        end
        @assert length(weights) == lambda "Number of weights must be $lambda"
        if length(metrics) == 0
            push!(metrics, AbsDiff(T(1e-12)))
        end
        new{T}(mu, lambda, c_1, c_c, c_mu, c_sigma, sigma0, c_m, weights, metrics)
    end
end
summary(m::CMAES) = "($(m.μ),$(m.λ))-CMA-ES"
show(io::IO,m::CMAES) = print(io, summary(m))
population_size(method::CMAES) = method.μ
default_options(method::CMAES) = (iterations=1500,)

mutable struct CMAESState{T,TI} <: AbstractOptimizerState
    N::Int
    μ_eff::T
    c_1::T
    c_c::T
    c_μ::T
    c_σ::T
    fitpop::Vector{T}
    C::Matrix{T}
    s::Vector{T}
    s_σ::Vector{T}
    σ::T
    wᵢ::Vector{T}
    d_σ::T
    parent::TI
    fittest::TI
    z::Matrix{T}
end
value(s::CMAESState) = first(s.fitpop)
minimizer(s::CMAESState) = s.fittest

"""Initialization of CMA-ES algorithm state"""
function initial_state(method::CMAES, options, objfun, population)
    @unpack μ,λ,c_1,c_c,c_μ,c_σ,σ₀,cₘ,wᵢ = method
    @assert μ < λ "Offspring population must be larger then parent population"

    T = typeof(value(objfun))
    individual = first(population)
    TI = typeof(individual)
    n = length(individual)
    α_cov = 2

    # set parameters
    μ_eff = 1/sum(w->w^2,wᵢ[1:μ])
    if !(1 ≤ μ_eff ≤ μ)
        # default parameters for weighted recombination
        w′ᵢ = T[log((λ+1)/2.0)-log(i) for i in 1:λ]
        μ = sum(w′ᵢ .> 0)
        μ_eff = sum(w′ᵢ[1:μ])^2/sum(w->w^2,w′ᵢ[1:μ])
        μ⁻_eff = sum(w′ᵢ[μ+1:end])^2/sum(w->w^2,w′ᵢ[μ+1:end])
        c_1 = isnan(c_1) ? α_cov/((n+1.3)^2+μ_eff) : c_1
        c_μ = isnan(c_μ) ? min(1 - c_1, α_cov*(μ_eff-2+1/μ_eff)/((n+2)^2+α_cov*μ_eff/2)) : c_μ
        c_c = isnan(c_c) ? (4 + μ_eff/n)/(n + 4 + 2*μ_eff/n) : c_c
        c_σ = isnan(c_σ) ? (μ_eff+2)/(n+μ_eff+5) : c_σ
        w⁺ = sum(w′ᵢ[w′ᵢ .≥ 0])
        wᵢ = [w/w⁺ for w in w′ᵢ]
        w⁻ = sum(wᵢ[wᵢ .< 0])
        wᵢ[μ+1:end] ./= -w⁻
        if wᵢ[end] < 0
            if c_μ > 0
                setweights!(wᵢ, μ, 1+c_1/c_μ)
                boundweights!(wᵢ, μ, (1-c_1-c_μ)/c_μ/n )
            end
            boundweights!(wᵢ, μ, 1+2*μ⁻_eff/(μ_eff+2))
        end
    else
        c_c = isnan(c_c) ? 1/sqrt(n) : c_c
        c_σ = isnan(c_σ) ? 1/sqrt(n) : c_σ
        c_μ = isnan(c_μ) ? min(μ_eff/n^2, 1-c_1) : c_μ
        c_1 = isnan(c_1) ? min(2/n^2, 1-c_μ) : c_1
    end

    @assert 1 ≤ μ_eff ≤ μ "μ_eff ∉ [1, μ]"
    @assert c_1 + c_μ ≤ 1 "c_1 + c_μ > 1"
    @assert c_σ < 1 "c_σ ≥ 1"
    @assert c_c ≤ 1 "c_c > 1"
    @assert 0 < wᵢ[1] ≤ 1 "0 < wᵢ[1] ≤ 1"
    @assert wᵢ[end] ≤ 0 "wᵢ[end] ≤ 0"
    @assert all(wᵢ[i] >= wᵢ[i+1] for i in 1:length(wᵢ)-1) "wᵢ are not monotonous"
    @assert wᵢ[μ] > 0 >= wᵢ[μ+1] "There are $μ positive weights"

    d_σ = 1 + 2 * max(0, sqrt((μ_eff-1)/(n+1))-1) + c_σ

    # setup initial state
    return CMAESState{T,TI}(n, μ_eff, c_1, c_c, c_μ, c_σ,
            fill(convert(T, Inf), μ),
            diagm(0=>ones(T,n)),
            zeros(T, n), zeros(T, n), method.σ₀, wᵢ, d_σ,
            copy(individual), copy(individual), zeros(T, n, λ) )
end

function update_state!(objfun, constraints, state::CMAESState{T,TI}, population::AbstractVector{TI},
                       method::CMAES, options, itr) where {T, TI}
    μ, λ, cₘ = method.μ, method.λ, method.cₘ
    N,μ_eff,σ,w,d_σ = state.N, state.μ_eff, state.σ, state.wᵢ, state.d_σ
    c_1,c_c,c_μ,c_σ = state.c_1, state.c_c, state.c_μ, state.c_σ
    weights = view(w, 1:μ)
    evaltype = options.parallelization

    parent = reshape(state.parent,N)
    parentshape = size(state.parent)

    randn!(options.rng, state.z)
    # y = zeros(T, N, λ)
    ȳ = zeros(T, N)
    offspring = Array{Vector{T}}(undef, λ)
    fitoff = fill(Inf, λ)

    B, D = try
        F = eigen!(Symmetric(state.C))
        F.vectors, Diagonal(sqrt.(max.(0,F.values)))
    catch ex
        @error "Break on eigendecomposition: $ex: $(state.C)"
        return true
    end

    C = σ * B * D
    for i in 1:λ
        # offspring are generated by transforming standard normally distributed random vectors using a transformation matrix
        offspring[i] = apply!(constraints, parent + C * state.z[:,i])
    end
    # calculate fitness of the population
    value!(objfun, fitoff, reshape.(offspring,parentshape...))
    # apply penalty to fitness
    penalty!(fitoff, constraints, offspring)

    # Select new parent population
    idx = sortperm(fitoff)
    for i in 1:μ
        o = offspring[idx[i]]
        population[i] = reshape(o,parentshape...)
        ȳ .+=(o.-parent).*(w[i]/σ)
        state.fitpop[i] = fitoff[idx[i]]
    end

    # recombination
    z_λ = view(state.z,:,idx)
    z̄ = vec(sum((@view z_λ[:,1:μ]).*weights', dims=2))
    parent += (cₘ*σ).*ȳ  #  forming recombinant parent for next generation
    # update evolution paths
    state.s_σ .= (1-c_σ).*state.s_σ .+ sqrt(μ_eff*c_σ*(2 - c_σ))*(B*z̄)
    h_σ = norm(state.s_σ)/N/(1-(1-c_σ)^(2*itr/λ)) < (2 + 4/(N+1))
    state.s = (1 - c_c)*state.s + (h_σ*sqrt(μ_eff*c_c*(2 - c_c)))*ȳ
    # perform rank-one update
    rank_1 = c_1.*state.s*state.s'
    # perform rank-μ update
    rank_μ = c_μ*sum( (w ≥ 0 ? 1 : N/norm(B*zi)^2)*w*(zi*zi') for (w,zi) in zip(w, eachcol(z_λ)) )
    # update covariance
    c1a = c_1 * (1 - (1-h_σ)*c_c*(2-c_c))
    state.C .= (1 - c1a - c_μ*sum(w)).*state.C + rank_1 + rank_μ
    # adapt step-size σ
    state.σ = σ * exp(min(1, (c_σ/d_σ)*(norm(state.s_σ)/N-1)/2))

    state.fittest = population[1]
    state.parent = reshape(parent,parentshape...)

    return false
end

function setweights!(w, μ, val)
    @assert w[μ+1] < 0 "There are no negative weights: $(w[μ-1:μ+1])"
    λ = length(w)
    if w[end] >= 0
        w[μ+1:λ] .= -abs(val)/(λ-μ)
    end
    fct = abs(val / sum(w[μ+1:end]))
    w[μ+1:end] .*= fct
    return w
end

function boundweights!(w, μ, val)
    sum(w[μ+1:end]) >= -abs(val) && return w
    @assert w[end] < 0 && w[μ+1] <= 0 "Negative weights are not set"
    fct = abs(val / sum(w[μ+1:end]))
    if fct < 1
        w[μ+1:end] .*= fct
    end
    return w
end


text:
module HydraulicFracturing

include("extra_file.jl")

end


text:
using Radio, Winston

# generate random 3 bit data modulate
data = rand( 0:7, 10000 )

# generate 10,000 random QPSK symbols
symbols = pskmod( data, 8 )
# create some gaussian noise and add it to the symbols
noise  = wgn( length( symbols ), 10, "dBm", 1.0, true )
signal = symbols .+ noise

constellation = plot_constellation( signal )
setattr( constellation, title = "8-PSK Modulation" )

display( constellation )


text:
@doc raw"""
    patricle_swarm(M, F)

perform the particle swarm optimization algorithm (PSO), starting with the initial particle positions ``x_0``[^Borckmans2010].
The aim of PSO is to find the particle position ``g`` on the `Manifold M` that solves
```math
\min_{x ∈\mathcal{M}} F(x).
```
To this end, a swarm of particles is moved around the `Manifold M` in the following manner.
For every particle ``k`` we compute the new particle velocities ``v_k^{(i)}`` in every step ``i`` of the algorithm by
```math
v_k^{(i)} = ω \, \operatorname{T}_{x_k^{(i)}\gets x_k^{(i-1)}}v_k^{(i-1)} + c \,  r_1  \operatorname{retr}_{x_k^{(i)}}^{-1}(p_k^{(i)}) + s \,  r_2 \operatorname{retr}_{x_k^{(i)}}^{-1}(g),
```
where ``x_k^{(i)}`` is the current particle position, ``ω`` denotes the inertia,
``c`` and ``s`` are a cognitive and a social weight, respectively,
``r_j``, ``j=1,2`` are random factors which are computed new for each particle and step,
``\operatorname{retr}^{-1}`` denotes an inverse retraction on the `Manifold` `M`, and
``\operatorname{T}`` is a vector transport.

Then the position of the particle is updated as
```math
x_k^{(i+1)} = \operatorname{retr}_{x_k^{(i)}}(v_k^{(i)}),
```
where ``\operatorname{retr}`` denotes a retraction on the `Manifold` `M`. At the end of each step for every particle, we set
```math
p_k^{(i+1)} = \begin{cases}
x_k^{(i+1)},  & \text{if } F(x_k^{(i+1)})<F(p_{k}^{(i)}),\\
p_{k}^{(i)}, & \text{else,}
\end{cases}
```
and
```math
g_k^{(i+1)} =\begin{cases}
p_k^{(i+1)},  & \text{if } F(p_k^{(i+1)})<F(g_{k}^{(i)}),\\
g_{k}^{(i)}, & \text{else,}
\end{cases}
```
i.e. ``p_k^{(i)}`` is the best known position for the particle ``k`` and ``g^{(i)}`` is the global best known position ever visited up to step ``i``.

[^Borckmans2010]:
    > P. B. Borckmans, M. Ishteva, P.-A. Absil, __A Modified Particle Swarm Optimization Algorithm for the Best Low Multilinear Rank Approximation of Higher-Order Tensors__,
    > In: Dorigo M. et al. (eds) Swarm Intelligence. ANTS 2010. Lecture Notes in Computer Science, vol 6234. Springer, Berlin, Heidelberg,
    > doi [10.1007/978-3-642-15461-4_2](https://doi.org/10.1007/978-3-642-15461-4_2)

# Input
* `M` – a manifold ``\mathcal M``
* `F` – a cost function ``F:\mathcal M→ℝ`` to minimize

# Optional
* `n` - (`100`) number of random initial positions of x0
* `x0` – the initial positions of each particle in the swarm ``x_k^{(0)} ∈ \mathcal M`` for ``k = 1, \dots, n``, per default these are n [`random_point`](@ref)s
* `velocity` – a set of tangent vectors (of type `AbstractVector{T}`) representing the velocities of the particles, per default a [`random_tangent`](@ref) per inital position
* `inertia` – (`0.65`) the inertia of the patricles
* `social_weight` – (`1.4`) a social weight factor
* `cognitive_weight` – (`1.4`) a cognitive weight factor
* `retraction_method` – (`ExponentialRetraction()`) a `retraction(M,x,ξ)` to use.
* `inverse_retraction_method` - (`LogarithmicInverseRetraction()`) an `inverse_retraction(M,x,y)` to use.
* `vector_transport_mthod` - (`ParallelTransport()`) a vector transport method to use.
* `stopping_criterion` – ([`StopWhenAny`](@ref)`(`[`StopAfterIteration`](@ref)`(500)`, [`StopWhenChangeLess`](@ref)`(10^{-4})))`
  a functor inheriting from [`StoppingCriterion`](@ref) indicating when to stop.
* `return_options` – (`false`) – if activated, the extended result, i.e. the
    complete [`Options`](@ref) are returned. This can be used to access recorded values.
    If set to false (default) just the optimal value `x_opt` if returned

...
and the ones that are passed to [`decorate_options`](@ref) for decorators.

# Output
* `g` – the resulting point of PSO
OR
* `options` - the options returned by the solver (see `return_options`)
"""
function particle_swarm(
    M::AbstractManifold,
    F::TF;
    n::Int=100,
    x0::AbstractVector=[random_point(M) for i in 1:n],
    kwargs...,
) where {TF}
    x_res = allocate.(x0)
    copyto!.(Ref(M), x_res, x0)
    return particle_swarm!(M, F; n=n, x0=x_res, kwargs...)
end
@doc raw"""
    patricle_swarm!(M, F; n=100, x0::AbstractVector=[random_point(M) for i in 1:n], kwargs...)

perform the particle swarm optimization algorithm (PSO), starting with the initial particle positions ``x_0``[^Borckmans2010]
in place of `x0`.

# Input
* `M` – a manifold ``\mathcal M``
* `F` – a cost function ``F:\mathcal M→ℝ`` to minimize

# Optional
* `n` - (`100`) number of random initial positions of x0
* `x0` – the initial positions of each particle in the swarm ``x_k^{(0)} ∈ \mathcal M`` for ``k = 1, \dots, n``, per default these are n [`random_point`](@ref)s

for more optional arguments, see [`particle_swarm`](@ref).
"""
function particle_swarm!(
    M::AbstractManifold,
    F::TF;
    n::Int=100,
    x0::AbstractVector=[random_point(M) for i in 1:n],
    velocity::AbstractVector=[random_tangent(M, y) for y in x0],
    inertia::Real=0.65,
    social_weight::Real=1.4,
    cognitive_weight::Real=1.4,
    stopping_criterion::StoppingCriterion=StopWhenAny(
        StopAfterIteration(200), StopWhenChangeLess(10.0^-4)
    ),
    retraction_method::AbstractRetractionMethod=ExponentialRetraction(),
    inverse_retraction_method::AbstractInverseRetractionMethod=LogarithmicInverseRetraction(),
    vector_transport_method::AbstractVectorTransportMethod=ParallelTransport(),
    return_options=false,
    kwargs..., #collect rest
) where {TF}
    p = CostProblem(M, F)
    o = ParticleSwarmOptions(
        x0,
        velocity,
        inertia,
        social_weight,
        cognitive_weight,
        stopping_criterion,
        retraction_method,
        inverse_retraction_method,
        vector_transport_method,
    )
    o = decorate_options(o; kwargs...)
    resultO = solve(p, o)
    if return_options
        return resultO
    else
        return get_solver_result(resultO)
    end
end

#
# Solver functions
#
function initialize_solver!(p::CostProblem, o::ParticleSwarmOptions)
    j = argmin([get_cost(p, y) for y in o.x])
    return o.g = deepcopy(o.x[j])
end
function step_solver!(p::CostProblem, o::ParticleSwarmOptions, iter)
    for i in 1:length(o.x)
        o.velocity[i] =
            o.inertia .* o.velocity[i] +
            o.cognitive_weight * rand(1) .*
            inverse_retract(p.M, o.x[i], o.p[i], o.inverse_retraction_method) +
            o.social_weight * rand(1) .*
            inverse_retract(p.M, o.x[i], o.g, o.inverse_retraction_method)
        xOld = o.x[i]
        o.x[i] = retract(p.M, o.x[i], o.velocity[i], o.retraction_method)
        o.velocity[i] = vector_transport_to(
            p.M, xOld, o.velocity[i], o.x[i], o.vector_transport_method
        )
        if get_cost(p, o.x[i]) < get_cost(p, o.p[i])
            o.p[i] = o.x[i]
            if get_cost(p, o.p[i]) < get_cost(p, o.g)
                o.g = o.p[i]
            end
        end
    end
end
get_solver_result(o::ParticleSwarmOptions) = o.g


text:
function lineargradient(color1::String, color2::String;
                        x0::Number = 0,
                        y0::Number = 0,
                        x2::Number = 0,
                        y2::Number = 1)

    JSFunction("""new echarts.graphic.LinearGradient($x0, $y0, $x2, $y2, [{offset: 0, color: '$color1'}, {offset: 1, color: '$color2'}])""")

end

function radialgradient(color1::String,
                        color2::String;
                        x::Number = 0.5,
                        y::Number = 0.5,
                        r::Number = 0.5)

    JSFunction("""new echarts.graphic.RadialGradient($x, $y, $r, [{offset: 0, color: '$color1'}, {offset: 1, color: '$color2'}])""")

end


text:
"""
    jl2wat(path::AbstractString)

Convert contents of a julia source code file to WebAssembly text.

# Examples
```julia
julia> using WebAssemblyText
julia> wat = jl2wat(\"example.jl\")
julia> println(wat)
```
"""
function jl2wat(path::AbstractString; debuginfo::Bool=false)
    isfile(path) || return error("No such file: $path")
    # str = open(f -> read(f, String), path)
    str = read(path, String)
    return jlstring2wat(str; debuginfo=debuginfo)
end

"""
    jlstring2wat(str::AbstractString)

Convert a string of julia source code to WebAssembly text.

# Examples
```julia
julia> using WebAssemblyText
julia> str="
hello(x) = 2.0*x
hello(1.0)
";
julia> wat = jlstring2wat(str);
julia> println(wat)
```
```wat
(module

(func \$hello (export "hello") (param \$x f32) (result f32)
( return ( f32.mul (f32.const 2.0) (local.get \$x) ) ))
)
```
"""
function jlstring2wat(str::AbstractString; debuginfo::Bool=false)
    try
        result = Base.eval(Evalscope, Meta.parse("begin $str end"))
    catch e
        return e
    end

    funcs, argtypes, imports = blockparse(str)
    SSAs = []
    WATs = []
    processed = Dict()

    # allow for nesting a(b(c())) up to some depth
    maxdepth = 100
    for _ = 1:maxdepth
        for func in keys(argtypes)
            haskey(processed, func) && continue
            ssa, wat = process(func, funcs, argtypes, imports; debuginfo=debuginfo)

            if func != Symbol("exports")
                push!(SSAs, ssa)
                push!(WATs, wat)
            end
            processed[func] = true
        end
    end

    # WATs = joinn([getimports(imports); WATs; getbuiltins(SSAs)...])
    WATs = joinn([getimports(imports); WATs; getallbuiltins()])
    memoryimport = """(memory (import "imports" "memory") 1)"""
    default_jsimports = raw"""(func $rand (import "imports" "rand") (result f32))
    (func $cos (import "imports" "cos") (param $a f32) (result f32))
    (func $log (import "imports" "log") (param $a f32) (result f32))
    (func $^ (import "imports" "^") (param $a f32) (param $b f32) (result f32))"""
    return """(module$memoryimport\n\n$default_jsimports\n$WATs\n)"""
end

getallbuiltins() = open(f -> read(f, String), joinpath(@__DIR__, "builtins.wat"))

"""
    process(func, funcs, argtypes)

Get a string with a self contained (func ) expression in .wat format

# Details
The main steps of translating a single function
- type infer func given argtypes[func]
- structure ssa
- update argtypes for other functions called in this function
- translate to wat and inline to a string
- wrap string in a function declaration
"""
function process(func, funcs, argtypes, imports; debuginfo::Bool=false)
    cinfo, Rtype = codeinfo(func, argtypes[func])
    ssa = structure(cinfo.code)

    # debuginfo && debugprint(ssa, cinfo, binfo)

    ssa = [restructure(cinfo, i, ssa, ssa[i]) for i = 1:length(ssa)]
    binfo = blockinfo(ssa)

    debuginfo && debugprint(ssa, cinfo, binfo)

    argtypes!(cinfo, argtypes, funcs, ssa)
    imports!(imports, cinfo, funcs, builtinfuncs, ssa)

    ssa = inlinessarefs(ssa)
    wat = [translate(i, cinfo, ssa[i]) for i = 1:length(ssa)]
    wat = [translate_gotos(binfo, i, wat[i]) for i = 1:length(wat)]

    wat = addparens.(wat)
    wat = addblocks(binfo, wat)
    wat = stringify(wat)
    decl = declaration(cinfo, func, argtypes[func], Rtype)
    wat = parenwrap(decl, wat)
    return ssa, wat
end

"""
    @code_wat expression

Macro for translating a single function without adding on imports and builtins.

# Examples

```julia
julia> hello(x) = 3.1*x
julia> @code_wat hello(1.2)

(func \$hello (export "hello") (param \$x f32) (result f32)
(return (f32.mul (f32.const 3.1) (local.get \$x))))
```
"""
macro code_wat(ex)
    :(code_wat($(esc(ex.args[1])), $(ex.args[2:end])))
end

function code_wat(func, args)
    types = typeof.([Base.eval(Evalscope, arg) for arg in args])
    imports = Dict()
    funcs = Dict(func => func)
    argtypes = Dict(func => types)

    ssa, wat = process(func, funcs, argtypes, imports)
    return println(wat)
end


text:
export iauNut06a
"""
IAU 2000A nutation with adjustments to match the IAU 2006
precession.

Given:
   date1,date2   double   TT as a 2-part Julian Date (Note 1)

Returned:
   dpsi,deps     double   nutation, luni-solar + planetary (Note 2)

Status:  canonical model.

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. The nutation components in longitude and obliquity are in radians
   and with respect to the mean equinox and ecliptic of date,
   IAU 2006 precession model (Hilton et al. 2006, Capitaine et al.
   2005).

   3. The function first computes the IAU 2000A nutation, then applies
   adjustments for (i) the consequences of the change in obliquity
   from the IAU 1980 ecliptic to the IAU 2006 ecliptic and (ii) the
   secular variation in the Earth's dynamical form factor J2.

   4. The present function provides classical nutation, complementing
   the IAU 2000 frame bias and IAU 2006 precession.  It delivers a
   pole which is at current epochs accurate to a few tens of
   microarcseconds, apart from the free core nutation.

Called:
   iauNut00a    nutation, IAU 2000A

References:

   Chapront, J., Chapront-Touze, M. & Francou, G. 2002,
   Astron.Astrophys. 387, 700

   Lieske, J.H., Lederle, T., Fricke, W. & Morando, B. 1977,
   Astron.Astrophys. 58, 1-16

   Mathews, P.M., Herring, T.A., Buffet, B.A. 2002, J.Geophys.Res.
   107, B4.  The MHB_2000 code itself was obtained on 9th September
   2002 from ftp//maia.usno.navy.mil/conv2000/chapter5/IAU2000A.

   Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
   Francou, G., Laskar, J. 1994, Astron.Astrophys. 282, 663-683

   Souchay, J., Loysel, B., Kinoshita, H., Folgueira, M. 1999,
   Astron.Astrophys.Supp.Ser. 135, 111

   Wallace, P.T., "Software for Implementing the IAU 2000
   Resolutions", in IERS Workshop 5.1 (2002)

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauNut06a(double date1, double date2, double *dpsi, double *deps)

function iauNut06a(date1::Real, date2::Real)
   # Allocate return values
   ref_dpsi = Ref{Float64}(0.0)
   ref_deps = Ref{Float64}(0.0)

   ccall((:iauNut06a, libsofa_c), Cvoid,
        (Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
        convert(Float64, date1), convert(Float64, date2),
        ref_dpsi, ref_deps)

   return ref_dpsi[], ref_deps[]
end


text:
module TestBasicTaskInfo

using Base.Test
using SimpleTasksTests.Utils.TestTasks

import SimpleTasks.Tasks.BasicTask
import JSON

function test_basic_info_no_id()
    info = make_valid_basic_info()
    dict = JSON.parse(JSON.json(info))
    delete!(dict, "id")
    BasicTask.Info(dict)
end

function test_basic_info_empty_id()
    info = make_valid_basic_info()
    dict = JSON.parse(JSON.json(info))
    dict["id"] = ""
    # fail ArgumentError: premature end of integer
    @test_throws ArgumentError BasicTask.Info(dict)
end

function test_basic_info_bad_id()
    info = make_valid_basic_info()
    dict = JSON.parse(JSON.json(info))
    dict["id"] = "a"
    #@test_throws ArgumentError: trying to convert "a"
    @test_throws ArgumentError BasicTask.Info(dict)
end

function test_basic_info_no_name()
    info = make_valid_basic_info()
    dict = JSON.parse(JSON.json(info))
    delete!(dict, "name")
    @test_throws KeyError BasicTask.Info(dict)
end

function test_basic_info_empty_name()
    info = make_valid_basic_info()
    dict = JSON.parse(JSON.json(info))
    dict["name"] = ""
    @test_throws ArgumentError BasicTask.Info(dict)
end

function test_basic_info_no_base_directory()
    info = make_valid_basic_info()
    dict = JSON.parse(JSON.json(info))
    delete!(dict, "base_directory")
    @test_throws KeyError BasicTask.Info(dict)
end

function test_basic_info_empty_base_directory()
    info = make_valid_basic_info()
    dict = JSON.parse(JSON.json(info))
    dict["base_directory"] = "   "
    @test_throws ArgumentError BasicTask.Info(dict)
end

function test_basic_info_no_inputs()
    info = make_valid_basic_info()
    dict = JSON.parse(JSON.json(info))
    delete!(dict, "inputs")
    @test_throws KeyError BasicTask.Info(dict)
end

function test_basic_info_empty_inputs()
    info = make_valid_basic_info()
    dict = JSON.parse(JSON.json(info))
    dict["inputs"] = []
    @test_throws ArgumentError BasicTask.Info(dict)
end

function test_basic_info_info_good()
    info = make_valid_basic_info()
    dict = JSON.parse(JSON.json(info))
    new_info = BasicTask.Info(dict)
    @test new_info != nothing
    @test new_info.id == TEST_ID
    @test new_info.name == TEST_TASK_NAME
    @test new_info.base_directory == TEST_BASE_DIRECTORY
    @test length(new_info.inputs) == length(TEST_INPUTS)
    for index in 1:length(TEST_INPUTS)
        @test new_info.inputs[index] == TEST_INPUTS[index]
    end
end

function __init__()
    test_basic_info_no_id()
    test_basic_info_empty_id()
    test_basic_info_bad_id()
    test_basic_info_no_name()
    test_basic_info_empty_name()
    test_basic_info_no_base_directory()
    test_basic_info_empty_base_directory()
    test_basic_info_no_inputs()
    test_basic_info_empty_inputs()

    test_basic_info_info_good()
end

end # module TestBasicTaskInfo

max_forks_repo_name: seung-lab/SimpleTasks


text:
@testset "entropy_thresholding" begin
    img = testimage("cameraman")
    edges, histogram = build_histogram(img,256)
    @test_throws ErrorException find_threshold(Entropy(),histogram,edges)
    @test find_threshold(Entropy(),histogram[1:256],edges) ≈ 0.76171875

    img = testimage("mandril_gray")
    edges, histogram = build_histogram(img,256)
    @test find_threshold(Entropy(),histogram[1:256],edges) ≈ 0.408517187461257

    # Same as above, but passing image directly.
    t = find_threshold(img, Entropy(); nbins = 256)
    @test t ≈ 0.408517187461257

    # There are some performance tweaks, make sure generic arrays are not broken.
    f(counts) = find_threshold(Entropy(),counts,edges)
    @test f(view(histogram, 1:256)) ==
          f(view(SVector{257}(collect(histogram)), 2:257)) ==
          f(SVector{256}(histogram[1:256])) ==
          f(histogram[1:256])
end

max_stars_repo_name: zygmuntszpak/HistogramThresholding

max_forks_repo_name: zygmuntszpak/HistogramThresholding


text:
# function that chooses the right index1, index2
# to add hybridization for Case F bad diamond
# to tree_example.jl
# warning: requires other functions not updated here
#          as this function will only run as part of
#          tests_5taxon.jl
# Claudia September 2015

# we need a different chooseEdgesGamma to control the edges1,2 chosen
# change index1, index2 according to the specific case
# warning: not real chooseEdgesGamma, the real one is in functions.jl
function chooseEdgesGamma(net::HybridNetwork)
warn("function chooseEdgesGamma is deterministic")
    index1 = 7;
    index2 = 6;
    edge1 = net.edge[index1];
    edge2 = net.edge[index2];
    gamma = rand()*0.5;
    return edge1, edge2, gamma
end

# index1=7, index2=6 => case F (bad diamond)
include("print_add.jl")
testCaseF(net)


text:
include("tabular_approximator.jl")
include("neural_network_approximator.jl")


text:
using Plots
using StatsBase
include("Dynamic.jl")

all(assn_r1[2] .≥ assn_r2[2]) ? println("✓ Rankwise dominance holds") : println("✗ Rankwise dominance fails")
diff_ = assn_r1[2] - assn_r2[2]
println("Sum-of-ranks improvement: ", sum(diff_))

dist_cmap = countmap(diff_)
rank_hist = [get(dist_cmap, i, 0) for i in 0:maximum(keys(dist_cmap))]

descr = "Two-round dynamic reassignment with $n students, $m schools, γ = $γ"

p = plot(0:maximum(keys(dist_cmap)), rank_hist, seriestype=:bar,
         c=:rebeccapurple,
         legend=false,
         title=descr, titlefontsize=11,
         xlabel="decrease in rank of assigned school", ylabel= "number of students")

display(p)

readline()

q = plot([assn_r1[2], [0, m + 1]], [assn_r2[2], [0, m + 1]],
         seriestype=[:scatter :line],
         lc=[nothing :dimgray],
         ls=[:auto :dash],
         c=[:indigo nothing],
         ma=[0.5 0],
         msc=nothing,
         legend=false,
         title=descr, titlefontsize=11,
         xlabel="rank in round 1", ylabel= "rank in round 2")

display(q)

readline()

# savefig(p, string("examples/dynamic/plot-hist.pdf"))
# savefig(p, string("examples/dynamic/plot-hist.png"))
#
# savefig(q, string("examples/dynamic/plot-scatter.pdf"))
# savefig(q, string("examples/dynamic/plot-scatter.png"))

max_stars_repo_name: maxkapur/DeferredAcceptance

max_issues_repo_name: maxkapur/DeferredAcceptance

max_forks_repo_name: maxkapur/DeferredAcceptance


text:
@testset "thermint" begin
    D = 5
    prior = MvNormal(0.5 * ones(D))
    likelihood = MvNormal(2.0 * ones(D))
    logprior(x) = logpdf(prior, x)
    loglikelihood(x) = logpdf(likelihood, x)
    alg = ThermInt(n_samples=5000)
    logZ = alg(logprior, loglikelihood, rand(prior))
    true_logZ = -0.5 * (logdet(cov(prior) + cov(likelihood)) + D * log(2π))
    @test logZ ≈ true_logZ atol=1e-1
    logZparallel =  alg(logprior, loglikelihood, rand(prior), TIParallelThreads(); progress=false)
    @test logZparallel ≈ logZ atol=1e-1
end


text:
# responsible for the interpretation of solutions
"""
An abstract type that requires to implement the basic functions needed for path following:
Homotopy `H`, its jacobian `J` (and `show` to display a GUI)
"""
abstract SystemCore

# usually displays a GUI
Base.show(::SystemCore) = error("implement!")

# return  R^(N+1)->R^N  /  R^(N+1)->R^(Nx(N+1))  functions
H(::SystemCore) = error("implement!")
J(::SystemCore) = error("implement!")

#TODO replace observer with required functions in system core?
#	interface should be clear...


text:
"""
    _pr_metafmt(level::Logging.LogLevel, _module, group, id, file, line)

MetaFormatter for ConsoleLogger for PR to adjust log message format
"""
function _pr_metafmt(level::Logging.LogLevel, _module, group, id, file, line)
    @nospecialize
    color = Logging.default_logcolor(level)
    prefix =
        "$(_module) | " *
        (level == Logging.Warn ? "Warning" : string(level)) *
        " ] :"
    suffix = ""
    Logging.Info <= level < Logging.Warn && return color, prefix, suffix
    _module !== nothing && (suffix *= "$(_module)")
    if file !== nothing
        _module !== nothing && (suffix *= " ")
        suffix *= Base.contractuser(file)
        if line !== nothing
            suffix *= ":$(isa(line, UnitRange) ? "$(first(line))-$(last(line))" : line)"
        end
    end
    !isempty(suffix) && (suffix = "@ " * suffix)

    return color, prefix, suffix
end

"""
    silence!()

Sets loglevel for PMD to :Error, silencing Info and Warn
"""
function silence!()
    return set_logging_level!(:Error)
end

"""
    reset_logging_level!()

Resets the log level to Info
"""
function reset_logging_level!()
    Logging.global_logger(_LOGGER)

    return
end

"""
    restore_global_logger!()

Restores the global logger to its default state (before PMD was loaded)
"""
function restore_global_logger!()
    Logging.global_logger(_DEFAULT_LOGGER)

    return
end

"""
    set_logging_level!(level::Symbol)

Sets the logging level for PR: :Info, :Warn, :Error, :Debug
"""
function set_logging_level!(level::Symbol)
    Logging.global_logger(_make_filtered_logger(getfield(Logging, level)))

    return
end

"""
    _make_filtered_logger(level::Logging.LogLevel)

Helper function to create the filtered logger for PMD
"""
function _make_filtered_logger(level)
    LoggingExtras.EarlyFilteredLogger(_LOGGER) do log
        if log._module == PolyhedralRelaxations && log.level < level
            return false
        else
            return true
        end
    end
end


text:
using Test
using Basics.AbstractDataStructure: DataStructure, ins!, del!, at, search, len, ordered_data_structure_p
using Basics.LinkedLists #: LinkedList, ins!, del!, at, search, len, ordered_data_structure_p
using Basics.SkipLists #: SkipList, ins!, del!, at, search, len, ordered_data_structure_p

function test_order_via_interface(ds::DataStructure{Int, Int})
    @test issorted([i for i in ds])
    inserted_zero = ins!(ds, 0)
    inserted_three = ins!(ds, 3)
    inserted_minus_two = ins!(ds, -2)
    inserted_zero = ins!(ds, 1)
    inserted_zero = ins!(ds, 0)
    @test issorted([i for i in ds])
    del!(ds, search(ds, 1)[1])
    @test issorted([i for i in ds])
end

function test_correctness_via_interface(ds::DataStructure{Int, Int})
    # Results when empty
    @test len(ds) == 0
    @test_throws KeyError at(ds, 1)
    @test search(ds, 1) == Int[]
    @test_throws KeyError del!(ds, 1)
    # Insert one 0
    first_inserted = ins!(ds, 0)
    @test first_inserted isa Int
    @test at(ds, first_inserted) == 0
    @test len(ds) == 1
    @test search(ds, 0) == [first_inserted]
    @test search(ds, 1) == Int[]
    @test_throws KeyError del!(ds, first_inserted + 1)
    # Delete it and run some tests
    @test del!(ds, first_inserted) == 0
    @test_throws KeyError at(ds, first_inserted)
    @test search(ds, 0) == Int[]
    @test len(ds) == 0
    # Insert 0, 3, -2
    inserted_zero = ins!(ds, 0)
    @test at(ds, inserted_zero) == 0
    inserted_three = ins!(ds, 3)
    @test at(ds, inserted_three) == 3
    inserted_minus_two = ins!(ds, -2)
    @test at(ds, inserted_minus_two) == -2
    @test search(ds, -2) == [inserted_minus_two]
    @test (inserted_zero, inserted_three, inserted_minus_two) isa Tuple{Int,Int,Int}
    @test len(ds) == 3
    @test search(ds, 1) == Int[]
    @test_throws KeyError del!(ds, search(ds, -2)[1] * search(ds, 0)[1] * search(ds, 3)[1])
    # Delete 3, check that 0 and -2 are left
    three = del!(ds, search(ds, 3)[1])
    @test three == 3
    @test len(ds) == 2
    @test at(ds, search(ds, -2)[1]) == -2
    @test at(ds, search(ds, 0)[1]) == 0
    @test search(ds, 3) == Int[]
    @test at(ds, search(ds, -2)[1]) == -2
    # Insert another -2, check that there are two of them
    inserted_second_minus_two = ins!(ds, -2)
    @test length(search(ds, -2)) == 2
    @test length(search(ds, 0)) == 1
    @test len(ds) == 3
    # Empty and run order tests
    if ordered_data_structure_p(ds)
        del!(ds, search(ds, -2)[1])
        del!(ds, search(ds, -2)[1])
        del!(ds, search(ds, 0)[1])
        @test len(ds) == 0
        test_order_via_interface(ds)
    end
end

@testset "Interface contracts" begin
    test_correctness_via_interface(LinkedList{Int}())
    test_correctness_via_interface(SkipList{Int}())
end


text:
"""
Test the peakedisotropicspectrum function.
"""
function testpeakedisotropicspectrum(dev::Device=CPU())
  n, L = 128, 2π
  grid = TwoDGrid(dev, n, L)
  k0, E0 = 6, 0.5
  qi = peakedisotropicspectrum(grid, k0, E0; allones=true)
  ρ, qhρ = FourierFlows.radialspectrum(rfft(qi) .* grid.invKrsq, grid)

  ρtest = ρ[ (ρ.>15.0) .& (ρ.<=17.5)]
  qhρtest = qhρ[ (ρ.>15.0) .& (ρ.<=17.5)]

  CUDA.@allowscalar isapprox(abs.(qhρtest)/abs(qhρtest[1]), (ρtest/ρtest[1]).^(-2), rtol=5e-3)
end

function testpeakedisotropicspectrum_rectangledomain()
  nx, ny = 32, 34
  Lx, Ly = 2π, 3π
  grid = TwoDGrid(nx, Lx, ny, Ly)
  k0, E0 = 6, 0.5
  qi = peakedisotropicspectrum(grid, k0, E0; allones=true)
end


text:
import Distributions

export Dists
const Dists = Distributions

function Measure(dist::Dists.Distribution{F,S}) where {F,S}
    X = Base.eltype(dist)
    DistributionMeasure{F,S,X}(dist)
end

struct DistributionMeasure{F, S, X} <: AbstractMeasure{X}
    dist :: Dists.Distribution{F, S}
end

function baseMeasure(μ::DistributionMeasure{F,S,X}) where {F, S <: Dists.Continuous, X}
    return Lebesgue(X)
end

export logdensity

function logdensity(μ::DistributionMeasure{F,S,X}, x::X) where {F, S, X}
    return Dists.logpdf(μ.dist, x)
end

export rand

function Base.rand(μ::DistributionMeasure{F,S,X}) where {F, S, X}
    return rand(μ.dist)
end

export Normal

function baseMeasure(μ::Dists.Distribution{Dists.Univariate,Dists.Continuous})
    T = eltype(μ)
    return Lebesgue(T)
end

function baseMeasure(μ::Dists.Distribution{Dists.Multivariate,Dists.Continuous})
    x = rand(μ)
    return Lebesgue(typeof(x))
end


text:
using SpecialFunctions
import NaNMath
import Calculus

# This implementation is essentially identical to the implementation in DualNumbers.jl

# force use of NaNMath functions in derivative calculations
function to_nanmath(x::Expr)
    if x.head == :call
        funsym = Expr(:.,:NaNMath,Base.Meta.quot(x.args[1]))
        return Expr(:call,funsym,[to_nanmath(z) for z in x.args[2:end]]...)
    else
        return Expr(:call,[to_nanmath(z) for z in x.args]...)
    end
end
to_nanmath(x) = x

for (funsym, exp) in Calculus.symbolic_derivatives_1arg()
    funsym == :exp && continue
    funsym == :abs2 && continue
    funsym == :inv && continue

    if isdefined(SpecialFunctions, funsym)
        @eval function SpecialFunctions.$(funsym)(z::DiffCheck{T,S,N}) where {T, S, N}
            x = value(z)
            xp = values(derivatives(z))
            DiffCheck{T, S}($(funsym)(x), StaticTuple{T,S,N}($exp .* xp))
        end
    elseif isdefined(Base, funsym)
        @eval function Base.$(funsym)(z::DiffCheck{T,S,N}) where {T, S, N}
            x = value(z)
            xp = values(derivatives(z))
            DiffCheck{T, S}($(funsym)(x), StaticTuple{T,S,N}($exp .* xp))
        end
    end
    # extend corresponding NaNMath methods
    if funsym in (:sin, :cos, :tan, :asin, :acos, :acosh, :atanh, :log, :log2, :log10,
          :lgamma, :log1p)
        funsym = Expr(:.,:NaNMath,Base.Meta.quot(funsym))
        @eval function $(funsym)(z::DiffCheck{T,S,N}) where {T, S, N}
            x = value(z)
            xp = values(derivatives(z))
            DiffCheck{T, S}($(funsym)(x), StaticTuple{T,S,N}($(to_nanmath(exp)) .*xp))
        end
    end
end

# only need to compute exp/cis once
Base.exp(z::DiffCheck{T,S,N}) where{T,S,N} = (expval = exp(value(z)); DiffCheck{T, S}(expval, StaticTuple{T,S,N}(expval .* values(derivatives(z)))))
Base.cis(z::DiffCheck{T,S,N}) where{T,S,N}= (cisval = cis(value(z)); DiffCheck{T, S}(cisval, StaticTuple{T,S,N}( (cisval*im) .* values(derivatives(z)))))
Base.exp10(z::DiffCheck{T,S,N}) where{T,S,N} = (y = exp10(value(x)); DiffCheck{T, S}(cisval, StaticTuple{T,S,N}( (y*log(10)) .* values(derivatives(z)))))

## TODO: should be generated in Calculus
Base.sinpi(z::DiffCheck{T,S,N}) where{T,S,N} = DiffCheck{T, S}(sinpi(value(z)),StaticTuple{T,S,N}(values(derivatives(z)).*(cospi(value(z))*π)))
Base.cospi(z::DiffCheck{T,S,N}) where{T,S,N} = DiffCheck{T, S}(cospi(value(z)),StaticTuple{T,S,N}(values(derivatives(z)).*(sinpi(value(z))*π)))

Base.checkindex(::Type{Bool}, inds::AbstractUnitRange, i::DiffCheck) = checkindex(Bool, inds, value(i))


text:
# Copyright 2020 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

# usage: plot_matrix_games_results.py [results directory]

using CSV
using DataFrames
using Plots
using Plots.PlotMeasures

gr()

fntsm = Plots.font("serif", pointsize=12)
fntlg = Plots.font("serif", pointsize=18)
default(titlefont=fntlg, guidefont=fntlg, tickfont=fntsm, legendfont=fntsm)

function do_plot(csv_path, threshold_for_best_fixed, y_axis_range)

  results = DataFrame(CSV.File(csv_path))
  restart_frequencies = [8, 32, 128, 512, 2048]  # Must match up with pdhg_matrix_games.py.

  # Find the fixed restart frequency for which median_residual_at_current first
  # drops below threshold_for_best_fixed.
  best_restart_frequency = nothing
  earliest_time = typemax(Int32)
  for restart_frequency in restart_frequencies
    filtered = filter(row -> row["method"] == "pdhg restart $(restart_frequency)", results)
    indices = findall(filtered[!, "median_residual_at_current"] .<= threshold_for_best_fixed)
    if length(indices) > 0 && indices[1] < earliest_time
      earliest_time = indices[1]
      best_restart_frequency = restart_frequency
    end
  end
  @show best_restart_frequency

  filtered = filter(row -> row["method"] == "pdhg", results)
  plot(filtered[!, "iteration_num"], filtered[!, "median_residual_at_current"],
       ribbon=(filtered[!, "median_residual_at_current"] .- filtered[!, "lower_range_at_current"], filtered[!, "upper_range_at_current"] .- filtered[!, "median_residual_at_current"]),
       label="No restarts", line=(1,:solid), right_margin = 4mm)

  filtered = filter(row -> row["method"] == "pdhg adaptive", results)
  plot!(filtered[!, "iteration_num"], filtered[!, "median_residual_at_avg"],
       ribbon=(filtered[!, "median_residual_at_avg"] .- filtered[!, "lower_range_at_avg"], filtered[!, "upper_range_at_avg"] .- filtered[!, "median_residual_at_avg"]),
       label="Our adaptive", line=(1,:solid))

  filtered = filter(row -> row["method"] == "pdhg restart $(best_restart_frequency)", results)
  plot!(filtered[!, "iteration_num"], filtered[!, "median_residual_at_avg"],
       ribbon=(filtered[!, "median_residual_at_avg"] .- filtered[!, "lower_range_at_avg"], filtered[!, "upper_range_at_avg"] .- filtered[!, "median_residual_at_avg"]),
       label="Best fixed ($(best_restart_frequency))", line=(1,:solid))

  # This plots all the fixed restarts.
  # for restart in restart_frequencies
  #   filtered = filter(row -> row["method"] == "pdhg restart $(restart)", results)
  #   plot!(filtered[!, "iteration_num"], filtered[!, "median_residual_at_avg"],
  #      ribbon=(filtered[!, "median_residual_at_avg"] .- filtered[!, "lower_range_at_avg"], filtered[!, "upper_range_at_avg"] .- filtered[!, "median_residual_at_avg"]),
  #      label="Fixed Restart ($(restart))", line=(1,:solid))
  # end

  yaxis!("Residual", :log10)
  if y_axis_range !== nothing
    ylims!(y_axis_range)
  end
  xaxis!("Iteration")
  savefig(replace(csv_path, ".csv" => ".pdf"))
end

do_plot(joinpath(ARGS[1], "uniform_game_100_100.csv"), 1e-4, nothing)
do_plot(joinpath(ARGS[1], "normal_game_100_100.csv"), 1e-6, (1e-10, 1.0))


text:
filter_ext(path, ext::AbstractString) = filter(fn -> splitext(fn)[2] == ext, readdir(path))
filter_ext(path, exts::Union{Tuple,AbstractArray}) =
    filter(fn -> splitext(fn)[2] in exts, readdir(path))
filter_ext(path, ext::Nothing) = readdir(path)

# Check that array order matches expectation
checkarrayorder(A, order::Tuple) = map(checkarrayorder, dims(A), order)
checkarrayorder(dim::Dimension, order::Order) =
    arrayorder(dim) == order || @warn "Array order for `$(DD.basetypeof(order))` is `$(arrayorder(dim))`, usually `$order`"

checkindexorder(A, order::Tuple) = map(checkindexorder, dims(A), order)
checkindexorder(dim::Dimension, order::Order) =
    indexorder(dim) == order || @warn "Array order for `$(DD.basetypeof(order))` is `$(indexorder(dim))`, usually `$order`"

cleankeys(keys) = Tuple(map(Symbol, keys))

noindex_to_sampled(A) = rebuild(A; dims=noindex_to_sampled(dims(A)))
function noindex_to_sampled(dims::DimTuple)
    map(dims) do d
        mode(d) isa NoIndex ? set(d, Sampled) : d
    end
end

function maybe_typemin_as_missingval(filename::String, A::AbstractGeoArray{T}) where T
    if ismissing(missingval(A))
        newmissingval = typemin(Missings.nonmissingtype(T))
        ext = splitext(filename)
        @warn "`missing` cant be written to $ext, typemin of $newmissingval used instead"
        replace_missing(A, newmissingval)
    else
        A
    end
end

# We often need to convert the locus and the mode in the same step,
# as doing it in the wrong order can give errors.
# function convert_locus_mode(M1::Type{<:IndexMode}, L1::Type{<:Locus}, dim::Dimension)
#     _convert(S1, L1, sampling(dim), locus(dim), span(dim), dim)
# end

# _convert(::Type{M1}, ::Type{L1}, mode::M2, l2::L2, span, dim) where {M1,M2<:M1,L1,L2<:L1} = dim
# _convert(::Type{M1}, ::Type{L1}, mode::M2, l2::L2, span, dim) where {M1,M2<:M1,L1,L2} = shiftlocus(L1(), dim)
# _convert(::Type{M1}, ::Type{L1}, mode::M2, l2::L2, span, dim) where {M1,M2<:M1,L1,L2} =
#     _convert_by_locus(M1, L1, mode, l2, span, dim)

# _convert_by_locus(M1, ::Type{Center}, mode, l2::Union{Start,End}, span, dim) =
#     _convert_by_mode(M1, dim)
# _convert_by_locus(M1, L1::Type{Union{Start,End}}, mode, l2::Center, span, dim) =
#     _convert_by_mode(M1, L1, dim)
# _convert_by_locus(M1, L1::Type{Start}, mode, l2::End, span, dim) =
#     convertmode(M1, shiftlocus(L1, dim))
# _convert_by_locus(M1, L1::Type{End}, mode, l2::Start, span, dim) =
#     convertmode(M1, shiftlocus(L1, dim))
# _convert_by_locus(M1, ::Type{L1}, mode, l2::L2, span, dim) where {L1,L2<:L1} =
#     convertmode(M1, dim)

# # Projected will have an accurate center point on an equal scale, but Mapped may not.
# # So we always shift the locus while in Projected mode to avoid errors.
# _convert_by_mode(::Type{Mapped}, dim) = convertmode(Mapped, shiftlocus(Center(), dim))
# _convert_by_mode(::Type{Projected}, dim) = shiftlocus(Center(), convertmode(Projected, dim))
#

_not_a_dimcol(data, dimcols::DimTuple) = _not_a_dimcol(data, map(DD.dim2key, dimcols))
_not_a_dimcol(data, dimcols::Tuple{Vararg{<:Pair}}) = _not_a_dimcol(data, map(last, dimcols))
function _not_a_dimcol(data, dimcols::Tuple{Vararg{Symbol}})
    names = Tables.columnnames(data)
    if length(names) == 0
        names = keys(first(Tables.rows(data)))
    end
    not_dim_keys = Tuple(k for k in names if !(k in dimcols))
    return  not_dim_keys
end


text:
for N in [Float64, Rational{Int}, Float32]
    # constructor
    b = BallInf(N[2, 2, 2], N(1))
    s = Singleton(N[0, 0])
    r = Dict(1 => N(4), 3 => N(0))
    r_none = Dict{Int, N}()  # no reset
    r_1 = Dict{Int, N}(1 => N(1))  # reset x1
    r_12 = Dict{Int, N}(1 => N(1), 2 => N(1))  # reset x1 and x2
    rm = ResetMap(b, r)
    # rm is the line segment in 3D from [4, 1, 0] to [4, 3, 0]

    # type
    @test isoperationtype(typeof(rm)) == true
    @test isconvextype(typeof(rm)) == true

    # dimension
    @test dim(rm) == 3

    # support vector
    @test σ(N[1, 1, 1], rm) == N[4, 3, 0]
    b1 = Ball1(N[1, 1], N(1))
    @test σ(N[1, 1], ResetMap(b1, Dict(1 => N(3)))) == N[3, 2]
    @test σ(N[1, 1], ResetMap(b1, Dict(2 => N(3)))) == N[2, 3]
    svec = σ(N[1, 1], ResetMap(s, r_1))
    @test svec == N[1, 0] && element(s) == N[0, 0]  # does not modify s

    # boundedness
    @test isbounded(rm)  # bounded set
    @test isbounded(ResetMap(Singleton(N[1, 2]), r_none))  # bounded set without resets
    @test !isbounded(ResetMap(Universe{N}(2), r_1))  # unbounded set without enough resets
    @test isbounded(ResetMap(Universe{N}(2), r_12))  # unbounded set with enough resets

    # an_element function
    an_element(rm)
    elem = an_element(ResetMap(s, r_1))
    @test elem == N[1, 0] && element(s) == N[0, 0]  # does not modify s

    # isempty
    @test !isempty(rm)

    # getter for affine map
    @test matrix(rm) == Diagonal(N[0, 1, 0])
    @test vector(rm) == sparsevec([1], N[4], 3)

    # constraints_list
    if test_suite_polyhedra
        p = HPolytope(constraints_list(rm))
        @test N[4, 1, 0] ∈ p && N[4, 3, 0] ∈ p && N[2, 2, 2] ∉ p
    end
    # constraints_list of a hyperrectangular set
    @test ispermutation(constraints_list(rm),
                        [HalfSpace(N[1, 0, 0], N(4)),
                         HalfSpace(N[-1, 0, 0], N(-4)),
                         HalfSpace(N[0, 1, 0], N(3)),
                         HalfSpace(N[0, -1, 0], N(-1)),
                         HalfSpace(N[0, 0, 1], N(0)),
                         HalfSpace(N[0, 0, -1], N(0))])

    # constraints list of a lazy set
    if test_suite_polyhedra && N == Float64
        rm = ResetMap(b, r)
        rm_id = ResetMap(Matrix(one(N)*I, 3, 3) * b, r)
        X = HPolytope(constraints_list(rm))
        Y = HPolytope(constraints_list(rm_id))
        @test isequivalent(X, Y)
    end
end

for N in [Float64]
    b = BallInf(N[2, 2, 2], N(1))
    r = Dict(1 => N(4), 3 => N(0))
    r_none = Dict{Int, N}()  # no reset
    r_1 = Dict{Int, N}(1 => N(1))  # reset x1
    r_12 = Dict{Int, N}(1 => N(1), 2 => N(1))  # reset x1 and x2
    rm = ResetMap(b, r)

    # boundedness
    P = HPolyhedron([HalfSpace(N[1, 1], N(1))])
    @test !isbounded(ResetMap(P, r_none))  # unbounded set without resets
    @test !isbounded(ResetMap(P, r_1))  # unbounded set without enough resets
    @test isbounded(ResetMap(P, r_12))  # unbounded set, but captured by resets
    @test !isbounded(ResetMap(HPolyhedron([HalfSpace(N[1, 1], N(1))]), r_none))  # unbounded set without resets
    @test !isbounded(ResetMap(HPolyhedron([HalfSpace(N[1, 1], N(1))]), r_1))  # unbounded set without enough resets
    @test isbounded(ResetMap(HPolyhedron([HalfSpace(N[1, 1], N(1))]), r_12))  # unbounded set, but captured by resets

    # intersection
    b2 = BallInf(N[4, 2, 0], N(1))
    for cap in [intersection(b2, rm), intersection(rm, b2)]
        @test N[4, 1, 0] ∈ cap && N[4, 3, 0] ∈ cap && N[3, 1, 0] ∉ cap &&
              N[3, 1, 2] ∉ cap
    end
end