text:
import BlackBoxOptim
import Optim
import NLopt

"Perform WCR minimization, through the SLiSe framework"
function wise( ws, z1, a1, G_user;
            G_in = 1,
            G_out = 2,
            outpath = nothing,
            lb = nothing,
            refine = true,

            # To disable parts of oriented restart
            gappingEnabled = true,
            weightsEnabled = true,
            lastIntEnabled = true
        )
    p      = length( z1 )
    x0     = SLiSeFilters.startingPointWS(ws) # the starting weight function
    G      = sqrt( G_user )                   # shifting for better filters
    mintol = Inf

    # Use bounds within minimization
    minalgo = "optim"
    if lb != nothing
        minalgo = "lbfgsb"
    end

    # Maximum iterations of minimizations
    maxiter          = 200
    every            = 10
    runs             = 20
    inner_max_eval   = 300
    inner_max_eval_2 = 1000

    # Strech variables for parameters in blackbox optimization
    weight_multiplier    = 10.
    interval_multplier   = 3.

    # Remeber set-up, to be able to change and reset it easily
    intv_free_prev       = copy(ws.intv_free)
    intv_free_count_prev = copy(ws.intv_free_count)
    wght_free_prev       = copy(ws.wght_free)
    wght_free_count_prev = copy(ws.wght_free_count)

    # Monitor progress
    prevWCR    = Inf # monitor improvement step by step
    prevItWcr  = Inf # monitor improvement iteration by iteration
    tol        = 1e-9
    tolReached = false

    # Objective function
    function f(value::Vector)
        c      = SLiSeFilters.updateWS!(ws, abs.(value))
        z3, a3 = SLiSeFilters.optimise( minalgo, z1, a1, c, lb )

        # flush output to console
        if VERSION >= v"0.7.0"
            flush(stdout)
        else
            flush(STDOUT)
        end

        return wcr( z3, a3, G, mintol = mintol, lowerGap = G );
    end

    # For oriented restart
    function callBlackBox(ws, z1, a1, prevWCR, searchRange, i; traceMode=:compact, evals=inner_max_eval, method=:adaptive_de_rand_1_bin_radiuslimited) # method is default method in BlackBoxOptim
        x0  = SLiSeFilters.startingPointWS(ws)
        res = BlackBoxOptim.bboptimize(f; SearchRange = searchRange, NumDimensions = length(searchRange), MaxFuncEvals = evals, Method=method, TraceMode=traceMode, TraceInterval=20.0);
        bf  = BlackBoxOptim.best_fitness(res)
        bs  = BlackBoxOptim.best_candidate(res)

        # Paramter Gap: Evaluation and reset
        if (bf < prevWCR)
            @printf "[%d] Setting to new min: %.16e\n" i bf
            prevWCR = bf

            c       = SLiSeFilters.updateWS!( ws, abs.(bs) )
            z1,a1   = SLiSeFilters.optimise( minalgo, z1, a1, c, lb )
        else
            SLiSeFilters.updateWS!(ws, x0) # to call modifier before resetting it
        end

        return z1,a1,prevWCR
    end

    for i = 1:runs
        println("Starting.")

        ws.wght_free       = falses(length(ws.wght_free))
        ws.intv_free       = falses(length(ws.intv_free))
        ws.intv_free_count = 0
        ws.wght_free_count = 0

        # Oriented restart/Minimize gap position
        if (gappingEnabled)
            ws.intv_free_count = 1

            @printf "Performing gap minimization [1/2]..\n"
            ws.intv_free[G_in] = true
            z1,a1,prevWCR = callBlackBox(ws, z1, a1, prevWCR, [(G_user, 1.)], i)
            ws.intv_free[G_in]  = false

            @printf "Performing gap minimization [2/2]..\n"
            ws.intv_free[G_out] = true
            z1,a1,prevWCR = callBlackBox(ws, z1, a1, prevWCR, [(1., G_user^-1)], i)
            ws.intv_free[G_out] = false

            ws.intv_free_count = 0
        end

        # Oriented restart/Loop through weights (except gaps)
        if (weightsEnabled)
            ws.wght_free_count  = 1

            for j = 1:length(ws.wght_free)
                @printf "Performing weight (%d) minimization..\n" j
                ws.wght_free[j] = true
                x0  = SLiSeFilters.startingPointWS(ws)
                z1,a1,prevWCR = callBlackBox(ws, z1, a1, prevWCR, [(x0[1]/weight_multiplier, weight_multiplier*x0[1])], i)
                ws.wght_free[j] = false
            end

            ws.wght_free[end]  = false
            ws.wght_free_count = 0
        end

        # Oriented restart/Special, dependent minimization for last interval
        if (lastIntEnabled)
            ws.wght_free_count = 1
            ws.intv_free_count = 1
            ws.wght_free[end]  = true
            ws.intv_free[end]  = true

            @printf "Performing last interval minimization..\n"
            x0  = SLiSeFilters.startingPointWS(ws)
            z1,a1,prevWCR = callBlackBox(ws, z1, a1, prevWCR, [(ws.intv[end-1], interval_multplier*ws.intv[end]), (x0[2]/weight_multiplier, weight_multiplier*x0[2])], i, evals=inner_max_eval_2, method=:resampling_memetic_search)
        end

        # Reset free variables
        ws.intv_free       = intv_free_prev
        ws.wght_free       = wght_free_prev
        ws.intv_free_count = intv_free_count_prev
        ws.wght_free_count = wght_free_count_prev

        # Do Nelder-Mead
        x0  = SLiSeFilters.startingPointWS(ws)
        res = Optim.optimize(f,
	                         x0,
	                         Optim.NelderMead( ),
	                         Optim.Options(iterations        = maxiter,
	                                       g_tol             = 1e-10, # =√(Σ(yᵢ-ȳ)²)/n, the conv crit for NelderMead
	                                       show_every        = every,
	                                       allow_f_increases = true, # in the end, minimizer is returned
	                                       show_trace        = true
	                                      )
	                        )

        # Evaluation and reset
        if (res.minimum < prevWCR)
            @printf "[%d] Setting to new min: %.16e\n" i res.minimum
            prevWCR = res.minimum

            c       = SLiSeFilters.updateWS!( ws, abs.(res.minimizer) )
            z1,a1   = SLiSeFilters.optimise( minalgo, z1, a1, c, lb )
        else
            SLiSeFilters.updateWS!(ws, x0) # reset search space
        end

        # Cache the results
        if (outpath != nothing)
            writedlm( string(outpath,"/wcr_",p,"_",G_user,"_lb",lb,"_filter_it",i,"_wise.dat"),([z1 a1]/G)',',') # comma-separated, for JULIA v0.7
            open( string(outpath,"/wcr_",p,"_",G_user,"_lb",lb,"_filter_it",i,"_wise.dat"), "a") do f # append
                writecsv(f, ([real(z1) imag(z1) real(a1) imag(a1)]/G)') # comma-separated, for JULIA v0.6
            end
            writedlm(string(outpath,"/wcr_",p,"_",G_user,"_lb",lb,"_weight_it",i,"_wise.dat"),[ws.intv ws.wght]',',')
        end

        # Validate change in WCR
        resd = abs( prevItWcr - prevWCR )/prevWCR
        @printf "[%d] CurrentMin: %.16e (internal)\t%.16e (real)" i prevWCR wcr( z1/G, a1/G, G_user ) # It holds that wcr( z1, a1, G, lowerGap = G ) = wcr( z1/G, a1/G, G_user, lowerGap = 1. )
        if (resd < tol)
            # If no progress in two subsequent runs, finish
            if (tolReached)
                @printf "\tEND. TOL reached.\n"
                break
            else
                @printf "\tTOL reached. Continuing..\n"
                tolReached = true
            end
        else
            tolReached = false
            if (i == runs)
                @printf "\tEND. MAXITER reached."
            end
            @printf "\n"
        end
        prevItWcr = prevWCR
    end

    if (!refine)
        return z1, a1
    end

    @printf "Performing refinement step..\n"
    return wiseRefine( z1, a1, G_user, outpath = outpath, lb = lb )
end

"Perform WCR minimization, without SLiSe"
function wiseRefine( z1, a1, G_user; lb = nothing, outpath = nothing )
    p       = length( z1 )
    maxiter = 100000
    every   = 100
    mintol  = Inf

    # Transformation from C to R
    toR(x)     = [real(x); imag(x)]
    fromR(x)   = x[1:Int(length(x)/2)] + im*x[Int(length(x)/2) + 1:end]
    x0         = toR([z1; a1])

    function g(x)
        # flush output to console
        if VERSION >= v"0.7.0"
            flush(stdout)
        else
            flush(STDOUT)
        end
        return wcr(split_(fromR(x))...,G_user, mintol = mintol )
    end

    count = 0
    function f2(value::Vector, grad::Vector)
        count::Int += 1

        y = wcr(split_(fromR(value))...,G_user, mintol = mintol )

        if (count % every == 0)
         @printf "%6d   %14e\n" count y
        end

        return y
    end

    # Do the optimization
    z1, a1 = nothing, nothing

    if (lb != nothing)
        n     = length(x0)
        lower = [ones(div(n,2))*-Inf;ones(div(n,4))*lb;ones(div(n,4))*-Inf]
        upper = [ones(n) * Inf]

        opt = NLopt.Opt(:LN_NELDERMEAD, length(x0))
        NLopt.lower_bounds!( opt, lower )
        NLopt.maxeval!(opt, maxiter)
        NLopt.ftol_rel!(opt, 1e-9)
        NLopt.min_objective!(opt, f2)

        # Perform NLOpt Minimization
        @printf "Evals    Function value    \n"
        @printf "------   --------------    \n"
        (minf,minx,reason) = NLopt.optimize(opt, x0)

        z1,a1 = split_(fromR(minx))
    else
        res = Optim.optimize(g,
                   x0,
                   Optim.NelderMead(),
                   Optim.Options(iterations        = maxiter,
                                 g_tol             = 1e-13, # =√(Σ(yᵢ-ȳ)²)/n, the conv crit for NelderMead
                                 show_every        = every,
                                 allow_f_increases = true, # in the end, minimizer is returned
                                 show_trace        = true
                                 )
                   )

        # Check for convergence
        if !Optim.converged(res)
            @printf "Attention: No convergence reached after %d steps!\n" maxiter
        end

        # Obtain the results
        z1,a1 = split_(fromR(res.minimizer))
    end

    # Output the results
    @printf "[Refined] CurrentMin: %.16e\n" wcr( z1, a1, G_user )

    # Save the results to disk
    if (outpath != nothing)
        writedlm( string(outpath,"/wcr_",p,"_",G_user,"_lb",lb,"_filter_wise_refined.dat"),([z1 a1])',',') # comma-separated, for JULIA v0.7
        open( string(outpath,"/wcr_",p,"_",G_user,"_lb",lb,"_filter_wise_refined.dat"), "a") do f # append
            writecsv(f, ([real(z1) imag(z1) real(a1) imag(a1)])') # comma-separated, for JULIA v0.6
        end
    end

    return z1,a1
end


text:
import Splittings: @Lie, @Strang, @TripleJump, @Order6
"""

    Implements split operators for linear pendulum

    Solve linear pendulum problem: ``\\frac{delta1}{dt} = v ``,

    ``\\frac{delta2}{dt} = - ω^2 x``. The exact solution is
    ``x(t)=  x(0)   cos(ωt) + \frac{v(0)}{ω}sin(ωt),```
    ``v(t)= -x(0) ω sin(ωt) + v(0)cos(ωt) ``

"""
function exact( t_final, ω, x0, v0 )
    x = x0 * cos( ω * t_final ) + (v0 / ω) * sin( ω * t_final )
    v = -x0 * ω * sin( ω* t_final ) + v0 * cos( ω * t_final )
    return x, v
end

function check_order(do_split_steps::Function,
		     steps_fine::Int64,
		     expected_order::Int )

    ω  = 2.0
    x0 = 1.0       # initial x  for order checking
    v0 = 2.0       # initial v  for order checking
    t_final = 1.0  # final time for order checking

    x_exact, v_exact = exact( t_final, ω, x0, v0 )

    # do iterations with smallest time step
    dt = t_final/steps_fine
    number_time_steps = steps_fine

    x, v = do_split_steps((x0, v0), dt, number_time_steps, ω)

    # compute  mean square error
    error0 = sqrt( (x - x_exact)^2 + (v - v_exact)^2 )

    # do iterations with middle time step
    dt = 2dt
    number_time_steps = number_time_steps ÷ 2

    x, v = do_split_steps((x0, v0), dt, number_time_steps, ω)

    # compute mean square error
    error1 = sqrt( (x - x_exact)^2 + (v - v_exact)^2 )

    # do iterations with largest time step
    dt = 2dt
    number_time_steps = number_time_steps ÷ 2

    x, v = do_split_steps((x0, v0), dt, number_time_steps, ω)

    # compute  mean square error
    error2 = sqrt( (x - x_exact)^2 + (v - v_exact)^2 )

    # compute order
    order1 = log(error1/error0)/log(2.0)
    order2 = log(error2/error1)/log(2.0)

    if (    (abs(order1-expected_order) > 5.e-2)
         || (abs(order2-expected_order) > 5.e-2))

       println( "error coarse = $error2")
       println( "error middle = $error1")
       println( "      order (coarse/middle) = $order2")
       println( "error fine   = $error0")
       println( "      order (middle/fine)   = $order1")

       return false

    end

    true

end

" First operator of splitting for linear pendulum"
function push_x!( x, v, dt )
    x .= x .+ v * dt
end

" Second operator of splitting for linear pendulum"
function push_v!( x, v, dt, ω)
    v .= v .- ω^2 * x * dt
end

@testset "Splitting Operators macros" begin

    function do_steps_lie_tv( start::Tuple{Float64,Float64},
    		   dt::Float64, nsteps::Int64, ω::Float64 )
	x = zeros(Float64, 1)
	v = zeros(Float64, 1)
	x[1], v[1] = start
        for i = 1:nsteps
            @Lie push_x!(x,v,dt) push_v!(x,v,dt,ω)
        end
	x[1], v[1]
    end

    @test check_order(do_steps_lie_tv, 200, 1)

    function do_steps_lie_vt( start::Tuple{Float64,Float64},
    		   dt::Float64, nsteps::Int64, ω::Float64 )
	x = zeros(Float64, 1)
	v = zeros(Float64, 1)
	x[1], v[1] = start
        for i = 1:nsteps
	    @Lie push_v!(x,v,dt,ω) push_x!(x,v,dt)
        end
	x[1], v[1]
    end

    @test check_order(do_steps_lie_vt, 200, 1)

    function do_steps_strang_tvt( start::Tuple{Float64,Float64},
    		   dt::Float64, nsteps::Int64, ω::Float64 )
	x = zeros(Float64, 1)
	v = zeros(Float64, 1)
	x[1], v[1] = start
        for i = 1:nsteps
            @Strang push_x!(x,v,dt) push_v!(x,v,dt,ω)
        end
	x[1], v[1]
    end

    @test check_order(do_steps_strang_tvt, 100, 2)

    function do_steps_strang_vtv( start::Tuple{Float64,Float64},
    		   dt::Float64, nsteps::Int64, ω::Float64 )
	x = zeros(Float64, 1)
	v = zeros(Float64, 1)
	x[1], v[1] = start
        for i = 1:nsteps
            @Strang push_v!(x,v,dt,ω) push_x!(x,v,dt)
        end
	x[1], v[1]
    end

    @test check_order(do_steps_strang_vtv, 100, 2)

    function do_steps_triple_jump_tvt( start::Tuple{Float64,Float64},
    		   dt::Float64, nsteps::Int64, ω::Float64 )
	x = zeros(Float64, 1)
	v = zeros(Float64, 1)
	x[1], v[1] = start
        for i = 1:nsteps
            @TripleJump push_x!(x,v,dt) push_v!(x,v,dt,ω)
        end
	x[1], v[1]
    end

    @test check_order(do_steps_triple_jump_tvt, 64, 4)

    function do_steps_triple_jump_vtv( start::Tuple{Float64,Float64},
    		   dt::Float64, nsteps::Int64, ω::Float64 )
	x = zeros(Float64, 1)
	v = zeros(Float64, 1)
	x[1], v[1] = start
        for i = 1:nsteps
            @TripleJump push_v!(x,v,dt,ω) push_x!(x,v,dt)
        end
	x[1], v[1]
    end

    @test check_order(do_steps_triple_jump_vtv, 64, 4)

    function do_steps_order6_tvt( start::Tuple{Float64,Float64},
    		   dt::Float64, nsteps::Int64, ω::Float64 )
	x = zeros(Float64, 1)
	v = zeros(Float64, 1)
	x[1], v[1] = start
        for i = 1:nsteps
            @Order6 push_x!(x,v,dt) push_v!(x,v,dt,ω)
        end
	x[1], v[1]
    end

    @test check_order(do_steps_order6_tvt, 20, 6)

    function do_steps_order6_vtv( start::Tuple{Float64,Float64},
    		   dt::Float64, nsteps::Int64, ω::Float64 )
	x = zeros(Float64, 1)
	v = zeros(Float64, 1)
	x[1], v[1] = start
        for i = 1:nsteps
            @Order6 push_v!(x,v,dt,ω) push_x!(x,v,dt)
        end
	x[1], v[1]
    end

    @test check_order(do_steps_order6_vtv, 20, 6)

end

ext: jl


text:
using WordCloud

WordCloud.wordcloud(𝑓::FeatureSet, args...; kwargs...) = wordcloud(string.(getnames(𝑓)), args...; kwargs...)
"""
We want to take a FeatureSet and and vector of weights for a PC and produce a word cloud
"""
function lowdimwordcloud(𝑓::FeatureSet, weights::AbstractVector; path="../Figures/$(rand(UInt16)).svg")
    @assert length(𝑓) == length(weights)
    wc = wordcloud(𝑓, weights; colors=:seaborn_dark, angles = 0, maskshape=squircle, outline=false, backgroundcolor = (1, 1, 1), font="Arial Medium", masksize=(1000, 400))
    generate!(wc)
    paint(wc, path)
    return wc
end

function lowdimwordcloud(𝑓::FeatureSet, M, PCs, args...; kwargs...)
    weights = sum(abs.(featureweights(M, PCs)), dims=2)[:]
    lowdimwordcloud(𝑓, weights, args...; kwargs...)
end

export lowdimwordcloud


text:
using MLBasedESC
using Flux
using LinearAlgebra
import DiffEqFlux
using DiffEqFlux: FastChain, FastDense
using ReverseDiff

const N = 2
const M⁻¹ = inv(diagm([0.1, 0.2]))
const G = [-1.0, 1.0]
const G⊥ = [1.0 1.0]

V(q) = [ 10.0*(cos(q[1]) - 1.0) ]
MLBasedESC.jacobian(::typeof(V), q) = [-10.0*sin(q[1]), zero(eltype(q))]

inmap(q,::Any=nothing) = [cos(q[1]); sin(q[1]); cos(q[2]); sin(q[2])]
function MLBasedESC.jacobian(::typeof(inmap), q,::Any=nothing)
    qbar = inmap(q)
    [
        -qbar[2] 0
         qbar[1] 0
         0 -qbar[4]
         0  qbar[3]
    ]
end

function runalltests()
    test1_chain()
    test1_fastchain()
    test1_poly()
    # test2_chain()
    # test2_fastchain()
    # test2_poly()
    test3_chain()
    test3_fastchain()
    test3_poly()
    test4_chain()
    test4_fastchain()
    test4_poly()
end

#==============================================================================
Test 1 ===> constant M; constant Md; neural-net Vd; no input mapping; no J2
==============================================================================#
function test1_chain()
    q = rand(2)
    Md⁻¹ = inv(diagm([0.1, 0.2]))
    Vd = Chain(
        Dense(2, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 1, square),
    )
    P = IDAPBCProblem(2,M⁻¹,Md⁻¹,V,Vd,G,G⊥)
    L1 = PDELossPotential(P)
    L2 = PDELossKinetic(P)
    L1(q)
    L2(q)
    θ = Flux.params(MLBasedESC.trainable(P)...)
    dL1(x) = Flux.gradient(()->L1(x), θ)
    dL2(x) = Flux.gradient(()->L2(x), θ)
    dL1(q)[ θ[rand(1:length(θ))] ]
    dL2(q)[ θ[rand(1:length(θ))] ]
end
function test1_fastchain()
    q = rand(2)
    Md⁻¹ = inv(diagm([0.1, 0.2]))
    Vd = FastChain(
        FastDense(2, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 1, square),
    )
    P = IDAPBCProblem(2,M⁻¹,Md⁻¹,V,Vd,G,G⊥)
    L1 = PDELossPotential(P)
    L2 = PDELossKinetic(P)
    θ = paramstack(P)
    L1(q,θ)
    L2(q,θ)
    dL1(q,ps) = ReverseDiff.gradient(_2->L1(q,_2), ps)
    dL2(q,ps) = ReverseDiff.gradient(_2->L2(q,_2), ps)
    unstack(P, dL1(q,θ))
    unstack(P, dL2(q,θ))
end
function test1_poly()
    q = rand(2)
    Md⁻¹ = inv(diagm([0.1, 0.2]))
    Vd = SOSPoly(2, 1:2)
    P = IDAPBCProblem(2,M⁻¹,Md⁻¹,V,Vd,G,G⊥)
    L1 = PDELossPotential(P)
    L2 = PDELossKinetic(P)
    θ = paramstack(P)
    L1(q,θ)
    L2(q,θ)
    dL1(q,ps) = ReverseDiff.gradient(_2->L1(q,_2), ps)
    dL2(q,ps) = ReverseDiff.gradient(_2->L2(q,_2), ps)
    unstack(P, dL1(q,θ))
    unstack(P, dL2(q,θ))
end

#==============================================================================
Test 2 ===> constant M; constant Md; neural-net Vd; no input mapping; yes J2
==============================================================================#
function test2_chain()
    q = rand(2)
    Md⁻¹ = inv(diagm([0.1, 0.2]))
    Vd = Chain(
        Dense(2, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 1, square),
    )
    J2 = Tuple(
        Chain(
            Dense(2, 5, elu),
            Dense(5, 4),
        ) |> makeskewsym for _=1:2
    )
    P = IDAPBCProblem(2,M⁻¹,Md⁻¹,V,Vd,J2,G,G⊥)
    L1 = PDELossPotential(P)
    L2 = PDELossKinetic(P)
    L1(q)
    L2(q)
    θ = Flux.params(MLBasedESC.trainable(P)...)
    dL1(x) = Flux.gradient(()->L1(x), θ)
    dL2(x) = Flux.gradient(()->L2(x), θ)
    dL1(q)[ θ[rand(1:length(θ))] ]
    dL2(q)[ θ[rand(1:length(θ))] ]
end
function test2_fastchain()
    q = rand(2)
    Md⁻¹ = inv(diagm([0.1, 0.2]))
    Vd = FastChain(
        FastDense(2, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 1, square),
    )
    J2 = Tuple(
        FastChain(
            FastDense(2, 5, elu),
            FastDense(5, 4),
        ) |> makeskewsym for _=1:2
    )
    P = IDAPBCProblem(2,M⁻¹,Md⁻¹,V,Vd,J2,G,G⊥)
    L1 = PDELossPotential(P)
    L2 = PDELossKinetic(P)
    θ = paramstack(P)
    L1(q,θ)
    L2(q,θ)
    dL1(q,ps) = ReverseDiff.gradient(_2->L1(q,_2), ps)
    dL2(q,ps) = ReverseDiff.gradient(_2->L2(q,_2), ps)
    unstack(P, dL1(q,θ))
    unstack(P, dL2(q,θ))
end
function test2_poly()
    q = rand(2)
    Md⁻¹ = inv(diagm([0.1, 0.2]))
    Vd = SOSPoly(2, 1:2)
    J2 = Tuple(
        FastChain(
            FastDense(2, 5, elu),
            FastDense(5, 4),
        ) |> makeskewsym for _=1:2
    )
    P = IDAPBCProblem(2,M⁻¹,Md⁻¹,V,Vd,J2,G,G⊥)
    L1 = PDELossPotential(P)
    L2 = PDELossKinetic(P)
    θ = paramstack(P)
    L1(q,θ)
    L2(q,θ)
    dL1(q,ps) = ReverseDiff.gradient(_2->L1(q,_2), ps)
    dL2(q,ps) = ReverseDiff.gradient(_2->L2(q,_2), ps)
    unstack(P, dL1(q,θ))
    unstack(P, dL2(q,θ))
end

#==============================================================================
Test 3 ===> constant M; neural-net Md; neural-net Vd; no input mapping; no J2
==============================================================================#
function test3_chain()
    q = rand(2)
    Md⁻¹ = Chain(
        Dense(2, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 4),
    ) |> makeposdef
    Vd = Chain(
        Dense(2, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 1, square),
    )
    P = IDAPBCProblem(2,M⁻¹,Md⁻¹,V,Vd,G,G⊥)
    L1 = PDELossPotential(P)
    L2 = PDELossKinetic(P)
    L1(q)
    L2(q)
    θ = Flux.params(MLBasedESC.trainable(P)...)
    dL1(x) = Flux.gradient(()->L1(x), θ)
    dL2(x) = Flux.gradient(()->L2(x), θ)
    dL1(q)[ θ[rand(1:length(θ))] ]
    dL2(q)[ θ[rand(1:length(θ))] ]
end
function test3_fastchain()
    q = rand(2)
    Md⁻¹ = FastChain(
        FastDense(2, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 4),
    ) |> makeposdef
    Vd = FastChain(
        FastDense(2, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 1, square),
    )
    P = IDAPBCProblem(2,M⁻¹,Md⁻¹,V,Vd,G,G⊥)
    L1 = PDELossPotential(P)
    L2 = PDELossKinetic(P)
    θ = paramstack(P)
    L1(q,θ)
    L2(q,θ)
    dL1(q,ps) = ReverseDiff.gradient(_2->L1(q,_2), ps)
    dL2(q,ps) = ReverseDiff.gradient(_2->L2(q,_2), ps)
    unstack(P, dL1(q,θ))
    unstack(P, dL2(q,θ))
end
function test3_poly()
    q = rand(2)
    Md⁻¹ = FastChain(
        FastDense(2, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 4),
    ) |> makeposdef
    Vd = SOSPoly(2, 1:2)
    P = IDAPBCProblem(2,M⁻¹,Md⁻¹,V,Vd,G,G⊥)
    L1 = PDELossPotential(P)
    L2 = PDELossKinetic(P)
    θ = paramstack(P)
    L1(q,θ)
    L2(q,θ)
    dL1(q,ps) = ReverseDiff.gradient(_2->L1(q,_2), ps)
    dL2(q,ps) = ReverseDiff.gradient(_2->L2(q,_2), ps)
    unstack(P, dL1(q,θ))
    unstack(P, dL2(q,θ))
end

#==============================================================================
Test 4 ===> constant M; neural-net Md; neural-net Vd; no input mapping; yes J2
==============================================================================#
function test4_chain()
    q = rand(2)
    Md⁻¹ = Chain(
        Dense(2, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 4),
    ) |> makeposdef
    Vd = Chain(
        Dense(2, 10, elu),
        Dense(10, 5, elu),
        Dense(5, 1, square),
    )
    J2 = Tuple(
        Chain(
            Dense(2, 5, elu),
            Dense(5, 4),
        ) |> makeskewsym for _=1:2
    )
    P = IDAPBCProblem(2,M⁻¹,Md⁻¹,V,Vd,J2,G,G⊥)
    L1 = PDELossPotential(P)
    L2 = PDELossKinetic(P)
    L1(q)
    L2(q)
    θ = Flux.params(MLBasedESC.trainable(P)...)
    dL1(x) = Flux.gradient(()->L1(x), θ)
    dL2(x) = Flux.gradient(()->L2(x), θ)
    dL1(q)[ θ[rand(1:length(θ))] ]
    dL2(q)[ θ[rand(1:length(θ))] ]
end
function test4_fastchain()
    q = rand(2)
    Md⁻¹ = FastChain(
        FastDense(2, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 4),
    ) |> makeposdef
    Vd = FastChain(
        FastDense(2, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 1, square),
    )
    J2 = Tuple(
        FastChain(
            FastDense(2, 5, elu),
            FastDense(5, 4),
        ) |> makeskewsym for _=1:2
    )
    P = IDAPBCProblem(2,M⁻¹,Md⁻¹,V,Vd,J2,G,G⊥)
    L1 = PDELossPotential(P)
    L2 = PDELossKinetic(P)
    θ = paramstack(P)
    L1(q,θ)
    L2(q,θ)
    dL1(q,ps) = ReverseDiff.gradient(_2->L1(q,_2), ps)
    dL2(q,ps) = ReverseDiff.gradient(_2->L2(q,_2), ps)
    unstack(P, dL1(q,θ))
    unstack(P, dL2(q,θ))
end
function test4_poly()
    q = rand(2)
    Md⁻¹ = FastChain(
        FastDense(2, 10, elu),
        FastDense(10, 5, elu),
        FastDense(5, 4),
    ) |> makeposdef
    Vd = SOSPoly(2,1:2)
    J2 = Tuple(
        FastChain(
            FastDense(2, 5, elu),
            FastDense(5, 4),
        ) |> makeskewsym for _=1:2
    )
    P = IDAPBCProblem(2,M⁻¹,Md⁻¹,V,Vd,J2,G,G⊥)
    L1 = PDELossPotential(P)
    L2 = PDELossKinetic(P)
    θ = paramstack(P)
    L1(q,θ)
    L2(q,θ)
    dL1(q,ps) = ReverseDiff.gradient(_2->L1(q,_2), ps)
    dL2(q,ps) = ReverseDiff.gradient(_2->L2(q,_2), ps)
    unstack(P, dL1(q,θ))
    unstack(P, dL2(q,θ))
end


text:
module SuperCowPowers

export cowsay, cowthink
export moo, moomoo, moomoomoo

include("Cowsay.jl")
using .Cowsay: cowsay, cowthink

include("AptMoo.jl")
using .AptMoo: moo, moomoo, moomoomoo

end # module


text:
using ComoniconTypes
using Test
using Faker

desc = Description(Faker.text())
@test endswith(desc.brief, ".") # brief should be a sentence

arg = Argument(;name="arg", type=Int)

macro test_show(mime, ex)
    Meta.isexpr(ex, :block) || error("expect begin ... end")
    ret = Expr(:block)
    for each in ex.args
        if Meta.isexpr(each, :call) && each.args[1] === :in
            @gensym buf object
            push!(ret.args, :($buf = IOBuffer()))
            push!(ret.args, :($object = $(each.args[3])))
            push!(ret.args, :(show($buf, $mime(), $object)))
            push!(ret.args, :(@test occursin($(each.args[2]), String(take!($buf)))))
        else
            push!(ret.args, each)
        end
    end
    return esc(ret)
end

@testset "convert(Description, ...)" begin
    @test Description(nothing) == Description()
    @test convert(Description, nothing) == Description()
    @test convert(Description, "nothing") == Description("nothing")
    @test convert(Description, split("abcd. efdasdas.", '.')[1]) == Description("abcd")
end

@test_show MIME"text/plain" begin
    "<arg>" in Argument(;name="arg")
    "[arg...]" in Argument(;name="arg", vararg=true)
    "<arg::Int64>" in Argument(;name="arg", type=Int)
    "--option-a" in Option(;sym=:option_a)
    "--option-a <hint>" in Option(;sym=:option_a, hint="hint")
    "--flag-a" in Flag(;sym=:flag_a)
end

leaf = LeafCommand(;
    fn=identity,
    name="leaf",
    args=[Argument(;name="arg", description=Faker.text())],
    flags=Dict(
        "flag-a" => Flag(;
            sym=:flag_a,
            description="flag a."
        ),
        "flag-b" => Flag(;
            sym=:flag_b,
            description="flag b."
        )
    ),
    description=Faker.text(),
)

@test_show MIME"text/plain" begin
    "  leaf <args> [options] [flags]" in leaf
    "Args\n\n" in leaf
    "  <arg>" in leaf
    "Flags\n\n" in leaf
end

@test_throws ErrorException NodeCommand(;name="abc", subcmds=Dict{String, Any}())

node = NodeCommand(;name="foo", subcmds=Dict("leaf"=>leaf))

@test_show MIME"text/plain" begin
    "  foo <command>" in node
    "Commands\n\n" in node
    "  leaf <arg>" in node
    "Flags\n\n" in node
    "-h, --help" in node
end


text:
using CurveProximityQueries
using ConvexBodyProximityQueries
using Test
using LinearAlgebra
using StaticArrays
using IntervalArithmetic
using Unitful: nm, μm, mm, cm

import CurveProximityQueries: differentiate, integrate

@testset "CurveProximityQueries" begin
    @testset "Constructors" begin
        # should return SVector when constructed without SVectors
        @test Bernstein([[0.0, 0.0], [1.0, 1.0]])(0.5) === @SVector([0.5, 0.5])

        # non-floats arguments should become floats to satisfy assumptions
        # elsewhere in the code (e.g. use of `eps(T)`)
        @test eltype(Bernstein([[0, 0], [1, 1]])) <: Float64

        # none of the following should error; some should be inferrable
        @test Bernstein([
                @SVector([0.0, 0.0]),
                @SVector([1.0, 1.0])
            ])(0.5) === @SVector([0.5, 0.5])
        @test @inferred(Bernstein(@SVector([
                @SVector([0.0, 0.0]),
                @SVector([1.0, 1.0])
            ])))(0.5) === @SVector([0.5, 0.5])
        @test @inferred(Bernstein(
                @SVector([0.0, 0.0]),
                @SVector([1.0, 1.0])
            ))(0.5) === @SVector([0.5, 0.5])

        # similar reasoning below, but allow for any StaticVector
        struct Point{T} <: StaticArrays.FieldVector{2, T}
            x::T
            y::T
        end
        StaticArrays.similar_type(::Type{P}, ::Type{T},
            ::StaticArrays.Size{(2,)}) where {P <: Point, T} = Point{T}

        @test Bernstein([Point(0.0, 0.0), Point(1.0, 1.0)])(0.5) ===
            Point(0.5, 0.5)
        @test @inferred(
            Bernstein(@SVector([Point(0.0, 0.0), Point(1.0, 1.0)])))(0.5) ===
                Point(0.5, 0.5)
        @test @inferred(Bernstein(Point(0.0, 0.0), Point(1.0, 1.0)))(0.5) ===
            Point(0.5, 0.5)
    end
    @testset "Bernstein Polynomials" begin
        B2 = rand(Bernstein{2,8})
        B3 = rand(Bernstein{3,5})
        @test typeof(B2) <: Bernstein{2,8}
        @test typeof(B3) <: Bernstein{3,5}
        @test sum(norm.((differentiate(integrate(B2.f)) - B2.f).control_points)) ≤ 1e-10
        @test sum(norm.((differentiate(integrate(B3.f)) - B3.f).control_points)) ≤ 1e-10
    end
    @testset "Curve - Polygon" begin
        obs = @point zeros(2)
        c = Bernstein([ [0.23603334566204692, 0.34651701419196046],
                        [0.3127069683360675, 0.00790928339056074],
                        [0.4886128300795012, 0.21096820215853596],
                        [0.951916339835734, 0.9999046588986136],
                        [0.25166218303197185, 0.9866663668987996],
                        [0.5557510873245723, 0.43710797460962514],
                        [0.42471785049513144, 0.773223048457377] ])
        @test abs(minimum_distance(obs, c) - 0.3683522584741768) ≤ 1e-5
        @test tolerance_verification(obs, c, 0.3) == true
        @test tolerance_verification(obs, c, 0.4) == false
        @test collision_detection(obs, c) == false
        @test abs(minimum_distance(c, obs) - 0.3683522584741768) ≤ 1e-5
        @test tolerance_verification(c, obs, 0.3) == true
        @test tolerance_verification(c, obs, 0.4) == false
        @test collision_detection(c, obs) == false

        obs = @point zeros(3)
        c = Bernstein([ [0.23603334566204692, 0.34651701419196046, 0.3127069683360675],
                        [0.00790928339056074, 0.4886128300795012, 0.21096820215853596],
                        [0.951916339835734, 0.9999046588986136, 0.25166218303197185],
                        [0.9866663668987996, 0.5557510873245723, 0.43710797460962514],
                        [0.42471785049513144, 0.773223048457377, 0.2811902322857298],
                        [0.20947237319807077, 0.25137920979222494, 0.02037486871266725],
                        [0.2877015122756894, 0.859512136087661, 0.07695088688120899],
                        [0.6403962459899388, 0.8735441302706854, 0.27858242002877853],
                        [0.7513126327861701, 0.6448833539420931, 0.07782644396003469],
                        [0.8481854810000327, 0.0856351682044918, 0.5532055454580578],
                        [0.46335024592359875, 0.18582130997265378, 0.11198087695816716] ])
        @test abs(minimum_distance(obs, c) - 0.511627527056288) ≤ 1e-5
        @test tolerance_verification(obs, c, 0.5) == true
        @test tolerance_verification(obs, c, 0.6) == false
        @test collision_detection(obs, c) == false
        @test abs(minimum_distance(c, obs) - 0.511627527056288) ≤ 1e-5
        @test tolerance_verification(c, obs, 0.5) == true
        @test tolerance_verification(c, obs, 0.6) == false
        @test collision_detection(c, obs) == false
    end
    @testset "Curve - Curve" begin
        c = Bernstein([ [0.23603334566204692, 0.34651701419196046],
                        [0.3127069683360675, 0.00790928339056074],
                        [0.4886128300795012, 0.21096820215853596],
                        [0.951916339835734, 0.9999046588986136],
                        [0.25166218303197185, 0.9866663668987996],
                        [0.5557510873245723, 0.43710797460962514],
                        [0.42471785049513144, 0.773223048457377] ])
        d = Bernstein([ [0.2811902322857298, 0.20947237319807077],
                        [0.25137920979222494, 0.02037486871266725],
                        [0.2877015122756894, 0.859512136087661] ])
        @test abs(minimum_distance(c, d)) ≤ 1e-5
        @test tolerance_verification(c, d, 0.3) == false
        @test collision_detection(c, d) == true

        obs = @point zeros(3)
        c = Bernstein([ [0.23603334566204692, 0.34651701419196046, 0.3127069683360675],
                        [0.00790928339056074, 0.4886128300795012, 0.21096820215853596],
                        [0.951916339835734, 0.9999046588986136, 0.25166218303197185],
                        [0.9866663668987996, 0.5557510873245723, 0.43710797460962514],
                        [0.42471785049513144, 0.773223048457377, 0.2811902322857298],
                        [0.20947237319807077, 0.25137920979222494, 0.02037486871266725],
                        [0.2877015122756894, 0.859512136087661, 0.07695088688120899],
                        [0.6403962459899388, 0.8735441302706854, 0.27858242002877853],
                        [0.7513126327861701, 0.6448833539420931, 0.07782644396003469],
                        [0.8481854810000327, 0.0856351682044918, 0.5532055454580578],
                        [0.46335024592359875, 0.18582130997265378, 0.11198087695816716] ])
        d = Bernstein([ [0.976311881619359, 0.051614620674327094, 0.5380295812064833],
                        [0.4556920516275036, 0.2793951106725605, 0.17824610354168602],
                        [0.5489828719625274, 0.37097066286146885, 0.8941659192657594],
                        [0.6480537482231894, 0.4170393538841062, 0.14456554241360564],
                        [0.6224031828206811, 0.8723344353741975, 0.5249746566167794] ])
        @test abs(minimum_distance(c, d) - 0.137873546917387) ≤ 1e-5
        @test tolerance_verification(c, d, 0.1) == true
        @test tolerance_verification(c, d, 0.3) == false
        @test collision_detection(c, d) == false
    end
    @testset "Dimensionful curves" begin
        b2 = Bernstein([[0.0, 1.0]μm, [1.0, 1.0]μm, [1.0, 0.0]μm])
        @test b2(0.0) == [0, 1000]nm
        @test b2(1.0) == [0.001, 0]mm
        @test arclength(b2) == arclength(Bernstein([[0,1], [1,1], [1,0]])) * μm
        @test length(CurveProximityQueries.cvxhull(b2, Interval(0.0, 1.0))) == 8

        b3 = Bernstein([[0.0, 0.0, 0.0]mm, [1.0, 0.0, 0.0]mm, [1.0, 0.0, 1.0]mm])
        c3 = Bernstein([[2.0, 1.0, 1.0]cm, [2.0, 3.0, 2.0]cm, [3.0, 1.0, 1.0]cm])
        @test length(CurveProximityQueries.cvxhull(b3, Interval(0.0, 1.0))) == 16
        @test abs(minimum_distance(b3, c3, atol=1e-8mm) - 2.328cm) ≤ 1e-4cm
        @test tolerance_verification(b3, c3, 2cm, atol=1e-8mm) == true
        @test collision_detection(b3, c3, atol=1e-8mm) == false
    end
    @testset "Support" begin
        @test CurveProximityQueries.support(Point(1.0, 2.0), Point(1.0, 1.0)) ===
            Point(1.0, 2.0)
        @test CurveProximityQueries.support(@SVector(
            [Point(0.0, 1.0), Point(1.0, 2.0)]), @SVector([1.0, 1.0])) ===
                Point(1.0, 2.0)
        @test CurveProximityQueries.support(@SVector(
            [@SVector([0.0, 1.0]), @SVector([1.0, 2.0])]),
                @SVector([1.0, 1.0])) === @SVector([1.0, 2.0])
        @test CurveProximityQueries.support(@SVector(
            [Point(0.0, 1.0), Point(1.0, 2.0), Point(3.0, 4.0)]), Point(1.0, 1.0)) ===
                Point(3.0, 4.0)

    end
end

max_stars_repo_name: arlk/CurveProximityQueries

max_issues_repo_name: arlk/CurveProximityQueries

max_forks_repo_name: arlk/CurveProximityQueries


text:
#=##############################################################################
# DESCRIPTION
    Frequency-domain processing tools

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Sept 2020
  * License   : MIT
=###############################################################################

"""
    `autospec(x, fs, ns, overlap; window=hanning, density=false)`

Calculates the autospectrum of the waveform `x` using the periodigram method.

# Arguments
* `x::Array`       : (Pa) waveform.
* `fs::Real`       : (Hz) sample rate.
* `ns::Int`        : Size of blocks.
* `overlap::Real`  : Overlap of blocks between 0 and 1.
* `window::Function`   : Windowing function.
* `density::Bool`  : Returns spectrual density if true.
"""
function autospec(x, fs::Real, ns::Int, overlap::Real; window=hanning,
                                                    density=false, pref=2e-5)

    N = ns                     # Number of samples going through the FFT
    M = floor(Int, N/2+1)      # Frequencies to keep
    Gxx = zeros(M)             # Build the mean square of each frequency
    nave = 0

    delta_freq = fs/(N-1)
    freqs = [(m-1)*delta_freq for m in 1:M]

    # Iterate over blocks
    for firsti in range(1, length(x), step=ceil(Int, ns*(1-overlap)))
        endi = firsti + ns - 1

        if endi > length(x); break; end;

        # Generate window factors
        ww = window.(1:ns, ns)
        W = mean(ww.*conj(ww))  # Used to scale the ASD for energy conservation

        # Correction factor for single-sided fft (this is a copy/paste from Dr. Gee's function)
        # aux = sqrt(2*delta_freq/N/fs/W)
        aux = sqrt(2)*(1/N)*sqrt(1/W)

        # Grab the block
        this_x = x[firsti:endi]

        # Enforce zero mean
        # this_x .-= mean(this_x)

        # Window the block
        this_x = ww.*this_x

        # Corrected single-sided FFT the block
        fftx = aux*FFTW.fft(this_x)[1:M]

        # Magnitude square of FFT
#         fftx2 = (norm.(fftx)).^2
#         fftx2 = (abs.(fftx)).^2
        fftx2 = conj(fftx).*fftx

        # Build the average
        Gxx .+= fftx2
        nave += 1
    end

    Gxx /= nave

    OASPL = 20*log10(sqrt(sum(Gxx))/pref)

    # Normalize by the frequency bin to get density
    if density
        Gxx /= delta_freq
    end

    return Gxx, freqs, OASPL
end

"""
    `crossspec(x, fs, ns, overlap; window=hanning, density=false)`

Calculates the cross-spectrum of the waveforms `x` and `y` using the periodigram
method. This function is adapted from the MATLAB function by Kent Gee and Alan
Wall:

# Arguments
* `x::Array`       : (Pa) waveform.
* `fs::Real`       : (Hz) sample rate.
* `ns::Int`        : Size of blocks.
* `overlap::Real`  : Overlap of blocks between 0 and 1.
* `window::Function`   : Windowing function.
* `density::Bool`  : Returns spectrual density if true.
"""
function crossspec(x, y, fs::Real, ns::Int, overlap::Real; window=hanning,
                                                    density=true, pref=2e-5)

    # Clip data sets
    # N = Int(2^floor(log2( min(length(x), length(y)) )))
    N = min(length(x), length(y))
    x1 = x[1:N]
    x2 = y[1:N]

    # Enforce zero mean
    x1 .-= mean(x1)
    x2 .-= mean(x2)

    # Generate window factors
    ww = window.(1:ns, ns)
    W = mean(ww.*conj(ww))  # Used to scale the ASD for energy conservation

    # SPLITS DATA INTO BLOCKS
    #  Divides total data set into blocks of length ns with overlap, and
    #  windowed.
    N = ns                     # Number of samples going through the FFT
    M = floor(Int, N/2+1)      # Frequencies to keep
    Gxy = zeros(Complex, M)             # Build the mean square of each frequency
    nave = 0

    # Correction factor for single-sided fft (this is a copy/paste from Dr. Gee's function)
    scale = 2/ns/fs/W

    delta_freq = fs/(N-1)
    freqs = [(m-1)*delta_freq for m in 1:M]

    # Iterate over blocks
    for firsti in range(1, length(x1), step=ceil(Int, ns*(1-overlap)))
        endi = firsti + ns - 1

        if endi > length(x1); break; end;

        # Grab the block and window it
        this_x1 = ww.*x1[firsti:endi]
        this_x2 = ww.*x2[firsti:endi]

        # Corrected single-sided FFT the block
        fftx1 = FFTW.fft(this_x1)[1:M]
        fftx2 = FFTW.fft(this_x2)[1:M]

        # Build the average
        Gxy .+= scale*conj(fftx1).*fftx2     # Units are Pa^2/Hz
        nave += 1
    end

    Gxy /= nave

    aux = sqrt(2)*(1/N)*sqrt(1/W)
    OASPL = 20*log10(sqrt( aux^2/scale *  sum(Gxy))/pref)

    if density==false # Multiple by the frequency bin
        Gxy *= delta_freq
    end

    return Gxy, freqs, OASPL
end

"""
    `hanning(n, N)`
Hann window of period N sampled at n
"""
hanning(n::Int, N::Int) = n>N || n<1 ? 0 : sin(pi*(n-1)/(N-1))^2

triangular(n::Int, N::Int) = n>N || n<1 ? 0 : 1 - abs((n - N/2) / (N/2))

sinw(n::Int, N::Int) = n>N || n<1 ? 0 : sin(pi*n/N)

nuttall(n::Int, N::Int) = n>N || n<1 ? 0 :
                            0.355768 - 0.487396*cos(2*pi*n/N) +
                            0.144232*cos(4*pi*n/N) - 0.012604*cos(6*pi*n/N)

gaussian(n::Int, N::Int) = n>N || n<1 ? 0 : exp(-1/2*( (n-N/2) / (0.5*N/2) )^2)

expw(n::Int, N::Int) = n>N || n<1 ? 0 : exp(-abs(n-N/2)* 1/((N/2)/(60/8.69)))

"""

# Arguments
* `ps::Array{Float64, 1}`     : (Pa) pressure waveform
* `fs::Int`                   : (Hz) sample rate
"""
function spl_spectrum_from_autospectrum(ps, fs; pref=20e-6, numblocks=10,
                                                    overlap=0.5, window=hanning)

    ns = ceil(Int, length(ps)/numblocks)                   # Block size

    Gxx, freqs, OASPL = autospec(ps, fs, ns, overlap; window=window,
                                                        density=true, pref=pref)

    ww = window.(1:ns, ns)
    wwmean = mean(ww)
    df = freqs[2]-freqs[1]

    # Convert to amplitude
    ampl = sqrt.(Gxx .* df ./ wwmean^2)

    # Convert to SPL
    spls = 20*log10.(ampl./pref)

    return freqs, spls, OASPL
end

"""

# Arguments
* `ps::Array{Float64, 1}`     : (Pa) pressure waveform
* `fs::Int`                   : (Hz) sample rate
"""
function spl_spectrum_from_crossspectrum(ps, fs; numblocks=10, optargs...)

    ns = ceil(Int, length(ps)/numblocks)                   # Block size

    return spl_spectrum_from_crossspectrum(ps, fs, ns; optargs...)

end

function spl_spectrum_from_crossspectrum(ps, fs, ns; optargs...)

    return spl_crossspectrum_from_Gxy(ps, ps, fs, ns; optargs...)
end

function spl_crossspectrum_from_Gxy(ps_x, ps_y, fs, ns; pref=20e-6,
                                                    overlap=0.5, window=hanning)

    Gxy, freqs, OASPL = crossspec(ps_x, ps_y, fs, ns, overlap; window=window,
                                                        density=true, pref=pref)

    ww = window.(1:ns, ns)
    wwmean = mean(ww)
    df = freqs[2]-freqs[1]

    # Convert to amplitude
    ampl = sqrt.(Gxy .* df ./ wwmean^2)

    # Convert to SPL
    spls = 20*log10.(ampl./pref)

    return freqs, spls, OASPL
end

"""

# Arguments
* `ps::Array{Float64, 1}`     : (Pa) pressure waveform
* `fs::Int`                   : (Hz) sample rate
"""
function spl_spectrum_from_fft(ps, fs; pref=20e-6)

    dt = 1/fs                      # (s) sample period

    T = (length(ps)-1)/fs          # (s) total sampling time
    N = length(ps)                 # Number of samples

    # Fourier Transform
    ps_fft = 1/N*FFTW.fft(ps)  # FFTW doesn't multiply by 1/N for whatever reason
    ps_fftshift = FFTW.fftshift(ps_fft)

    # Define double sided frequency array
    dfreq = 1/T                 # (Hz) frequency step
    freqmax = floor(N/2)*dfreq  # (Hz) for even N, this is the Nyquist rate.
                                #      for odd N, this is 0.5*df less than the Nyquist rate

    # The DFT is N-periodic. The N+1th sample (if you went there) is the
    # same as the first sample (the Zero frequency bin).

    if mod(N, 2)==0 # Case N is even
    #=
        For even N, fmax is the nyquist rate.
        For an even N, the first sample is the 0 frequency and the (N/2
        +1)th sample is the sample at the Nyquist rate (by whatever standard that
        defines the FFT, the Nyquist rate is placed at the end of the negative
        side of the double sided array, so the positive side only goes up one bin
        less than the nyquist.
        The Nth sample is then the sample right before the 0 frequency
        (or equivalently the sampling frequency) again. because of periodicity
        the Nth sample is the same as the first negative frequency, the N-1th
        sample is the same as the second negative frequency, and so forth.
    =#
        freqs_raw = -freqmax:dfreq:(freqmax-dfreq)

    else # Case N is odd
    #=
        For odd N fmax is not the Nyquist, its half a bin below the Nyquist
        for odd N the first sample is the zero frequency and the floor(N/2)th
        sample represents the frequcny half a bin below the Nyquist. The
        floor(N/2)+1 to Nth frequencies then, by periodicity, are the same as the
        -floor(N/2)th ascending to the sample right before the zero frequency
        again.
    =#
        freqs_raw = -freqmax:dfreq:freqmax

    end

    freqs_unshifted = (0:(N-1))*dfreq
    freqs_singlesided = 0:dfreq:freqmax

    # (multiply the complex values by sqrt(2), so that when it's
    #  squared it's the same as being multiplied by 2, because we
    #  lost half the energy by only going single sided)
    ps_fft_singlesided = sqrt(2)*ps_fft[1:(floor(Int, N/2)+1)];

    spls = 20*log10.(abs.(ps_fft_singlesided)/pref)

    OASPL_freq_singlesided = 10*log10(sum(abs.(ps_fft_singlesided).^2)./pref^2)

    return freqs_singlesided, spls, OASPL_freq_singlesided
end

# """
# Convert the spectral density Gxx from narrow band to one-third octave band.
# """
# function narrow2onethird_Gxx(freqs, Gxx)
#
#     df = freqs[2]-freqs[1]                                # (Hz) Bin width
#
#     fc = [1, 1.25, 1.6, 2.0, 2.5, 3.15, 4, 5, 6.3, 8]     # (Hz) band center frequencies
#     fc = vcat(fc, fc*10, fc*100, fc*1000, fc*10000, 1e5)
#
#     flow = [0.88, 1.13, 1.414, 1.76, 2.25, 2.825, 3.53, 4.4, 5.65, 7.07] # (Hz) lower bound of 1/3-oct bins
#     flow = vcat(flow, flow*10, flow*100, flow*1000, flow*10000, 0.88e5)
#
#     fhi = [1.13, 1.414, 1.76, 2.25, 2.825, 3.53, 4.4, 5.65, 7.07, 8.8] # (Hz) upper bound of 1/3-oct bins
#     fhi = vcat(fhi, fhi*10, fhi*100, fhi*1000, fhi*10000, 1.13e5)
#
#     freqs_oto = fc                                        # (Hz) One-third octave frequency bins
#     Gxx_oto = zeros(eltype(Gxx), length(freqs_oto))       # One-third octave autospectrum
#
#     for i in 1:length(freqs)
#
#         # Find oto bin
#         otoi = nothing
#
#         for j in 1:length(freqs_oto)
#             if flow[j] <= freqs[i] &&  freqs[i] < fhi[j]
#                 otoi = j
#                 break
#             end
#         end
#
#         if otoi != nothing
#             Gxx_oto[otoi] += Gxx[i]*df
#         end
#     end
#
#     return freqs_oto, Gxx_oto
# end

# """
# # Convert spectral SPL from narrow band to one-third octave band.
# # """
# function narrow2onethird_spl_DEPRE(freqs, spls; pref=20e-6)
#
#     # Convert SPL to amplitude
#     ampl = pref * 10 .^ (spls ./ 20)
#
#     # Convert amplitude to pseudo-density
#     df = freqs[2]-freqs[1]
#     Gxx = ampl.^2 ./ df
#
#     # Convert narrow band to one-third octave band
#     freqs_oto, Gxx_oto = narrow2onethird_Gxx(freqs, Gxx)
#
#     # 1/3 octave band bounds
#     fc = freqs_oto
#     flow = vcat(fc[2:end]/sqrt(2), fc[end-2]*sqrt(2))
#     fhi = vcat(fc[3]/sqrt(2), fc[1:end-1]*sqrt(2))
#
#     # Convert pseudo-density to amplitude
#     # dfs = vcat(freqs_oto[2]-freqs_oto[1], [freqs_oto[i+1] - freqs_oto[i] for i in 1:length(freqs_oto)-1])
#     # ampl_oto = sqrt.(Gxx_oto .* dfs)
#     # dfs = [freqs_oto[i+1] - freqs_oto[i] for i in 1:length(freqs_oto)-1]
#     # ampl_oto = sqrt.(Gxx_oto[1:end-1] .* dfs)
#
#     dfs = fhi .- flow
#     ampl_oto = sqrt.(Gxx_oto .* dfs)
#
#     # Convert amplitude to SPL
#     spls_oto = 20*log10.(ampl_oto./pref)
#
#     # return freqs_oto[1:end-1], spls_oto
#     return freqs_oto, spls_oto
# end

"""
Convert spectral SPL from narrow band to one-third octave band.
"""
function narrow2onethird_spl(freqs, spls)

    spls_oto = zeros(eltype(spls), length(fc_oto)) # One-third octave autospectrum

    # Convert SPLs to pressure
    # and distribute all narrow band bins onto OTO bins
    otoi = 1
    for (i, freq) in enumerate(freqs) # Iterate over narrow bins

        # Find oto bin
        if freq < flow_oto[1]
            otoi = 1
        elseif freq >= fhi_oto[end]
            otoi = length(fc_oto)
        else
            for j in otoi:length(fc_oto)
                if flow_oto[j] <= freq &&  freq < fhi_oto[j]
                    otoi = j
                    break
                end
                if j==length(fc_oto)
                    error("LOGIC ERROR: OTO bin not found!")
                end
            end
        end

        # println("$i\t$otoi\t$(freq)\t$(fc_oto[otoi])")

        # Convert SPL to pressure and add to this OTO bin
        spls_oto[otoi] += 10^(spls[i]/10)
    end

    # Convert OTO pressure to OTO SPL
    for i in 1:length(spls_oto)
        if spls_oto[i] == 0
            spls_oto[i] = -55
        else
            spls_oto[i] = 10*log10(spls_oto[i])
        end
    end

    return deepcopy(fc_oto), spls_oto
end

# """
# Convert spectral SPL from narrow band to one-third octave band.
# """
# function narrow2onethird_spl_DEPRE(freqs, spls; pref=20e-6)
#
#     # Convert SPL to amplitude
#     ampl = pref * 10 .^ (spls ./ 20)
#
#     # Convert amplitude to pseudo-density
#     df = freqs[2]-freqs[1]
#     Gxx = ampl.^2 ./ df
#
#     # Convert narrow band to one-third octave band
#     freqs_oto, Gxx_oto = narrow2onethird_Gxx(freqs, Gxx)
#
#     # 1/3 octave band bounds
#     fc = freqs_oto
#     flow = vcat(fc[2:end]/sqrt(2), fc[end-2]*sqrt(2))
#     fhi = vcat(fc[3]/sqrt(2), fc[1:end-1]*sqrt(2))
#
#     # Convert pseudo-density to amplitude
#     # dfs = vcat(freqs_oto[2]-freqs_oto[1], [freqs_oto[i+1] - freqs_oto[i] for i in 1:length(freqs_oto)-1])
#     # ampl_oto = sqrt.(Gxx_oto .* dfs)
#     # dfs = [freqs_oto[i+1] - freqs_oto[i] for i in 1:length(freqs_oto)-1]
#     # ampl_oto = sqrt.(Gxx_oto[1:end-1] .* dfs)
#
#     dfs = fhi .- flow
#     ampl_oto = sqrt.(Gxx_oto .* dfs)
#
#     # Convert amplitude to SPL
#     spls_oto = 20*log10.(ampl_oto./pref)
#
#     # return freqs_oto[1:end-1], spls_oto
#     return freqs_oto, spls_oto
# end

# """
# Convert spectral SPL from narrow band to one-third octave band.
# """
# function onethird2narrow_spl_DEPRE(freqs_oto, df, spls_oto; pref=20e-6)
#
#
#     # Convert SPL to OTO amplitude
#     ampl_oto = pref * 10 .^ (spls_oto ./ 20)
#
#     # 1/3 octave band bounds
#     fc = freqs_oto
#     flow = vcat(fc[2:end]/sqrt(2), fc[end-2]*sqrt(2))
#     fhi = vcat(fc[3]/sqrt(2), fc[1:end-1]*sqrt(2))
#
#     # Convert amplitude to pseudo-density
#     dfs_oto = fhi .- flow
#     Gxx_oto = ampl_oto.^2 ./ dfs_oto
#
#     # Convert one-third octave band to narrow band
#     freqs, Gxx = onethird2narrow_Gxx(freqs_oto, Gxx_oto, df)
#
#     # Convert pseudo-density to amplitude
#     ampl = sqrt.(Gxx .* df)
#
#     # Convert amplitude to SPL
#     spls = 20*log10.(ampl ./ pref)
#
#     return freqs, spls
# end
#
#
# """
# Convert the spectral density Gxx from one-third octave band to narrow band.
# """
# function onethird2narrow_Gxx(freqs_oto, Gxx_oto, df)
#
#     fc = freqs_oto
#     flow = vcat(fc[2:end]/sqrt(2), fc[end-2]*sqrt(2))
#     fhi = vcat(fc[3]/sqrt(2), fc[1:end-1]*sqrt(2))
#
#     freqs = collect(flow[1]:df:fhi[end])
#     Gxx = zeros(eltype(Gxx_oto), length(freqs))
#
#     for otoi in 1:length(fc)
#
#         # Calculate average spectral density in this interval
#         df_oto = fhi[otoi]-flow[otoi]
#         this_Gxx = Gxx_oto[otoi]/df_oto
#         # this_Gxx = Gxx_oto[otoi]
#         # this_Gxx = Gxx_oto[otoi]*df_oto/df
#         # this_Gxx = Gxx_oto[otoi]*df/df_oto
#
#         # Find corresponding narrow band interval
#         loi = findfirst(f -> f>=flow[otoi], freqs)
#         upi = findlast(f -> f<fhi[otoi], freqs)
#         upi = upi==nothing ? length(freqs) : upi
#
#         # Distribute average spectral density into the narrow band bins
#         Gxx[loi:upi] .= this_Gxx
#
#     end
#
#     return freqs, Gxx
# end

"""
Convert spectral SPL from narrow band to one-third octave band.
"""
function onethird2narrow_spl(freqs_oto, df, spls_oto; smooth=true)

    # Find interval of frequencies
    fc_min = freqs_oto[findfirst(spl -> spl>-55, spls_oto)]
    fc_max = freqs_oto[findlast(spl -> spl>-55, spls_oto)]
    flow_mini = findlast(f -> f<=fc_min, flow_oto)
    fhi_maxi = findfirst(f -> f>fc_max, fhi_oto)

    # # Chop the bins that are smaller than the narrow band frequency
    # # This leads to losing from small energy, but makes things simpler
    # while fhi_oto[flow_mini] - flow_oto[flow_mini] < df && flow_mini < length(fc_oto)
    #     flow_mini += 1
    # end

    flow_min = flow_oto[flow_mini]
    fhi_max = fhi_oto[fhi_maxi]

    freqs = collect(flow_min:df:fhi_max)
    spls = zeros(eltype(spls_oto), length(freqs))

    for otoi in 1:length(freqs_oto)

        # Find corresponding narrow band interval
        loi = findfirst(f -> f>=flow_oto[flow_mini + otoi], freqs)
        upi = findlast(f -> f<fhi_oto[flow_mini + otoi], freqs)
        loi = loi==nothing ? length(freqs) : loi
        upi = upi==nothing ? 1 : upi

        # println("$otoi\t$(flow_oto[otoi])\t$(fhi_oto[otoi])\t$loi\t$upi\t$(freqs[loi])\t$(freqs[upi])")
        # println("$(upi-loi+1)\t$(log10(upi-loi+1))\t$(spls_oto[otoi] - log10(upi-loi+1))")

        # Distribute this OTO spl bin into the narrow band SPL bins
        spls[loi:upi] .+= spls_oto[otoi] - 10*log10(upi-loi+1)

    end

    # Smooth curve up through spline if requested
    if smooth
        spls_oto_dist = zeros(eltype(spls_oto), length(spls_oto))

        # Find distributed OTO SPL values
        for (i, freq_oto) in enumerate(freqs_oto)

            # Range where this OTO freq could be
            rng = searchsorted(freqs, freq_oto)
            rng = length(rng) >= 1 ? rng : rng.stop:rng.start

            # Pick the first valid position in that range
            freq_i = rng[findfirst( i -> 1 <= i <= length(spls), rng)]

            spls_oto_dist[i] = spls[freq_i]
        end

        # Interpolate the distributed OTO SPL through spline
        spls[:] .= math.akima(freqs_oto, spls_oto_dist, freqs)
    end

    # Bring empty bins down to minimum dB
    for i in 1:length(spls)
        if spls[i] == 0 || spls[i] < -55
            spls[i] = -55
        end
    end

    return freqs, spls
end

max_stars_repo_name: byuflowlab/FLOWNoise

max_issues_repo_name: byuflowlab/FLOWNoise

max_forks_repo_name: byuflowlab/FLOWNoise


text:
import DataStructures
import SimpleTraits
import Random
import Base.push!

include("index_util.jl")

struct MinHash
    salts :: Vector{UInt}
    n_hashes :: Int
    """
    Construct minhash hash functions.
    # Arguments
    - `threshold` the jaccard similarity above which documents are probably returned to be similar.
    """
    function MinHash(max_n_hashes=200, rng=Random.MersenneTwister())
        new([rand(rng, UInt) for _ in 1:max_n_hashes], max_n_hashes)
    end
end

SimpleTraits.@traitfn function fingerprint(mh:: MinHash, shngls::::SimpleTraits.BaseTraits.IsIterator)
    [minimum((hash(s, salt) for s in shngls)) for salt in mh.salts]
end

function fingerprint_all(mh::MinHash, shingles)
    n_hashes = length(mh.salts)
    n_docs = length(shingles)
    all_fingerprints = Array{UInt, 2}(undef,n_hashes,n_docs)
    Threads.@threads for i in 1:n_docs
        all_fingerprints[:,i] = fingerprint(mh, shingles[i])
    end
    [all_fingerprints[:,i] for i in 1:n_docs]
end

"""
    Compute MinHash fingerprint for string using the `shingle` function.
"""
function fingerprint(mh:: MinHash, str::AbstractString; shingle_size=3)::Vector{UInt}
    fingerprint(mh, shingle(str; size=shingle_size))
end

"""
    Cut a given string into shingles of size `size`.
    ```jldoctest
    julia> shingle("abcd"; size = 2)
    ["ab","bc","cd"]
    ```
"""
function shingle(s::AbstractString; size = 3)::Vector{String}
    if length(s) <= size
        return [s]
    end
    shingles = Vector{String}()
    for i in 1:(length(s) - size + 1)
        push!(shingles, s[thisind(s, i): thisind(s, i + size - 1)])
    end
    unique(shingles)
end

function estimate_jaccard(a::Vector{UInt}, b::Vector{UInt})
    sum(a .== b) / length(a)
end

mutable struct MinHashIndex <: LSHIndex
    threshold :: Float64
    tables :: Vector{DataStructures.DefaultDict{Vector{UInt}, Vector{Int}}}
    bands :: Int
    rows :: Int
    current_index :: Int
    max_n_hashes :: Int
    function MinHashIndex(;fingerprinter = MinHash(), threshold=0.5)
        max_n_hashes = fingerprinter.n_hashes
        bands, rows = _get_partition(threshold, max_n_hashes)
        tables = [DataStructures.DefaultDict{Vector{UInt}, Vector{Int}}(Vector{Int}) for _ in 1:bands]
        new{Vector{UInt}}(threshold, tables, bands, rows, 1, max_n_hashes)
    end
end

"""
Push new fingerprint to index.
"""
function push!(mhind:: MinHashIndex, fingerprint::Vector{Int}) where {FingerprintType, HashVectorType}
    index = mhind.current_index
    hashes = _bands(mhind, fingerprint)
    for (fingerprint_part, table) in zip(hashes, mhind.tables)
        push!(table[fingerprint_part], index)
    end
    mhind.current_index += 1;
end

function _bands(mhind, fingerprint::Vector{Int})
    (fingerprint[(i - 1) * mhind.rows + 1 : min(i * mhind.rows, mhind.max_n_hashes)] for i in 1:length(mhind.tables))
end

"""
Find probably similar indices in index.
"""
function find_similar(mhind:: MinHashIndex{HashVectorType},
    fingerprint::HashVectorType)::Vector{Int} where {FingerprintType, HashVectorType}
    indices = Vector{Vector{Int}}()
    bands = _bands(mhind, fingerprint)
    for (band, table) in zip(bands, mhind.tables)
        push!(indices, table[band])
    end
    unique(vcat(indices...))
end


text:
using ImageFiltering

img = rand(50,50)
alg = Algorithm.FIR
factor = ImageFiltering.centered(ones(13))
factored_ker = kernelfactors((factor, factor))

# works
imfilter(img, factored_ker, "symmetric")

# breaks
imfilter(img, factored_ker, "symmetric", alg)


text:
const OperatorVector = Vector{Operator{GenericBasis{Vector{Int64}}, GenericBasis{Vector{Int64}}, Matrix{ComplexF64}}}
const OperatorVectorArray = Array{ComplexF64, 3}

function operator_recast(rho::OperatorVector)::Array
    N = length(rho)
    M, K = size(rho[1].data)
    rho_array = zeros(eltype(rho[1].data), N, M, K)
    for t_i in 1:N
        rho_array[t_i, :, :] = rho[t_i].data
    end
    return rho_array
end

function operator_recast(rho::OperatorVectorArray)::Array
    return rho
end

######

function interaction_pic_to_schroedinger_pic(rho_int::OperatorVectorArray, tspan::Array, agg::Aggregate)
    rho_sch = deepcopy(rho_int)
    Ham_sys = agg.operators.Ham_sys
    for t_i in 1:length(tspan)
        t = tspan[t_i]
        U_op = evolutionOperator(Ham_sys, t)
        rho_sch[t_i, :, :] = U_op.data * rho_int[t_i, :, :] * (U_op').data
    end
    return rho_sch
end

function interaction_pic_to_schroedinger_pic(rho_int::OperatorVector, tspan::Array, agg::Aggregate)
    rho_array = operator_recast(rho_int)
    return interaction_pic_to_schroedinger_pic(rho_array, tspan, agg)
end

function schroedinger_pic_to_interaction_pic(rho_sch::OperatorVectorArray, tspan::Array, agg::Aggregate)
    rho_int = deepcopy(rho_sch)
    Ham_sys = agg.operators.Ham_sys
    for t_i in 1:length(tspan)
        t = tspan[t_i]
        U_op = evolutionOperator(Ham_sys, t)
        rho_int[t_i, :, :] = (U_op').data * rho_sch[t_i, :, :] * U_op.data
    end
    return rho_int
end

function schroedinger_pic_to_interaction_pic(rho_sch::OperatorVector, tspan::Array, agg::Aggregate)
    rho_array = operator_recast(rho_sch)
    return schroedinger_pic_to_interaction_pic(rho_array, tspan, agg)
end

#####

function local_st_to_exciton_st(rho_local::OperatorVectorArray, agg::Aggregate)
    N = size(rho_local, 1)
    rho_exciton = deepcopy(rho_local)
    Ham_sys = agg.operators.Ham_sys
    Ham_sys_lambda, Ham_sys_S = eigen(Ham_sys.data)
    Ham_sys_Sinv = inv(Ham_sys_S)
    # Ham_sys_lambda = diagm(Ham_sys_lambda)

    for t_i in 1:N
        rho_exciton[t_i, :, :] = Ham_sys_Sinv * rho_local[t_i, :, :] * Ham_sys_S
    end
    return rho_exciton
end

function local_st_to_exciton_st(rho_local::OperatorVector, agg::Aggregate)
    rho_array = operator_recast(rho_local)
    return local_st_to_exciton_st(rho_array, agg)
end

function exciton_st_to_local_st(rho_exciton::OperatorVectorArray, agg::Aggregate)
    N = size(rho_exciton, 1)
    rho_local = deepcopy(rho_exciton)
    Ham_sys = agg.operators.Ham_sys
    Ham_sys_lambda, Ham_sys_S = eigen(Ham_sys.data)
    Ham_sys_Sinv = inv(Ham_sys_S)
    # Ham_sys_lambda = diagm(Ham_sys_lambda)

    for t_i in 1:N
        rho_local[t_i, :, :] = Ham_sys_S * rho_exciton[t_i, :, :] * Ham_sys_Sinv
    end
    return rho_local
end

function exciton_st_to_local_st(rho_exciton::OperatorVector, agg::Aggregate)
    rho_array = operator_recast(rho_exciton)
    return exciton_st_to_local_st(rho_array, agg)
end

###

function tspan_cm_to_fs(tspan_cm)
    c_const = 299792458.0
    pi_const = 3.141592653589793
    c = 2.0*pi_const*c_const*1e-13
    tspan_fs = map((x) -> x / c, tspan_cm)
    return tspan_fs
end

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems

max_forks_repo_name: detrin/OpenQuantumSystems


text:
using Test
using DelimitedFiles
using Statistics
using openBF

@testset "openBF.jl" begin

    #unit tests
    println("Test initialise.jl functions")
    include("test_initialise.jl")

    println("Test boundary_conditions.jl functions")
    include("test_boundary_conditions.jl")

    #integration tests
    println("Test networks")
    include("test_single-artery.jl")
    include("test_conjunction.jl")
    include("test_bifurcation.jl")
    include("test_aspirator.jl")
    include("test_tapering.jl")
end


text:
module DataGeneration
include("./einsum3.jl")
import ..Einsum.@einsumnum
include("./cluster_test_simods.jl")
import .ClusterTests
using LinearAlgebra
using Random
using Distributions
using LightGraphs

"""
Inputs:
mu: d-by-k, each column j is the coordinate of centroid for Gaussian j;
sigma: 1-by-k, each entry j is the standard deviation of Gaussian j (we
assume the covariance matrix of component j is sigma_j * identity matrix;
w: 1-by-k, mixture probabilities;
n: number of samples
Outputs:
a: d-by-n, data of SON;
C_theta: n-by-1, the value of entry j is the assignment of j to its
Gaussian component i, the value is 0 if point j is > 0.5 * theta * min pairwise distance
among the means, away from its mean;
Vm: the list of indices of points which are within theta standard
deviation from its component mean
"""

function  gm_generator(mu, sigma, w, n, theta)

    rng = MersenneTwister(345)
    d, k = size(mu)
    compID = zeros(Int, n)
    mu_vec = zeros(d, n)
    sigma_vec = zeros(n)
    weight = zeros(n)

    x = rand(rng, n)
    for j = 1 : k
        if j == 1
            compID[x .<= w[1]] .= 1
        else
            compID[(x .> sum(w[1 : j - 1])) .&  (x .<= sum(w[1 : j]))] .= j
        end
        mu_vec[:, compID .== j] .= repeat(mu[:,j], 1, sum(compID .== j))
        sigma_vec[compID .== j] .= sigma[j]
    end
    y = randn(rng, (d, n))
    a = y .* repeat(sigma_vec', d, 1) .+ mu_vec
    Id = Matrix{Float64}(I,2,2)
    distr = MvNormal([0.0,0.0], Id)

    # a = repeat(a, 1, 2)
    # mu_vec = repeat(mu_vec, 1, 2)
    # sigma_vec = repeat(sigma_vec, 1, 2)
    # compID = repeat(compID, 2, 1)
    # weight = 50 .* pdf(distr, y)
    # weight = repeat(weight, 2, 1)
    weight = 100 .* pdf(distr, y)
    weight = weight[:]
    n = n

    C_theta = copy(compID)
    for j = 1: n
        if norm(a[:,j] - mu_vec[:, j]) > theta * sigma_vec[j]
            C_theta[j] = 0
        end
    end
    Vm = findall(C_theta .!= 0 )

    return a, compID, Vm, C_theta, weight
end

function twoGaussianspirals(N, std, noise)
    rng = MersenneTwister(345)
    N1 = floor(N/2)
    N2 = N-N1
    N1 = convert(Int64, N1)
    N2 = convert(Int64, N2)

    u1 = randn(rng, N1, 1)
    u2 = randn(rng, N2, 1)
    gamma = std

    u1 = randn(rng, N1, 1)
    d1 = [cos.(gamma .* u1) .+ rand(rng, N1, 1) .* noise sin.(gamma .* u1) .+ rand(rng, N1, 1) .* noise]

    u2 = randn(rng, N2, 1)
    d2 = [-cos.(gamma .* u2) .+ rand(rng, N2, 1) .* noise sin.(gamma .* u2) .+ rand(rng, N2, 1) .* noise .+ 1]

    a = transpose(vcat(d1, d2))
    compID = vcat(ones(Int, N1, 1), 2 .* ones(Int, N2, 1))

    return a, compID
end

function twohalfmoons(N, noise, std, amt)
    rng = MersenneTwister(345)
    N1 = floor(N/2)
    N2 = N-N1
    N1 = convert(Int64, N1)
    N2 = convert(Int64, N2)

    u1 = rand(rng, N1, 1)
    u2 = rand(rng, N2, 1)

    gamma = pi
    u1 = rand(rng, N1, 1) .- 0.5
    d1 = [cos.(gamma .* u1) .+ (rand(rng, N1, 1) .- 0.5) .* noise sin.(gamma .* u1) .+ (rand(rng, N1, 1) .- 0.5) .* noise]

    u2 = rand(rng, N2, 1) .- 0.5
    d2 = [-cos.(gamma .* u2) .+ (rand(rng, N2, 1) .- 0.5) .* noise sin.(gamma .* u2) .+ (rand(rng, N2, 1) .- 0.5) .* noise .+ 1]

    distr = Normal(0, std)
    f1 = amt .* pdf.(distr, gamma .* u1)
    f2 = amt .* pdf.(distr, gamma .* u2)

    a = transpose(vcat(d1, d2))
    weight = vcat(f1, f2)
    weight = weight[:]
    compID = vcat(ones(Int, N1, 1), 2 .* ones(Int, N2, 1))
    return a, compID, weight
end

function compute_randind(compID, C_admm)
    n = length(compID)

    inda = 0.0
    indb = 0.0
    indc = 0.0
    indd = 0.0

    inconclusive = findall(isequal(0), C_admm)
    num_inconclusive = length(inconclusive)
    wrongpair_inconclusive = num_inconclusive * (n - num_inconclusive) + num_inconclusive * (num_inconclusive-1) / 2

    conclusive = setdiff(1:n, inconclusive)
    num_conclusive = n - num_inconclusive

    j = 0
    k = 0
    for p = 1:num_conclusive
        for q = (j+1):num_conclusive
            j = conclusive[p]
            k = conclusive[q]
            Sj = compID[j]
            Sk = compID[k]
            Xj = C_admm[j]
            Xk = C_admm[k]

            if Sj == Sk
                if Xj == Xk
                    inda = inda + 1.0
                else
                    indc = indc + 1.0
                end
            else
                if Xj == Xk
                    indd = indd + 1.0
                else
                    indb = indb + 1.0
                end
            end
        end
    end
    Randind = (inda + indb)/(inda + indb + indc + indd + wrongpair_inconclusive);
    return Randind, inda, indb, indc, indd
end

function compute_weighted_randind(compID, C_admm, weight)
    n = length(compID)

    inda = 0.0
    indb = 0.0
    indc = 0.0
    indd = 0.0

    inconclusive = findall(isequal(0), C_admm)
    num_inconclusive = length(inconclusive)
    wrongpair_inconclusive = num_inconclusive * (n - num_inconclusive) + num_inconclusive * (num_inconclusive-1) / 2

    conclusive = setdiff(1:n, inconclusive)
    num_conclusive = n - num_inconclusive

    j = 0
    k = 0

    total_weight = 0.0
    for p = 1:n
        for q = (p+1):n
            total_weight = total_weight + weight[p] * weight[q]
        end
    end

    for p = 1:num_conclusive
        for q = (j+1):num_conclusive
            j = conclusive[p]
            k = conclusive[q]
            Sj = compID[j]
            Sk = compID[k]
            Xj = C_admm[j]
            Xk = C_admm[k]
            weightjk = weight[j] * weight[k]

            if Sj == Sk
                if Xj == Xk
                    inda = inda + weightjk
                else
                    indc = indc + weightjk
                end
            else
                if Xj == Xk
                    indd = indd + weightjk
                else
                    indb = indb + weightjk
                end
            end
        end
    end
    weighted_Randind = (inda + indb)/total_weight
    return weighted_Randind, inda, indb, indc, indd
end

function count_correctly_clustered_general(k, C_admm, compID)

    # For m=1:k, select the cluster among
    # from C_admm that has the
    # largest number of nodes from V_m
    # This is called mapclust (mapping from numbering
    # in C_theta to number in C_admm).
    mapclust = zeros(Int,k)
    n = length(compID)
    conclusive_index = findall(C_admm .!= 0)
    num_conclusive = length(conclusive_index)
    # println(num_conclusive)
    num_inconclusive = n - num_conclusive
    compID = compID[conclusive_index]
    C_admm = C_admm[conclusive_index]
    for m = 1 : k
        V_m = findall(compID .== m)
        u = zeros(Int, n)
        for i = 1 : length(V_m)
            ci = C_admm[V_m[i]]
            u[ci] = u[ci] + 1
        end
        scrap, maxoverlapm2 = findmax(u)
        mapclust[m] = maxoverlapm2
    end

    # Now score the clustering as follow:
    # Compute numvm_clust, numvm_total,
    # numall_clust, num_distinct_clust

    num_clust = 0
    for i = 1 : num_conclusive
        num_clust += (C_admm[i] == mapclust[compID[i]]) ? 1 : 0
    end
    num_distinct_clust = length(unique(mapclust))
    return num_clust, num_distinct_clust, num_inconclusive
end

"""
n = 80    # number of points
d = 2     # dimension
k = 2    # number of Gaussians
# sigma = 1.0    # standard deviation is sigma * sqrt(2)
theta = 0.7    # we only count Rand index for data points that are located at 0.5 * theta * min mean distance away from the respective means
lambdamulseq = 0.1:0.1:0.8    # the multiple of the lambda value which is supposed to recover Gaussians exactly from the theorem
meandistvec = 1:6    # the range of mean distance, mean distance is the ratio between the pairwise distance of Gaussian means and standard deviation
maxit = 10000    # maximum iteration
outfilename = "test_lambda3.txt"

"""
function run_test_lambda_rho_meandistance(n, d, k, theta, lambdamulseq, meandistvec, maxit, outfilename)
    open(outfilename, "w") do oh
        println(oh, "n = ", n, " d = ", d, " k = ", k, " theta = ", theta, "\n")
        sigma = sqrt(2)
        for meandist in meandistvec
            ee = Matrix{Float64}(meandist * I,d,d)
            mu = ee[1:k, :]
            # println(meandist)
            a, compID, Vm, C_theta, weight = gm_generator(mu,
                                                  sigma * ones(k),
                                                  ones(k) / k,
                                                  n,
                                                  theta)

            best_Randind_vec = Vector{Float64}()
            best_Randind_vec_s = Vector{Float64}()
            best_lambda_vec = Vector{Float64}()
            best_lambda_vec_s = Vector{Float64}()
            average_weight = sum(weight) / n
            theta_ineq = 0.7 * 0.5 * meandist
            chicum = cdf(Chi(d), theta_ineq)
            lambdamin = 2 * theta_ineq * sigma * k / (chicum * n * average_weight * 100) # the extra term 5 in the denominator is to ensure lambdamin < lambdamax
            # println(lambdamin)
            lambdamax = meandist * sqrt(2) / (2 * (n - 1) * average_weight * 0.1) # the extra term 0.2 in the denominator is to ensure lambdamin < lambdamax
            # println(lambdamax)
            println(oh, "------")
            println(oh, "lambdamin = ", lambdamin, " lambda_max = ", lambdamax)
            @assert lambdamin<lambdamax
            lambdamid = (lambdamin + lambdamax) / 2
            lambdavec = lambdamulseq * lambdamid
            Randind_vec = Vector{Float64}()
            Randind_vec_s = Vector{Float64}()
            num_point_s = 0
            for lambda in lambdavec
                if lambda <= 0
                    continue
                end
                X_s, Y_s, Z_s, s_s, u_s, t_s, Alpha_s, Beta_s, gamma_s, Optimality_cond, gap, clusterID, itcount,
                    clustering_cond, min_strict_comp = ClusterTests.admm_Chi_Lange_wtest(a, lambda, weight, maxit)
                num_clust, num_distinct_clust, num_inconclusive = count_correctly_clustered_general(k, clusterID, compID)
                Randind, inda, indb, indc, indd = compute_randind(compID, clusterID)
                Randind_s, inda_s, indb_s, indc_s, indd_s = compute_randind(compID[Vm], clusterID[Vm])
                push!(Randind_vec, Randind)
                push!(Randind_vec_s, Randind_s)
                num_point_s = length(Vm)
                println(oh, "lambda = ", lambda, " duality gap = ", gap,
                    " \n  num_clust = ", num_clust, " numvm_total = ",
                    n, " num_distinct_clust = ", num_distinct_clust, " num_inconclusive = ", num_inconclusive,
                    "\n clusters = ", clusterID)
                println(oh, "  itcount = ", itcount, " maxit = ", maxit)
                println(oh, "  clustering condition is ", clustering_cond)
                println(oh, "Rand index = ", Randind, ";")
                println(oh, "The smaller dataset has ", num_point_s, " points;")
                    println(oh, "Rand index for a smaller data set = ", Randind_s, ";")
                if length(unique(clusterID)) == 1 && num_inconclusive == 0
                    println(oh, "We have reached the max lambda!")
                    break
                end
            end
            bestrand_val, bestrand_location = findmax(Randind_vec)
            push!(best_Randind_vec, bestrand_val)
            push!(best_lambda_vec, lambdavec[bestrand_location])
            bestrand_val_s, bestrand_location_s = findmax(Randind_vec_s)
            push!(best_Randind_vec_s, bestrand_val_s)
            push!(best_lambda_vec_s, lambdavec[bestrand_location_s])
            println(oh, "The means are located ", meandist, " standard deviation away from each other. \n")
            println(oh, "The lambda values with best Rand index are ", best_lambda_vec,
                "\n The corresponding Rand indices are ", best_Randind_vec,
                "\n The lambda values with best Rand index of a smaller data set are ", best_lambda_vec_s,
                "\n The corresponding Rand indices of smaller data set are ", best_Randind_vec_s)
            println(oh, "The lambda values are ", lambdavec,
                "\n The Rand indices are ", Randind_vec,
                "\n The Rand idices for a smaller data set are ", Randind_vec_s)
            println(oh, "=======================================")
        end
    end
    nothing
end

"""
half moon dataset
"""

function run_test_lambda_halfmoon(n, std, noise, Gaussian, lambdamulseq, amtvec, outfilename, maxit)

    open(outfilename, "w") do oh

        k = 2
        numlambda = length(lambdamulseq)
        for amt in amtvec
            if Gaussian == true
                a, compID = twoGaussianspirals(n, std, noise)
            else
                a, compID, weight = twohalfmoons(n, noise, std, amt)
            end

            lambdamid = 0.05/amt
            lambdavec = lambdamulseq * lambdamid

            println(oh, "the data is:")
            println(oh, a)
            println(oh, "the weights are as follows:")
            println(oh, weight)
            println(oh, "------")
            println(oh, "n = ", n)
            # println(oh, "lambdamin = ", lambdamin, " lambda_max = ", lambdamax)
            println(oh, "the lambda vector is ", lambdavec)
            #weight = ones(n)
            #weight=weight_generator(a, rho, amt)
            mu_vec = zeros(numlambda)
            Randind_vec = zeros(numlambda)
            weighted_Randind_vec = zeros(numlambda)
            num_clust_vec = zeros(Int, numlambda,)
            num_distinct_clust_vec = zeros(Int, numlambda,)
            num_distinct_clust_vec = zeros(Int, numlambda,)
            itcount_all = zeros(Int, numlambda)
            clusteringcond_all = zeros(Bool, numlambda)
            min_strict_comp_all = zeros(numlambda)
            println(oh, "===============")
            println(oh, "The boosted amount is ", amt)
            for j = 1:numlambda
                lambda = lambdavec[j]
                X_s, Y_s, Z_s, s_s, u_s, t_s, Alpha_s, Beta_s, gamma_s, Optimality_cond, gap, clusterID, itcount,
                            clustering_cond, min_strict_comp = ClusterTests.admm_Chi_Lange_wtest(a, lambda, weight, maxit)
                num_clust, num_distinct_clust, num_inconclusive = count_correctly_clustered_general(k, clusterID, compID)
                Randind, inda, indb, indc, indd = compute_randind(compID, clusterID)
                weighted_Randind, weighted_inda, weighted_indb, weighted_indc, weighted_indd = compute_weighted_randind(compID, clusterID, weight)

                println(oh, "lambda = ", lambda, " duality gap = ", gap)
                println(oh, "Rand index = ", Randind, "weighted Rand index = ", weighted_Randind, ";")
                println(oh, "Same half moon, same cluster = ", inda, "; different half moons, different clusters = ", indb, ";")
                println(oh, "same half moon, different clusters = ", indc, "; different half moons, same cluster = ", indd)
                print(oh, " num_clust = ", num_clust, " numvm_total = ",
                        num_clust, " num_distinct_clust = ", num_distinct_clust,
                        "\n clustering = ", clusterID,
                        "\n clustering condition = ", clustering_cond)
                println(oh, "  itcount = ", itcount, " maxit = ", maxit)
                println(oh, "-------------------------------------")

                #println(oh, "X=", X_s)

                mu_vec[j] = gap
                Randind_vec[j] = Randind
                weighted_Randind_vec[j] = weighted_Randind
                num_clust_vec[j] = num_clust
                num_distinct_clust_vec[j] = num_distinct_clust
                itcount_all[j] = itcount
                clusteringcond_all[j] = clustering_cond
                min_strict_comp_all[j] = min_strict_comp

                if length(unique(clusterID)) == 1 && num_inconclusive == 0
                    println(oh, "We have reached the max lambda!")
                    break
                end
            end
            println(oh, "========================================")
            println(oh, "lambda = ", lambdavec, ";")
            println(oh, "Rand index = ", Randind_vec, ";")
            println(oh, "weighted Rand index = ", weighted_Randind_vec, ";")
            println(oh, "duality_gap_mu = ", mu_vec, ";")
            println(oh, "number_of_clustered_points = ", num_clust_vec, ";")
            println(oh, "number_of_distinct_clusters = ", num_distinct_clust_vec, ";")
            println(oh, "itcount_all = ", itcount_all, ";")
            println(oh, "clusteringcond_all = ",
                    [Int(clusteringcond_all[j]) for j = 1 : numlambda], ";")
            println(oh, "min_strict_comp_all = ", min_strict_comp_all, ";")
        end

    end
    nothing
end

end

max_stars_repo_name: taotolojiang/Certifying-clusters-SON

max_issues_repo_name: taotolojiang/Certifying-clusters-SON

max_forks_repo_name: taotolojiang/Certifying-clusters-SON


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
function add_existing_boiler_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
    r = Dict{String, Any}()

    # TODO all of these time series assume hourly time steps
    # TODO we convert MMBTU_TO_KWH from user inputs to the model, and then back to mmbtu in outputs: why not stay in mmbtu?
	r["year_one_boiler_fuel_consumption_mmbtu_per_hr"] =
        round.(value.(m[:dvFuelUsage]["ExistingBoiler", ts] for ts in p.time_steps) / MMBTU_TO_KWH, digits=3)
    r["year_one_boiler_fuel_consumption_mmbtu"] = round(sum(r["year_one_boiler_fuel_consumption_mmbtu_per_hr"]), digits=3)

	r["year_one_boiler_thermal_production_mmbtu_per_hr"] =
        round.(value.(m[:dvThermalProduction]["ExistingBoiler", ts] for ts in p.time_steps) / MMBTU_TO_KWH, digits=3)
	r["year_one_boiler_thermal_production_mmbtu"] = round(sum(r["year_one_boiler_thermal_production_mmbtu_per_hr"]), digits=3)

	# @expression(m, BoilerToHotTES[ts in p.time_steps],
	# 	m[:dvProductionToStorage]["HotTES","ExistingBoiler",ts])

	# r["boiler_thermal_to_tes_series"] = round.(value.(BoilerToHotTES), digits=3)

    # if !isempty(p.SteamTurbineTechs)
    #     @expression(m, BoilerToSteamTurbine[ts in p.time_steps], m[:dvThermalToSteamTurbine]["ExistingBoiler",ts])
    #     r["boiler_thermal_to_steamturbine_series"] = round.(value.(BoilerToSteamTurbine), digits=3)
    # else
    #     BoilerToSteamTurbine = zeros(p.time_stepsCount)
    #     r["boiler_thermal_to_steamturbine_series"] = round.(BoilerToSteamTurbine, digits=3)
    # end

	# @expression(m, BoilerToLoad[ts in p.time_steps],
	# 	m[:dvThermalProduction]["ExistingBoiler",ts] - BoilerToHotTES[ts] - BoilerToSteamTurbine[ts])
	# r["boiler_thermal_to_load_series"] = round.(value.(BoilerToLoad), digits=3)

	@expression(m, TotalBoilerFuelCharges,
		p.pwf_fuel["ExistingBoiler"] * p.hours_per_timestep * sum(
            p.s.existing_boiler.fuel_cost_series[ts] * m[:dvFuelUsage]["ExistingBoiler", ts] for ts in p.time_steps
        )
    )
	r["total_boiler_fuel_cost"] = round(value(TotalBoilerFuelCharges) * (1 - p.s.financial.offtaker_tax_pct), digits=3)
	r["year_one_boiler_fuel_cost"] = round(value(TotalBoilerFuelCharges) / p.pwf_fuel["ExistingBoiler"], digits=3)

    d["ExistingBoiler"] = r
	nothing
end

max_stars_repo_name: NREL/REoptLite

max_issues_repo_name: NREL/REoptLite

max_forks_repo_name: NREL/REoptLite


text:
using AdditionalUnits
using Test
using Unitful: Unitful, @u_str, J, W, hr

@testset "AdditionalUnits.jl" begin
    @testset "Power Related Unit" begin
        @test 1u"Wh" == 1u"W*hr"
    end
end


text:
# 加载补丁
include("lib/JSON2Ex.jl")

max_stars_repo_name: huaxk/rlservice

max_issues_repo_name: huaxk/rlservice

max_forks_repo_name: huaxk/rlservice


text:
"""
Callback to stop optimization if the function value is below some threshold.
"""
immutable FunctionStop{T}
    val::T
end

function (cb::FunctionStop)(os::OptimizationState)
    os.value ≤ cb.val
end

function (cb::FunctionStop)(tr::OptimizationTrace)
    tr[end].value ≤ cb.val
end

abstract type StopTolerance end;

"""
Tell the optimization procedure to exit when we hit `|g(x)|≤ tol * |g(x_0)|`.

For absolute, rather than relative, tolerance, set `g0norm = 1.0`.
"""
immutable GradientTolerance{T} <: StopTolerance
    # TODO: The g0norm is a bit superfluous
    tol::T
    g0norm::T
end

"""
Tell the optimization procedure to exit when we hit `f(x) ≤ f_L + tol * (f(x_0) - f_L)`,
where `f_L` is the minimum value of `f`.

For absolute, rather than relative, tolerance, set `f0 = 1.0 + fL`.
"""
immutable FunctionTolerance{T} <: StopTolerance
    # TODO: f0 and fL are superflous
    tol::T
    f0::T
    fL::T
end

"""
If we don't know the value, tell the optimization procedure that we want to find
the minimum function value achieved across a list of solvers,
by running them all to within a given gradient norm tolerance and number of iterations.
"""
immutable UnknownFunctionTolerance{T} <: StopTolerance
    # TODO: Maybe we just run these problems to high tolerance before
    #       starting this process, and then store the approximate minimum value.
    gtol::T
    iterations::Int
end

"Log of results from a given optimization run."
immutable OptimizationRun{T,Tf}
    iterations::Int
    fevals::Int
    gevals::Int
    hevals::Int
    cputime::T
    fval::Tf
    f0::Tf
    gnorm::Tf
    g0norm::Tf
    solvername::String
    problemname::String
    success::Bool
end

function OptimizationRun(iterations, fevals, gevals, hevals, cputime, fval,
                         f0, gnorm, g0norm, solvername, problemname,
                         ft::FunctionTolerance)
    success = fval ≤ ft.fL + ft.tol*(ft.f0 - ft.fL)
    OptimizationRun(iterations, fevals, gevals, hevals, cputime, fval, f0,
                    gnorm, g0norm, solvername, problemname, success)
end

function OptimizationRun(iterations, fevals, gevals, hevals, cputime, fval,
                         f0, gnorm, g0norm, solvername, problemname,
                         gt::GradientTolerance)
    success = gnorm ≤ gt.tol*gt.g0norm
    OptimizationRun(iterations, fevals, gevals, hevals, cputime, fval, f0,
                    gnorm, g0norm, solvername, problemname, success)
end

"Collection of options used create performance measures."
immutable TestSetup{T<:Real, Ts<:AbstractVector, Tn<:AbstractVector, Tl<:Real}
    solvers::Ts
    solvernames::Tn
    stoptype::Symbol
    stoptol::T
    timelog::Int # Number of calls to optimize to record timings. timelog = 0 calls optimize once, but does not record the timing
    maxiter::Int
    timelimit::Tl # Soft time limit for each optimize call
end


text:
import DSP.conv, RHEOS.convn, BenchmarkTools

# realistic test data sizes for RHEOS workflow
x = collect(0.0:0.1:1000.0)
y = exp.(-x/2.0)
z = [i<500.0 ? i : 500.0 - i for i in x]

# convolve them yeh

@btime conv($z, $y)
@btime convn($z, $y)

max_stars_repo_name: JuliaRheology/RheoBenchAndTest

max_issues_repo_name: JuliaRheology/RheoBenchAndTest

max_forks_repo_name: JuliaRheology/RheoBenchAndTest


text:
# interp1.jl

using MIRT: interp1

using Plots: scatter, plot!, gui
using Interpolations: Flat, Line, Periodic, Reflect, Throw
using Test: @test, @test_throws, @inferred

@test (@inferred interp1(1:4, 2:5, 1.5)) == 2.5
@test (@inferred interp1(1:2, [3im,4im], 1.5)) == 3.5im
x = LinRange(-1, 1, 21)
y = sin.(x*π)
xi = LinRange(-1, 1, 101) * 2
@test_throws BoundsError interp1(x, y, xi, extrap=nothing)
@test_throws BoundsError interp1(x, y, xi, extrap=Throw())
scatter(x, y, label="")

y0 = @inferred interp1(x, y, xi) # 0 extrapolation
plot!(xi, y0, label="0", ylim=[-2,2])
for ex in (Flat(), Line(), Periodic(), Reflect())
	yy = @inferred interp1(x, y, xi, extrap=ex)
	plot!(xi, yy, label="$ex")
end
gui()


text:
export max_a_posteriori, MAP, map_prob

import DataFrames: DataFrame, mapcols!

#####################
# Circuit MAP/MPE evaluation
#####################

"""
    max_a_posteriori(root::ProbCircuit, data::Union{Bool,Missing}...)
    max_a_posteriori(root::ProbCircuit, data::Union{Vector{<:Union{Bool,Missing}},CuVector{UInt8}})
    max_a_posteriori(circuit::ProbCircuit, data::DataFrame)
    max_a_posteriori(pbc::ParamBitCircuit, data; Float=Float32)

Evaluate maximum a-posteriori state of the circuit for given input(s).

Outputs the states, and the corresponding probabilities (in log domain).
"""
max_a_posteriori(root::ProbCircuit, data::Union{Bool,Missing}...) =
    max_a_posteriori(root, collect(Union{Bool,Missing}, data))

max_a_posteriori(root::ProbCircuit, data::Union{Vector{<:Union{Bool,Missing}},CuVector{UInt8}}) = begin
    map, prob = max_a_posteriori(root, DataFrame(reshape(data, 1, :), :auto))
    example(map, 1), prob[1]
end

max_a_posteriori(circuit::ProbCircuit, data::DataFrame; Float=Float32) =
    max_a_posteriori(same_device(ParamBitCircuit(circuit, data), data), data; Float)

function max_a_posteriori(pbc::ParamBitCircuit, data; Float=Float32)
    @assert isgpu(data) == isgpu(pbc) "ParamBitCircuit and data need to be on the same device"
    values = map_prob_all(pbc, data; Float)
    map_state =  map_down(pbc, data, values)
    @assert isgpu(values) == isgpu(map_state)
    map_state, values[:,end]
end

"""
Maximum a-posteriori queries
"""
const MAP = max_a_posteriori

"""
Mode of the distribution
"""
const mode = max_a_posteriori

#####################
# Circuit marginal evaluation
#####################

"""
    map_prob(root::ProbCircuit, data::Union{Real,Missing}...)
    map_prob(root::ProbCircuit, data::Union{Vector{Union{Bool,Missing}},CuVector{UInt8}})
    map_prob(circuit::ProbCircuit, data::DataFrame)
    map_prob(circuit::ParamBitCircuit, data::DataFrame)

The most likely world probability that agrees with the provided data"

Missing values should be denoted by `missing` in the data.
"""
map_prob(root::ProbCircuit, data::Union{Real,Missing}...) =
    map_prob(root, collect(Union{Bool,Missing}, data))

map_prob(root::ProbCircuit, data::Union{Vector{Union{Bool,Missing}},CuVector{UInt8}}) =
    map_prob(root, DataFrame(reshape(data, 1, :), :auto))[1]

map_prob(circuit::ProbCircuit, data::DataFrame) =
    map_prob(same_device(ParamBitCircuit(circuit, data), data) , data)

map_prob(circuit::ParamBitCircuit, data::DataFrame) =
    map_prob_all(circuit, data)[:,end]

#####################
# MAP Probability (upward pass)
#####################

"The most likely world that agrees with the provided data for all nodes"
map_prob_all(circuit::ProbCircuit, data::DataFrame) =
    map_prob_all(same_device(ParamBitCircuit(circuit, data), data), data)

function map_prob_all(circuit::ParamBitCircuit, data, reuse=nothing; Float=Float32)
    @assert num_features(data) == num_features(circuit)
    # here we can reuse `init_marginal` from the marginal computation because
    # our currently supported leafs all have MAP state `1`
    values = init_marginal(data, reuse, num_nodes(circuit); Float)
    map_prob_layers(circuit, values)
    return values
end

# upward pass helpers on CPU

"Compute MAP probabilities on the CPU (SIMD & multi-threaded)"
function map_prob_layers(circuit::ParamBitCircuit, values::Matrix)
    bc::BitCircuit = circuit.bitcircuit
    els = bc.elements
    pars = circuit.params
    for layer in bc.layers[2:end]
        Threads.@threads for dec_id in layer
            j = @inbounds bc.nodes[1,dec_id]
            els_end = @inbounds bc.nodes[2,dec_id]
            if j == els_end
                # first element value is just marginal value
                assign_marginal(values, dec_id, els[2,j], els[3,j], pars[j])
                j += 1
            else
                assign_map_prob(values, dec_id, els[2,j], els[3,j], els[2,j+1], els[3,j+1], pars[j], pars[j+1])
                j += 2
            end
            while j <= els_end
                accum_map_prob(values, dec_id, els[2,j], els[3,j], pars[j])
                j += 1
            end
        end
    end
end

assign_map_prob(v::Matrix{<:AbstractFloat}, i, e1p, e1s, e2p, e2s, p1, p2) = begin
    @avx for j=1:size(v,1)
        @inbounds x = v[j,e1p] + v[j,e1s] + p1
        @inbounds y = v[j,e2p] + v[j,e2s] + p2
        @inbounds v[j,i] = max(x, y)
    end
end

accum_map_prob(v::Matrix{<:AbstractFloat}, i, e1p, e1s, p1) = begin
    @avx for j=1:size(v,1)
        @inbounds x = v[j,i]
        @inbounds y = v[j,e1p] + v[j,e1s] + p1
        @inbounds v[j,i] = max(x, y)
    end
end

# upward pass helpers on GPU

"Compute marginals on the GPU"
function map_prob_layers(circuit::ParamBitCircuit, values::CuMatrix;  dec_per_thread = 8, log2_threads_per_block = 8)
    circuit = to_gpu(circuit)
    bc = circuit.bitcircuit
    CUDA.@sync for layer in bc.layers[2:end]
        num_examples = size(values, 1)
        num_decision_sets = length(layer)/dec_per_thread

        kernel = @cuda name="map_prob_layers_cuda" launch=false map_prob_layers_cuda(layer, bc.nodes, bc.elements, circuit.params, values)
        config = launch_configuration(kernel.fun)
        threads, blocks = balance_threads_2d(num_examples, num_decision_sets, config.threads)
        kernel(layer, bc.nodes, bc.elements, circuit.params, values
            ; threads, blocks)
    end
end

"CUDA kernel for circuit evaluation"
function map_prob_layers_cuda(layer, nodes, elements, params, values)
    index_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    index_y = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    stride_x = blockDim().x * gridDim().x
    stride_y = blockDim().y * gridDim().y
    for j = index_x:stride_x:size(values,1)
        for i = index_y:stride_y:length(layer)
            decision_id = @inbounds layer[i]
            k = @inbounds nodes[1,decision_id]
            els_end = @inbounds nodes[2,decision_id]
            @inbounds x = values[j, elements[2,k]] + values[j, elements[3,k]] + params[k]
            while k < els_end
                k += 1
                @inbounds y = values[j, elements[2,k]] + values[j, elements[3,k]] + params[k]
                x = max(x, y)
            end
            values[j, decision_id] = x
        end
    end
    return nothing
end

#####################
# MAP state (downward pass)
#####################

# helper function

"Find the MAP child value and node id of a given decision node"
function map_child(params, nodes, elements, ex_id, dec_id, values)
    @inbounds els_start = nodes[1,dec_id]
    @inbounds els_end = nodes[2,dec_id]
    pr_opt = typemin(eltype(values))
    j_opt = 1
    for j = els_start:els_end
        @inbounds prime = elements[2,j]
        @inbounds sub = elements[3,j]
        @inbounds pr = values[ex_id, prime] + values[ex_id, sub] + params[j]
        if pr > pr_opt
            pr_opt = pr
            j_opt = j
        end
    end
    @inbounds return elements[2,j_opt], elements[3,j_opt]
end

# CPU code

function map_down(pbc, data, values::Array)
    state = zeros(Bool, num_examples(data), num_features(data))
    Threads.@threads for ex_id = 1:size(state,1)
        map_rec(num_leafs(pbc), params(pbc), nodes(pbc), elements(pbc), ex_id, num_nodes(pbc), values, state)
    end
    df = DataFrame(state, :auto)
    mapcols!(c -> BitVector(c), df)
    return df
end

function map_rec(nl, params, nodes, elements, ex_id, dec_id, values, state)
    if isleafgate(nl, dec_id)
        if isliteralgate(nl, dec_id)
            l = literal(nl, dec_id)
            @inbounds state[ex_id, lit2var(l)] = (l > 0)
        end
    else
        prime, sub = map_child(params, nodes, elements, ex_id, dec_id, values)
        map_rec(nl, params, nodes, elements, ex_id, prime, values, state)
        map_rec(nl, params, nodes, elements, ex_id, sub,   values, state)
    end
end

# GPU code

function map_down(pbc, data, values::CuArray)
    state = CUDA.zeros(Bool, num_examples(data), num_features(data))
    stack = CUDA.zeros(Int32, num_examples(data), num_features(data)+3)
    @inbounds stack[:,1] .= 1 # start with 1 dec_id in the stack
    @inbounds stack[:,2] .= num_nodes(pbc) # start with the root in the stack

    CUDA.@sync begin
        kernel = @cuda name="map_cuda_kernel" launch=false map_cuda_kernel(num_leafs(pbc), params(pbc), nodes(pbc), elements(pbc), values, state, stack)
        config = launch_configuration(kernel.fun)
        threads = config.threads
        blocks = cld(size(state,1), threads)
        kernel(num_leafs(pbc), params(pbc), nodes(pbc), elements(pbc), values, state, stack
            ; threads, blocks)
    end
    # do the conversion to a CuBitVector on the CPU...
    df = DataFrame(to_cpu(state), :auto)
    mapcols!(c -> to_gpu(BitVector(c)), df)
    return df
end

function map_cuda_kernel(nl, params, nodes, elements, values, state, stack)
    index_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride_x = blockDim().x * gridDim().x
    for ex_id = index_x:stride_x:size(state,1)
        dec_id = pop_cuda!(stack, ex_id)
        while dec_id > zero(eltype(stack))
            if isleafgate(nl, dec_id)
                if isliteralgate(nl, dec_id)
                    l = literal(nl, dec_id)
                    var = lit2var(l)
                    @inbounds state[ex_id, var] = (l > 0)
                end
            else
                prime, sub = map_child(params, nodes, elements, ex_id, dec_id, values)
                push_cuda!(stack, prime, ex_id)
                push_cuda!(stack, sub, ex_id)
            end
            dec_id = pop_cuda!(stack, ex_id)
        end
    end
    return nothing
end


text:
using Base: @deprecate

@deprecate by{T,N}(ta::TimeArray{T,N}, t::Int; period::Function=day) when(ta, period, t)
@deprecate by{T,N}(ta::TimeArray{T,N}, t::String; period::Function=day) when(ta, period, t)

@deprecate to(ta::TimeArray, y::Int, m::Int, d::Int) to(ta, Date(y, m, d))
@deprecate from(ta::TimeArray, y::Int, m::Int, d::Int) from(ta, Date(y, m, d))

@deprecate findall(ta::TimeArray) find(ta)

@deprecate collapse{T,N,D}(ta::TimeArray{T,N,D}, timestamp::Function; period::Function=week) collapse(ta, period, timestamp)


text:
module BasicStats

    include("xlogx.jl")
    include("avg.jl")
    include("index.jl")
    export xlogx, avg, FisherIndex, normalization
end


text:
using Stars
using Test

function nan_equal(x, y)
    ind = findall(.!(isnan.(y)))
    x[ind] == y[ind]
end

@testset "Stars" begin
    cd(dirname(@__FILE__)) do
        # include("get_testdata.jl")
        # include("test_geoutils.jl")
        include("test-geoarray.jl")
        include("test-st_as_sf.jl")
        include("test-st_Ops.jl")

        include("test-gdalinfo.jl")
        include("test-gdal_polygonize.jl")

        include("test-readGDAL.jl")
        include("test-st_coords.jl")

        include("test-sp.jl")
        # include("test_interpolate.jl")
        # include("test_crs.jl")
        # include("test_operations.jl")
    end
end