text:
function dh_simulation(LP::AbstractLProblem{T}, p::Float64, nTrials::Int64, alg::OrdinaryDiffEqAlgorithm, logger::Logger{T1,T2} = Logger{:no, :cmd}(), interval::Int64 = 1, offset1::Int64 = 0; links = (:all,:all), axes = [:x, :y, :z], kwargs...) where {T<:PureClassicalApprox,T1,T2}
    assert(nTrials > 0)
    RNGarray = initRNG(nprocs())
    path_part = @sprintf("Dis Tm%.2f dt%.3f del%d Obs%s Pr%.3f/", LP.Tmax, LP.tstep, LP.delimiter, get_string(LP.rules.str_vec), p)
    filenames = Vector{String}()
    A = LP.A
    l = LP.cb.affect!.save_func.len
    tspan = (0.0, LP.Tmax*LP.delimiter)
    t_means = collect(0.0:LP.tstep:(LP.Tmax*LP.delimiter))
    saved_values = SavedValues(t_means, [zeros(l) for i in eachindex(t_means)])
    for thread = (offset1+1):(offset1+nprocs())
        push!(filenames, get_string(nTrials, thread, alg))
    end
    log_o, save_os = set_Logging(LP, nTrials, logger, A, path_part, filenames)
    cf_vals0 = LP.cf_vals
    @sync for i in workers()
        @spawnat i begin
                       global cf_vals = cf_vals0
                       global RNG = RNGarray[i]
                   end
    end
    iter = get_status(log_o)
    offset = get_offset(log_o)
    if iter > 0
        load_data!(cf_vals0, save_os[1])
        for i in workers()
            @spawnat i load_data!(cf_vals, save_os[i])
        end
    end
    t0 = time()-offset
    while iter<nTrials
        t1 = time()
        iter += 1
        ARHS = build_RHS_function(A)
        rules, OSET = set_CorrFuncs(A, links, axes)
        @sync begin
            for i in workers()
                @spawnat i dh_run(iter, A, ARHS, OSET, rules, t_means, saved_values, tspan, p, cf_vals, alg, save_os[i], RNG; kwargs...)
            end
            @spawnat 1 dh_run(iter, A, ARHS, OSET, rules, t_means, saved_values, tspan, p, cf_vals0, alg, save_os[1], RNGarray[1]; kwargs...)
        end
        t2 = time()
        output(log_o, @sprintf("%10d %20.5f %20.5f %10d", iter, t2-t1, (t2-t0)/iter, nTrials - iter))
    end
end

function dh_simulation(LP::AbstractLProblem{T}, p::Float64, nTrials::Int64, alg::OrdinaryDiffEqAlgorithm, logger::Logger{T1,T2} = Logger{:no, :cmd}(), interval::Int64 = 1, offset1::Int64 = 0; links = (:all,:allq), axes = [:x, :y, :z], kwargs...) where {T<:HybridApprox,T1,T2}
    assert(nTrials > 0)
    RNGarray = initRNG(nprocs()+1)
    path_part = @sprintf("Dis Tm%.2f dt%.3f del%d Obs%s Pr%.3f/", LP.Tmax, LP.tstep, LP.delimiter, get_string(LP.rules.str_vec), p)
    filenames = Vector{String}()
    A = HybridApprox(LP.A.L, copy(LP.A.q_spins), LP.A.name)
    l = LP.cb.affect!.save_func.len
    tspan = (0.0, LP.Tmax*LP.delimiter)
    t_means = collect(0.0:LP.tstep:(LP.Tmax*LP.delimiter))
    saved_values = SavedValues(t_means, [zeros(l) for i in eachindex(t_means)])
    for thread = (offset1+1):(offset1+nprocs())
        push!(filenames, get_string(nTrials, thread, alg))
    end
    log_o, save_os = set_Logging(LP, nTrials, logger, LP.A, path_part, filenames)
    cf_vals0 = LP.cf_vals
    @sync for i in workers()
        @spawnat i begin
                       global cf_vals = cf_vals0
                       global RNG = RNGarray[i]
                   end
    end
    iter = get_status(log_o)
    offset = get_offset(log_o)
    if iter > 0
        load_data!(cf_vals0, save_os[1])
        for i in workers()
            @spawnat i load_data!(cf_vals, save_os[i])
        end
    end
    t0 = time()-offset
    while iter<nTrials
        t1 = time()
        iter += 1
        q_spins = randsubseq(RNGarray[end], LP.A.q_spins, 1-p)
        A.q_spins = q_spins
        A.q_num = length(q_spins)
        ARHS = build_RHS_function(A)
        rules, OSET = set_CorrFuncs(A, links, axes)
        @sync begin
            for i in workers()
                @spawnat i dh_run(iter, A, ARHS, OSET, rules, t_means, saved_values, tspan, p, cf_vals, alg, save_os[i], RNG; kwargs...)
            end
            @spawnat 1 dh_run(iter, A, ARHS, OSET, rules, t_means, saved_values, tspan, p, cf_vals0, alg, save_os[1], RNGarray[1]; kwargs...)
        end
        t2 = time()
        output(log_o, @sprintf("%10d %20.5f %20.5f %10d", iter, t2-t1, (t2-t0)/iter, nTrials - iter))
    end
end

function dh_simulation(LP::AbstractLProblem{Approx}, p::Float64, nTrials::Int64, alg::OrdinaryDiffEqAlgorithm, logger::Logger{T1,T2} = Logger{:no, :cmd}(), interval::Int64 = 1, offset1::Int64 = 0; links = (:all,:allq), axes = [:x, :y, :z], kwargs...) where {Approx <: CompositeApproximation{<:HybridApprox, <:PureClassicalApprox}, T1,T2}
    A = LP.A
    Tmax = LP.Tmax
    delimiter = LP.delimiter
    tstep = LP.tstep
    A_temp = CompositeApproximation(HybridApprox(A.A1.L, copy(A.A1.q_spins), A.A1.name), A.A2, A.gamma)
    #rules, OSET = set_CorrFuncs(A, links, axes)
    l = LP.cb.affect!.save_func.len
    tspan = (0.0, Tmax*delimiter)
    t_means = collect(0.0:tstep:(Tmax*delimiter))
    t_cors = collect(0.0:tstep:Tmax)
    saved_values = SavedValues(t_means, [zeros(l) for i in eachindex(t_means)])
    ##############################
    assert(nTrials > 0)
    RNGarray = initRNG(nprocs()+1)
    path_part = @sprintf("Dis Tm%.2f dt%.3f del%d Obs%s Pr%.3f/", Tmax, tstep, delimiter, get_string(LP.rules.str_vec), p)
    filenames = Vector{String}()
    for thread = (offset1+1):(offset1+nprocs())
        push!(filenames, get_string(nTrials, thread, alg))
    end
    log_o, save_os = set_Logging(LP, nTrials, logger, A, path_part, filenames)
    cf_vals0 = LP.cf_vals
    @sync for i in workers()
        @spawnat i begin
                       global cf_vals = cf_vals0
                       global RNG = RNGarray[i]
                   end
    end
    iter = get_status(log_o)
    offset = get_offset(log_o)
    if iter > 0
        load_data!(cf_vals0, save_os[1])
        for i in workers()
            @spawnat i load_data!(cf_vals, save_os[i])
        end
    end
    t0 = time()-offset
    ##############################
    while iter<nTrials
        t1 = time()
        iter += 1
        q_spins = randsubseq(RNGarray[end], A.A1.q_spins, 1-p)
        A_temp.A1.q_spins = q_spins
        A_temp.A1.q_num = length(q_spins)
        print(length(q_spins), "\n")
        ARHS = build_RHS_function(A_temp)
        rules, OSET = set_CorrFuncs(A_temp.A1, links, axes)
        @sync begin
            for i in workers()
                @spawnat i dh_run(iter, A_temp, ARHS, OSET, rules, t_means, saved_values, tspan, p, cf_vals, alg, save_os[i], RNG; kwargs...)
            end
            @spawnat 1 dh_run(iter, A_temp, ARHS, OSET, rules, t_means, saved_values, tspan, p, cf_vals0, alg, save_os[1], RNGarray[1]; kwargs...)
        end
        t2 = time()
        output(log_o, @sprintf("%10d %20.5f %20.5f %10d", iter, t2-t1, (t2-t0)/iter, nTrials - iter))
    end
end

function dh_run(iter::Int64, A::abstractApproximation, ARHS::AbstractRHSFunction, OSET::ObservablesSet, rules::ConvRules, t_means::Vector{Float64}, saved_values::SavedValues, tspan::NTuple{2, Float64}, p::Float64, cf_vals::T, alg::OrdinaryDiffEqAlgorithm, save_o::Saver, RNG::AbstractRNG; abstol_i = 1e-14, reltol_i = 1e-6, kwargs...) where T <: Union{CFVals, CFVals0}
    u0 = randomDState(A, p, RNG)
    #print(u0, "\n")
    prob = ODEProblem(ARHS, u0, tspan)
    cb = SavingCallback(OSET, saved_values, u0; saveat = t_means)
    solve(prob, alg; save_everystep = false, callback = cb, abstol = abstol_i, reltol = reltol_i, kwargs...)
    calculateCorrFunc!(A, saved_values, rules, cf_vals, iter)
    save_data(save_o, cf_vals, iter)
    #print(cf_vals, "\n")
end

function randomDState(A::PureClassicalApprox, p::Float64, RNG::AbstractRNG = Base.GLOBAL_RNG)
    cl_state = randomState(A, RNG)
    discarded = randsubseq(RNG, collect(1:A.L.tot_spin_num), p)
    for i in discarded
        cl_state[i,:] .= 0.0
    end
    return cl_state
end

function randomDState(A::HybridApprox, p::Float64, RNG::AbstractRNG = Base.GLOBAL_RNG)
    Dh = 2^A.q_num
    cl_state = random_cl_state(A.cl_num, RNG)
    discarded = randsubseq(RNG, collect(1:A.cl_num), p)
    for i in discarded
        cl_state[i,:] .= 0.0
    end
    return ArrayPartition(random_q_state(Dh, RNG), cl_state)
end

function randomDState(A::CompositeApproximation{T1,T2}, p::Float64, RNG::AbstractRNG = Base.GLOBAL_RNG) where {T1<:HybridApprox, T2<:PureClassicalApprox}
    Dh = 2^A.A1.q_num
    cl_state = random_cl_state(A.A1.cl_num, RNG)
    discarded = randsubseq(RNG, collect(1:A.A1.cl_num), p)
    for i in discarded
        cl_state[i,:] .= 0.0
    end
    ArrayPartition(random_q_state(Dh, RNG), cl_state, randomState(A.A2, RNG))
end

function getrsubseq(A::AbstractArray{T}, p::Float64, RNG::AbstractRNG = Base.GLOBAL_RNG) where {T}
    B = Vector{T}()
    for i in eachindex(A)
        if rand(RNG) < p
            push!(B, A[i])
        end
    end
    return B
end

export dh_simulation


text:
using PyPlot
mutable struct Single_Stack
    L::Float64
    N::Int64
    t::Float64

    ρ::Array{Float64,1}
    w::Array{Float64,1}
    ρθ::Array{Float64,1}

    _grav::Float64

    Π::Function

end

# compute tendency
function spatial_residual(ss::Single_Stack, ρ, w, ρθ; t = 0.0)

    N, L = ss.N, ss.L
    Δz = L/N
    _grav = ss._grav
    Π = ss.Π

    ρ_res, w_res, ρθ_res = copy(ρ), copy(w), copy(ρθ)

    for i = 1:N
        # ρ₊ = (i == N ? ρ[i] : (ρ[i] + ρ[i+1])/2.0)
        # ρ₋ = (i == 1 ? ρ[i] : (ρ[i] + ρ[i-1])/2.0)

        ρ₊ = (i == N ? 0.0 : (ρ[i] + ρ[i+1])/2.0)
        ρ₋ = (i == 1 ? 0.0 : (ρ[i] + ρ[i-1])/2.0)

        ρ_res[i] = -(w[i+1]ρ₊ - w[i]ρ₋)/Δz

        # ρθ₊ = (i == N ? ρθ[i] : (ρθ[i] + ρθ[i+1])/2.0)
        # ρθ₋ = (i == 1 ? ρθ[i] : (ρθ[i] + ρθ[i-1])/2.0)

        ρθ₊ = (i == N ? 0.0 : (ρθ[i] + ρθ[i+1])/2.0)
        ρθ₋ = (i == 1 ? 0.0 : (ρθ[i] + ρθ[i-1])/2.0)

        ρθ_res[i] = -(w[i+1]ρθ₊ - w[i]ρθ₋)/Δz
    end

    w_res[1] = w_res[N+1] = 0
    for i = 2:N
        w_res[i] = -1/2*(ρθ[i-1]/ρ[i-1] + ρθ[i]/ρ[i]) *(Π(ρθ[i]) - Π(ρθ[i-1]))/Δz - _grav
    end

    return ρ_res,w_res, ρθ_res
end

# update with RK4
function temporal_update!(ss::Single_Stack, Δt::Float64)
    ρ, w, ρθ = ss.ρ, ss.w, ss.ρθ

    ρ_res1, w_res1, ρθ_res1 = spatial_residual(ss::Single_Stack, ρ, w, ρθ; t = ss.t)
    ρ_res2, w_res2, ρθ_res2 = spatial_residual(ss::Single_Stack, ρ + Δt/2*ρ_res1, w + Δt/2*w_res1, ρθ + Δt/2*ρθ_res1; t = ss.t+Δt/2)
    ρ_res3, w_res3, ρθ_res3 = spatial_residual(ss::Single_Stack, ρ + Δt/2*ρ_res2, w + Δt/2*w_res2, ρθ + Δt/2*ρθ_res2; t = ss.t+Δt/2)
    ρ_res4, w_res4, ρθ_res4 = spatial_residual(ss::Single_Stack, ρ + Δt*ρ_res3, w + Δt*w_res3, ρθ + Δt*ρθ_res3; t = ss.t+Δt)

    ρ  .+=  Δt/6*(ρ_res1  + 2ρ_res2  + 2ρ_res3  + ρ_res4)
    w  .+=  Δt/6*(w_res1  + 2w_res2  + 2w_res3  + w_res4)
    ρθ .+=  Δt/6*(ρθ_res1 + 2ρθ_res2 + 2ρθ_res3 + ρθ_res4)

    ss.t += Δt
    @info "ρ ∈ ",   minimum(ρ), maximum(ρ)
    @info "w ∈ ",   minimum(w), maximum(w)
    @info "ρθ ∈ ",  minimum(ρθ), maximum(ρθ)

end

function compute_wave_speed(ρ, w, ρθ, _R_m, _C_p, γ, _MSLP)
    p =  ( ρθ * _R_m / _MSLP^(_R_m/_C_p) ).^(γ)
    c = sqrt.(γ*p./ρ)

    return maximum(w) + maximum(c)
end

# Decaying temperature profile
function DecayingTemperatureProfile(z::Float64; T_virt_surf, T_min_ref, _R_d, _grav, _MSLP)
    # Scale height for surface temperature
    H_sfc = _R_d * T_virt_surf / _grav
    H_t = H_sfc

    z′ = z / H_t
    tanh_z′ = tanh(z′)

    ΔTv =  T_virt_surf -  T_min_ref
    Tv =  T_virt_surf - ΔTv * tanh_z′

    ΔTv′ = ΔTv / T_virt_surf
    p = -H_t * (z′ + ΔTv′ * (log(1 - ΔTv′ * tanh_z′) - log(1 + tanh_z′) + z′))
    p /= H_sfc * (1 - ΔTv′^2)
    p = _MSLP * exp(p)
    ρ = p/(_R_d*Tv)
    return (Tv, p, ρ)
end

# update discrete_hydrostatic_balance!
function discrete_hydrostatic_balance!(ρ, w, ρθ, Δz::Float64, _grav::Float64, Π::Function)
    for i = 1:length(ρ)-1
        ρ[i+1] = ρθ[i+1]/(-2Δz*_grav/(Π(ρθ[i+1]) - Π(ρθ[i])) - ρθ[i]/ρ[i])
    end
end

###################################### Start
ndays = 10

# whether use discrete hydrostatic balance correction
DHB = true

N = 30
L = 30e3
Δz = L/N

_MSLP = 1e5
_grav = 9.8
_R_m = 287.058
γ = 1.4
_C_p = _R_m * γ / (γ - 1)
_C_v = _R_m / (γ - 1)
_R_d = _R_m

Π = (ρθ) -> _C_p*(_R_d*ρθ/_MSLP)^(_R_m/_C_v)

###################################### Initialization
ρ = zeros(Float64, N)
w = zeros(Float64, N+1)
ρθ = zeros(Float64, N)

# set the bottom velocity
# w[1] = 1.0

for i = 1:N
    z = (i-0.5)*L/N
    Tvᵢ, pᵢ, ρᵢ = DecayingTemperatureProfile(z; T_virt_surf = 280.0, T_min_ref = 230.0, _R_d = _R_d, _grav = _grav, _MSLP = _MSLP)
    ρ[i]  = ρᵢ
    ρθ[i] = ρᵢ*Tvᵢ*(_MSLP/pᵢ)^(_R_m/_C_p)
end

if DHB
    discrete_hydrostatic_balance!(ρ, w, ρθ, Δz, _grav, Π)
end

ss = Single_Stack(L, N, 0.0, copy(ρ), copy(w), copy(ρθ), _grav, Π)
max_wave_speed = compute_wave_speed(ρ, w, ρθ, _R_m, _C_p, γ, _MSLP)

@info "maximum Δt (cfl = 1) = ",  Δz / max_wave_speed

###################################### Advance
Δt = 3.0
Nt = Int64(round(86400  * ndays / Δt))
t_end
for i = 1:Nt
    temporal_update!(ss, Δt)
end
T = Nt * Δt

zz = Array(LinRange(Δz/2, L - Δz/2, N))
zz_h = Array(LinRange(0, L, N+1))

fig, (ax1, ax2, ax3) = PyPlot.subplots(ncols = 3, nrows=1, sharex=false, sharey=false, figsize=(12,6))
ax1.plot(ρ, zz,  "-o", fillstyle = "none", label = "T=0")
ax1.plot(ss.ρ, zz, "-", fillstyle = "none", label = "T=$(ndays) days")
ax1.legend()
ax1.set_ylabel("ρ")

ax2.plot(w, zz_h, "-o", fillstyle = "none", label = "T=0")
ax2.plot(ss.w, zz_h,  "-", fillstyle = "none", label = "T=$(ndays) days")
ax2.legend()
ax2.set_ylabel("w")

ax3.plot(ρθ, zz,  "-o", fillstyle = "none", label = "T=0")
ax3.plot(ss.ρθ, zz, "-", fillstyle = "none", label = "T=$(ndays) days")
ax3.legend()
ax3.set_ylabel("ρθ")

fig.tight_layout()
fig.savefig("HB-Theta-$(ndays)-HB-$(DHB).png")

ext: jl


text:
using ConservationLawsParticles
using RecursiveArrayTools
using Test
using TestSetExtensions

@testset ExtendedTestSet "ConservationLawsParticles.jl" begin
    @testset "densities" begin
        include("test_densities.jl")
    end

    @testset "velocities" begin
        include("test_velocities.jl")
    end

    @testset "utils" begin
        @testset "empty_like" begin
            @test length(@test_deprecated empty_like([1,2,3])) == 3
            @test typeof(@test_deprecated empty_like([1,2,3])) == Vector{Int}
        end
        @testset "@time_independent" begin
            @time_independent V_time_indep(x, y) = x + y
            @test V_time_indep(2, 3) == 5
            @test V_time_indep(0, 2, 3) == 5
            @test V_time_indep(:time, 2, 3) == 5
            @test V_time_indep(nothing, 2, 3) == 5
        end
    end
end


text:
using TensND, LinearAlgebra, SymPy, Tensors, OMEinsum, Rotations
sympy.init_printing(use_unicode=true)

CS = coorsys_cylindrical() ; r, θ, z = getcoords(CS) ; 𝐞ʳ, 𝐞ᶿ, 𝐞ᶻ = unitvec(CS) ;
ℬ = get_normalized_basis(CS)
@set_coorsys CS

ξʳ, ξᶻ = SymFunction("ξʳ, ξᶻ", real = true)
𝛏 = ξʳ(r,z)*𝐞ʳ+ξᶻ(r,z)*𝐞ᶻ
𝛜 = SYMGRAD(𝛏)


text:
using ModelingToolkit, LinearAlgebra, Test

@parameters t k[1:13]
@variables A(t) B(t) C(t) D(t)
rxs = [Reaction(k[1], nothing, [A]),            # 0 -> A
       Reaction(k[2], [B], nothing),            # B -> 0
       Reaction(k[3],[A],[C]),                  # A -> C
       Reaction(k[4], [C], [A,B]),              # C -> A + B
       Reaction(k[5], [C], [A], [1], [2]),      # C -> A + A
       Reaction(k[6], [A,B], [C]),              # A + B -> C
       Reaction(k[7], [B], [A], [2], [1]),      # 2B -> A
       Reaction(k[8], [A,B], [A,C]),            # A + B -> A + C
       Reaction(k[9], [A,B], [C,D]),            # A + B -> C + D
       Reaction(k[10], [A], [C,D], [2], [1,1]), # 2A -> C + D
       Reaction(k[11], [A], [A,B], [2], [1,1]), # 2A -> A + B
       Reaction(k[12], [A,B,C], [C,D], [1,3,4], [2, 3]), # A+3B+4C -> 2C + 3D
       Reaction(k[13]*A/(2+A), [A], nothing; only_use_rate=true)  # A -> 0 with custom rate
       ]
rs = ReactionSystem(rxs,t,[A,B,C,D],k)
odesys = convert(ODESystem,rs)
sdesys = convert(SDESystem,rs)

# hard coded ODE rhs
function oderhs(u,k,t)
       A = u[1]; B = u[2]; C = u[3]; D = u[4];
       du = zeros(eltype(u),4)
       du[1] = k[1] - k[3]*A + k[4]*C + 2*k[5]*C - k[6]*A*B + k[7]*B^2/2 - k[9]*A*B - k[10]*A^2 - k[11]*A^2/2 - k[12]*A*B^3*C^4/144 - k[13]*A/(2+A)
       du[2] = -k[2]*B + k[4]*C - k[6]*A*B - k[7]*B^2 - k[8]*A*B - k[9]*A*B + k[11]*A^2/2 - 3*k[12]*A*B^3*C^4/144
       du[3] = k[3]*A - k[4]*C - k[5]*C  + k[6]*A*B + k[8]*A*B + k[9]*A*B + k[10]*A^2/2 - 2*k[12]*A*B^3*C^4/144
       du[4] = k[9]*A*B + k[10]*A^2/2 + 3*k[12]*A*B^3*C^4/144
       du
end

# sde noise coefs
function sdenoise(u,k,t)
       A = u[1]; B = u[2]; C = u[3]; D = u[4];
       G = zeros(eltype(u),length(k),length(u))
       z = zero(eltype(u))

       G = [sqrt(k[1]) z  z z;
            z -sqrt(k[2]*B) z z;
            -sqrt(k[3]*A) z sqrt(k[3]*A) z;
            sqrt(k[4]*C) sqrt(k[4]*C) -sqrt(k[4]*C) z;
            2*sqrt(k[5]*C) z -sqrt(k[5]*C) z;
            -sqrt(k[6]*A*B) -sqrt(k[6]*A*B) sqrt(k[6]*A*B) z;
            sqrt(k[7]*B^2/2) -2*sqrt(k[7]*B^2/2) z z;
            z -sqrt(k[8]*A*B) sqrt(k[8]*A*B) z;
            -sqrt(k[9]*A*B) -sqrt(k[9]*A*B) sqrt(k[9]*A*B) sqrt(k[9]*A*B);
            -2*sqrt(k[10]*A^2/2) z sqrt(k[10]*A^2/2) sqrt(k[10]*A^2/2);
            -sqrt(k[11]*A^2/2) sqrt(k[11]*A^2/2) z z;
            -sqrt(k[12]*A*B^3*C^4/144) -3*sqrt(k[12]*A*B^3*C^4/144) -2*sqrt(k[12]*A*B^3*C^4/144) 3*sqrt(k[12]*A*B^3*C^4/144);
            -sqrt(k[13]*A/(2+A)) z z z]'

       return G
end

# test by evaluating drift and diffusion terms
p  = rand(length(k))
u  = rand(length(k))
t  = 0.
du = oderhs(u,p,t)
G  = sdenoise(u,p,t)
sdesys = convert(SDESystem,rs)
sf = SDEFunction{false}(sdesys, states(rs), parameters(rs))
du2 = sf.f(u,p,t)
@test norm(du-du2) < 100*eps()
G2 = sf.g(u,p,t)
@test norm(G-G2) < 100*eps()


text:
# Character Set Encoding support
#
# Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones
# Licensed under MIT License, see LICENSE.md

@api public CSE, "@cse"
@api public! basecse
@api develop cse_types

struct CSE{CS, ENC}  end

"""List of installed character set encodings"""
const cse_types = []

CSE(cs, e)  = CSE{CharSet(cs), Encoding(e)}()

macro cse(cs, e)
    :(CSE{CharSet{$(quotesym(cs)), $(quotesym(e))}()})
end

const _CSE{U} = Union{CharSet{U}, Encoding{U}} where {U}

print(io::IO, ::S) where {S<:_CSE{U}} where {U} = print(io, U)

show(io::IO, ::Type{CSE{CS,E}}) where {S,T,CS<:CharSet{S},E<:Encoding{T}} =
    print(io, "CSE{", string(S), ", ", string(T), "}")
print(io::IO, ::T) where {S,U,CS<:CharSet{S},E<:Encoding{U},T<:CSE{CS,E}} =
    (show(io, T); print(io, "()"))

# Definition of built-in CSEs (Character Set Encodings)

codeunit(::Type{<:CSE}) = UInt8
basecse(::Type{C}) where {C<:CSE} = C

for lst in cse_info
    nam, cu = lst
    cse = symstr(nam, "CSE")
    if length(lst) > 2
        csnam = symstr(lst[3], "CharSet")
        enc = cu === UInt8 ? :UTF8Encoding : :NativeUTF16
    else
        csnam = symstr(nam, "CharSet")
        enc = cu === UInt8 ? :Native1Byte : cu === UInt16 ? :Native2Byte : :Native4Byte
    end
    @eval const $(symstr(nam, "CSE")) = CSE{$csnam, $enc}
    @eval show(io::IO, ::Type{$cse}) = print(io, $(String(cse)))
    cu === UInt8 || (@eval codeunit(::Type{$cse}) = $cu)
    @eval push!(cse_types, $cse)
    if String(nam)[1] == '_'
        @eval basecse(::Type{$cse}) = $(symstr(String(nam)[2:end], "CSE"))
        @eval @api develop $cse
    else
        @eval @api public $cse
    end
end

# Various useful groups of character set types

# These should be done via traits
const Binary_CSEs   = Union{Text1CSE, BinaryCSE}
const Latin_CSEs    = Union{LatinCSE, _LatinCSE}
const UTF8_CSEs     = Union{UTF8CSE,  RawUTF8CSE}
const UCS2_CSEs     = Union{UCS2CSE,  _UCS2CSE}
const UTF32_CSEs    = Union{UTF32CSE, _UTF32CSE}
const SubSet_CSEs   = Union{_LatinCSE, _UCS2CSE, _UTF32CSE}

const Byte_CSEs     = Union{ASCIICSE, Binary_CSEs, Latin_CSEs, UTF8_CSEs} # 8-bit code units
const Word_CSEs     = Union{Text2CSE, UCS2CSE, _UCS2CSE, UTF16CSE} # 16-bit code units
const Quad_CSEs     = Union{Text4CSE, UTF32CSE, _UTF32CSE}         # 32-bit code units

@api develop Binary_CSEs, Latin_CSEs, UTF8_CSEs, UCS2_CSEs, UTF32_CSEs, SubSet_CSEs,
             Byte_CSEs, Word_CSEs, Quad_CSEs

cse(::Type{<:AbstractString}) = RawUTF8CSE     # allows invalid sequences
cse(::Type{<:SubString{T}}) where {T} = basecse(T)
cse(::T) where {T<:AbstractString} = cse(T)

basecse(::Type{T}) where {T<:AbstractString} = basecse(cse(T))
basecse(::T) where {T<:AbstractString} = basecse(cse(T))

# Get charset based on CSE
charset(::Type{<:CSE{CS}}) where {CS<:CharSet} = CS
charset(::Type{T}) where {T<:AbstractString} = charset(cse(T))
charset(::T) where {T<:AbstractString} = charset(cse(T))

# Get encoding based on CSE
encoding(::Type{<:CSE{CS,E}}) where {CS<:CharSet,E<:Encoding} = E

# Promotion rules for character set encodings

promote_rule(::Type{Text2CSE}, ::Type{Text1CSE}) = Text2CSE
promote_rule(::Type{Text4CSE}, ::Type{Text1CSE}) = Text4CSE
promote_rule(::Type{Text4CSE}, ::Type{Text2CSE}) = Text4CSE

promote_rule(::Type{T}, ::Type{ASCIICSE}) where {T<:CSE} = T
promote_rule(::Type{T}, ::Type{<:Latin_CSEs}
             ) where {T<:Union{UTF8CSE,UTF16CSE,UCS2_CSEs,UTF32_CSEs}} = T

promote_rule(::Type{T}, ::Type{_LatinCSE}) where {T<:Union{ASCIICSE,LatinCSE}} = LatinCSE

promote_rule(::Type{T}, ::Type{_UCS2CSE}) where {T<:Union{ASCIICSE,Latin_CSEs,UCS2CSE}} = UCS2CSE
promote_rule(::Type{T}, ::Type{_UTF32CSE}) where {T<:CSE} = UTF32CSE
promote_rule(::Type{T}, ::Type{UTF32CSE}) where {T<:UCS2_CSEs} = UTF32CSE


text:
function part1(input, number_of_days = "0")
    instructions = readlines(input)
    black_tiles = get_initial_black_tiles(instructions)
    return length(black_tiles)
end

function part2(input, number_of_days = "100")
    instructions = readlines(input)
    black_tiles = get_initial_black_tiles(instructions)
    for _ = 1:parse(Int, number_of_days)
        black_tiles = flip_tiles(black_tiles)
    end
    return length(black_tiles)
end

function get_initial_black_tiles(instructions)
    black_tiles = Set{Tuple{Int, Int}}()
    for instruction in instructions
        location = parse_instruction(instruction)
        if location in black_tiles
            delete!(black_tiles, location)
        else
            push!(black_tiles, location)
        end
    end
    return black_tiles
end

function parse_instruction(instruction)
    e, ne = 0, 0
    i = 1
    while i <= length(instruction)
        if instruction[i] == 'w'
            e -= 1
            i += 1
        elseif instruction[i] == 'e'
            e += 1
            i += 1
        elseif instruction[i] == 'n'
            if instruction[i + 1] == 'w'
                e -= 1
                ne += 1
                i += 2
            elseif instruction[i + 1] == 'e'
                ne += 1
                i += 2
            end
        elseif instruction[i] == 's'
            if instruction[i + 1] == 'w'
                ne -= 1
                i += 2
            elseif instruction[i + 1] == 'e'
                e += 1
                ne -= 1
                i += 2
            end
        end
    end
    return e, ne
end

function flip_tiles(tiles)
    neighbor_counts = Dict{Tuple{Int, Int}, Int}()
    for tile in tiles
        neighbor_counts[tile] = get(neighbor_counts, tile, 0) + 1
        for Δ in ((1, 0), (-1, 0), (-1, 1), (0, 1), (0, -1), (1, -1))
            neighbor_counts[tile .+ Δ] = get(neighbor_counts, tile .+ Δ, 0) + 2
        end
    end
    new_tiles = Set{Tuple{Int, Int}}()
    for (tile, neighbor_count) in neighbor_counts
        if 3 <= neighbor_count <= 5
            push!(new_tiles, tile)
        end
    end
    return new_tiles
end


text:
using Base:func_for_method_checked
import MPI

include("../types/types.jl")
include("../setting/constants.jl")
include("../utils/readfiles.jl")
include("../utils/kdtree.jl")
include("../utils/parse_commandline_ppm.jl")

function array_split_mpi(input_array, num_split, rank)
    div_result = div(length(input_array), num_split)
    mod_result = mod(length(input_array), num_split)
    startindex = nothing
    endindex = nothing
    if rank <= mod_result
        startindex = (div_result + 1) * (rank - 1) + 1
        endindex = startindex + div_result
    else
        startindex = (div_result + 1) * mod_result + div_result * (rank - mod_result - 1) + 1
        endindex = startindex + div_result - 1
    end
    return input_array[startindex:endindex]
end

function generate_profile_points(latnpts, lonnpts, vnpts, lon1, lat1, dep1, lon2, lat2, dep2, comm, latproc, lonproc)
    # MPI parameters
    rank = MPI.Comm_rank(comm)
    size_mpi = MPI.Comm_size(comm)

    # get ranges for the three directions
    rank_lat = rank % latproc + 1
    rank_lon = div(rank, latproc) + 1
    coor_lat = array_split_mpi(1:latnpts, latproc, rank_lat)
    coor_lon = array_split_mpi(1:lonnpts, lonproc, rank_lon)

    # init rθϕ_new
    ngll_new_this_rank = length(coor_lat) * length(coor_lon) * vnpts

    deplatlon_new_this_rank = zeros(Float64, 3, ngll_new_this_rank)
    # fill in deplatlon_new of this rank
    latnpts_this_rank = length(coor_lat)
    lonnpts_this_rank = length(coor_lon)
    for (vindex, dep) in enumerate(range(dep1, stop=dep2, length=vnpts))
        for (latindex, lat) in enumerate(range(lat1, stop=lat2, length=latnpts)[coor_lat])
            for (lonindex, lon) in enumerate(range(lon1, stop=lon2, length=lonnpts)[coor_lon])
                id = (vindex - 1) * latnpts_this_rank * lonnpts_this_rank + (latindex - 1) * lonnpts_this_rank + lonindex
                deplatlon_new_this_rank[:,id] = [dep,lat,lon]
            end
        end
    end

    xyz_new = zeros(Float64, 3, ngll_new_this_rank)
    # convert dep,lat,lon to x,y,z
    for id in 1:ngll_new_this_rank
        (dep, lat, lon) = deplatlon_new_this_rank[:,id]
        R_EARTH_KM = 6371.e0

        # normalized r
        r = (R_EARTH_KM - dep) / R_EARTH_KM
        θ = 90 - lat
        ϕ = lon

        # get x,y,z
        z = r * cosd(θ)
        h = r * sind(θ)
        x = h * cosd(ϕ)
        y = h * sind(ϕ)

        xyz_new[:,id] = [x,y,z]
    end

    return xyz_new
end

function run_interp(command_args::Dict{String,Any}, comm::MPI.Comm)
    # MPI
    rank = MPI.Comm_rank(comm)
    isroot = (rank == 0)

    # * init some variables
    nproc_old = command_args["nproc_old"]
    old_mesh_dir = command_args["old_mesh_dir"]
    old_model_dir = command_args["old_model_dir"]
    model_tags = command_args["model_tags"]
    output_file = command_args["output_file"]
    region = command_args["region"]
    npts = command_args["npts"]
    nproc = command_args["nproc"]
    (lon1, lat1, lon2, lat2, dep1, dep2) = parse.(Float64, split(region, "/"))
    (latnpts, lonnpts, vnpts) = parse.(Int64, split(npts, "/"))
    (latproc, lonproc) = parse.(Int64, split(nproc, "/"))

    # * init some variables
    mesh_old = sem_mesh_data()

    # * parse model tags
    model_names = String.(split(model_tags, ","))
    nmodel = length(model_names)

    if isroot
        @info "# nmodel=$(nmodel)"
        @info "# model names=$(model_names)"
    end

    # * interpolate old mesh/model onto new mesh

    # typical element size at surface for old mesh
    # read from constants.jl for old mesh
    typical_size = deg2rad(max(ANGULAR_WIDTH_XI_IN_DEGREES_VAL / NEX_XI_VAL, ANGULAR_WIDTH_ETA_IN_DEGREES_VAL / NEX_ETA_VAL)) * R_UNIT_SPHERE
    max_search_dist = 10.0 * typical_size
    max_misloc = typical_size / 4.0

    # get xyz_new
    xyz_new = generate_profile_points(latnpts, lonnpts, vnpts, lon1, lat1, dep1, lon2, lat2, dep2, comm, latproc, lonproc)

    ngll_new_this_rank = size(xyz_new)[2]
    idoubling_new = zeros(Int64, ngll_new_this_rank)
    idoubling_new .= IFLAG_DUMMY

    # * initialize variables for interpolation
    location_1slice = Array{sem_mesh_location}(undef, ngll_new_this_rank)
    stat_final = zeros(Int64, ngll_new_this_rank)
    misloc_final = zeros(Float64, ngll_new_this_rank)

    stat_final .= -1
    misloc_final .= HUGEVAL

    model_interp_this_rank = zeros(nmodel, ngll_new_this_rank)
    # make model_interp_this_rank as 99999 to remove the poitns that hasn't been interpolated
    model_interp_this_rank .= 9999999.0

    # * loop each slices of the old mesh
    flag = true

    # * loop for all points in xyz_new
    for iproc_old = 0:nproc_old - 1
        @info "[$(rank)]# iproc=$(iproc_old)"
        # read old mesh slice
        mesh_old = sem_mesh_read(old_mesh_dir, iproc_old)
        nspec_old = mesh_old.nspec
        # test if the new and old mesh slices are separated apart
        min_dist = HUGEVAL
        for ispec = 1:nspec_old
            iglob = mesh_old.ibool[MIDX,MIDY,MIDZ,ispec]
            old_x = mesh_old.xyz_glob[1,iglob]
            old_y = mesh_old.xyz_glob[2,iglob]
            old_z = mesh_old.xyz_glob[3,iglob]
            dist_this_spec = sqrt(minimum(@. (xyz_new[1,:] - old_x)^2 + (xyz_new[2,:] - old_y)^2 + (xyz_new[3,:] - old_z)^2))
            min_dist = min(min_dist, dist_this_spec)
        end
        if min_dist > max_search_dist
            # @info "[$rank]#$(iproc_old) slices too far away, skip"
            continue
        end

        # read old model
        model_gll_old = zeros(Float64, nmodel, NGLLX, NGLLY, NGLLZ, nspec_old)
        sem_io_read_gll_file_n!(old_model_dir, iproc_old, model_names, nmodel, model_gll_old)

        # locate points in this mesh slice
        nnearest = 10
        location_1slice = sem_mesh_locate_kdtree2!(mesh_old, ngll_new_this_rank, xyz_new, idoubling_new, nnearest, max_search_dist, max_misloc, iproc_old)

        for igll = 1:ngll_new_this_rank
            if (stat_final[igll] == 1 && location_1slice[igll].stat == 1)
                @info "[$rank]# multi-located, $(xyz_new[:,igll])"
                continue
            end
            # for point located inside one element in the first time or closer to one element than located before
            if location_1slice[igll].stat == 1 || (location_1slice[igll].stat == 0 && location_1slice[igll].misloc < misloc_final[igll])
                for imodel = 1:nmodel
                    model_interp_this_rank[imodel,igll] = sum(location_1slice[igll].lagrange .* model_gll_old[imodel,:,:,:,location_1slice[igll].eid])
                end
                stat_final[igll] = location_1slice[igll].stat
                misloc_final[igll] = location_1slice[igll].misloc
            end
        end
    end

    # gather all ngll_new_this_rank to ngll_new

    # * write out gll files for this new mesh slice
    # get latnpts_this_rank and lonnpts_this_rank
    rank = MPI.Comm_rank(comm)
    size_mpi = MPI.Comm_size(comm)

    rank_lat = rank % latproc + 1
    rank_lon = div(rank, latproc) + 1
    coor_lat = array_split_mpi(1:latnpts, latproc, rank_lat)
    coor_lon = array_split_mpi(1:lonnpts, lonproc, rank_lon)

    latnpts_this_rank = length(coor_lat)
    lonnpts_this_rank = length(coor_lon)

    # output files
    if isroot
        run(`rm -rf $output_file`)
        run(`mkdir -p $output_file`)
    end
    MPI.Barrier(comm)

    open(output_file * "/$rank", "w") do io
        for (vindex, dep) in enumerate(range(dep1, stop=dep2, length=vnpts))
            for (latindex, lat) in enumerate(range(lat1, stop=lat2, length=latnpts)[coor_lat])
                for (lonindex, lon) in enumerate(range(lon1, stop=lon2, length=lonnpts)[coor_lon])
                    id = (vindex - 1) * latnpts_this_rank * lonnpts_this_rank + (latindex - 1) * lonnpts_this_rank + lonindex
                    thesize = size(model_interp_this_rank)[1]
                    if thesize == 1
                        write(io, "$lon $lat $dep $(model_interp_this_rank[1,id]) \n")
                    else
                        out_string = "$lon $lat $dep "
                        for item in model_interp_this_rank[:,id]
                            out_string *= "$item "
                        end
                        out_string *= "\n"
                        write(io, out_string)
                    end
                end
            end
        end
    end

end

function main()
    MPI.Init()
    comm = MPI.COMM_WORLD

    command_args = parse_commandline()
    run_interp(command_args, comm)

    MPI.Finalize()
end

main()

max_stars_repo_name: ziyixi/seisflow

max_issues_repo_name: ziyixi/seisflow

max_forks_repo_name: ziyixi/seisflow


text:
import CuArrays: cu
import YaoArrayRegister: _measure, measure, measure!, measure_collapseto!, measure_remove!
import YaoBase: batch_normalize!
import Yao: expect

export cpu, cu, GPUReg

cu(reg::ArrayReg{B}) where B = ArrayReg{B}(CuArray(reg.state))
cpu(reg::ArrayReg{B}) where B = ArrayReg{B}(collect(reg.state))
const GPUReg{B, T, MT} = ArrayReg{B, T, MT} where MT<:GPUArray

function batch_normalize!(s::CuSubArr, p::Real=2)
    p!=2 && throw(ArgumentError("p must be 2!"))
    s./=norm2(s, dims=1)
    s
end

@inline function tri2ij(l::Int)
    i = ceil(Int, sqrt(2*l+0.25)-0.5)
    j = l-i*(i-1)÷2
    i+1,j
end

############### MEASURE ##################
measure(rng::AbstractRNG, ::ComputationalBasis, reg::GPUReg{1}, ::AllLocs; nshots::Int=1) = _measure(rng::AbstractRNG, reg |> probs |> Vector, nshots)
# TODO: optimize the batch dimension using parallel sampling
function measure(rng::AbstractRNG, ::ComputationalBasis, reg::GPUReg{B}, ::AllLocs; nshots::Int=1) where B
    regm = reg |> rank3
    pl = dropdims(mapreduce(abs2, +, regm, dims=2), dims=2)
    _measure(rng::AbstractRNG, pl |> Matrix, nshots)
end

function measure_remove!(rng::AbstractRNG, ::ComputationalBasis, reg::GPUReg{B}, ::AllLocs) where B
    regm = reg |> rank3
    nregm = similar(regm, 1<<nremain(reg), B)
    pl = dropdims(mapreduce(abs2, +, regm, dims=2), dims=2)
    pl_cpu = pl |> Matrix
    res_cpu = map(ib->_measure(rng::AbstractRNG, view(pl_cpu, :, ib), 1)[], 1:B)
    res = CuArray(res_cpu)
    @inline function kernel(nregm, regm, res, pl)
        state = (blockIdx().x-1) * blockDim().x + threadIdx().x
        if state <= length(nregm)
            i,j = GPUArrays.gpu_ind2sub(nregm, state)
            r = Int(res[j])+1
            @inbounds nregm[i,j] = regm[r,i,j]/CUDAnative.sqrt(pl[r, j])
        end
        return
    end
    X, Y = cudiv(length(nregm))
    @cuda threads=X blocks=Y kernel(nregm, regm, res, pl)
    reg.state = reshape(nregm,1,:)
    res
end

function measure!(rng::AbstractRNG, ::ComputationalBasis, reg::GPUReg{B, T}, ::AllLocs) where {B, T}
    regm = reg |> rank3
    pl = dropdims(mapreduce(abs2, +, regm, dims=2), dims=2)
    pl_cpu = pl |> Matrix
    res_cpu = map(ib->_measure(rng::AbstractRNG, view(pl_cpu, :, ib), 1)[], 1:B)
    res = CuArray(res_cpu)

    @inline function kernel(regm, res, pl)
        state = (blockIdx().x-1) * blockDim().x + threadIdx().x
        if state <= length(regm)
            k,i,j = GPUArrays.gpu_ind2sub(regm, state)
            @inbounds rind = Int(res[j]) + 1
            @inbounds regm[k,i,j] = k==rind ? regm[k,i,j]/CUDAnative.sqrt(pl[k, j]) : T(0)
        end
        return
    end

    X, Y = cudiv(length(regm))
    @cuda threads=X blocks=Y kernel(regm, res, pl)
    res
end

function measure_collapseto!(rng::AbstractRNG, ::ComputationalBasis, reg::GPUReg{B, T}, ::AllLocs; config=0) where {B, T}
    regm = reg |> rank3
    pl = dropdims(mapreduce(abs2, +, regm, dims=2), dims=2)
    pl_cpu = pl |> Matrix
    res_cpu = map(ib->_measure(rng::AbstractRNG, view(pl_cpu, :, ib), 1)[], 1:B)
    res = CuArray(res_cpu)

    @inline function kernel(regm, res, pl, val)
        state = (blockIdx().x-1) * blockDim().x + threadIdx().x
        if state <= length(regm)
            k,i,j = GPUArrays.gpu_ind2sub(regm, state)
            @inbounds rind = Int(res[j]) + 1
            @inbounds k==val+1 && (regm[k,i,j] = regm[rind,i,j]/CUDAnative.sqrt(pl[rind, j]))
            CuArrays.sync_threads()
            @inbounds k!=val+1 && (regm[k,i,j] = 0)
        end
        return
    end

    X, Y = cudiv(length(regm))
    @cuda threads=X blocks=Y kernel(regm, res, pl, config)
    res
end

import YaoArrayRegister: insert_qubits!, join
function batched_kron(A::Union{CuArray{T1, 3}, Adjoint{<:Any, <:CuArray{T1, 3}}}, B::Union{CuArray{T2, 3}, Adjoint{<:Any, <:CuArray{T2, 3}}}) where {T1 ,T2}
    res = CuArrays.zeros(promote_type(T1,T2), size(A,1)*size(B, 1), size(A,2)*size(B,2), size(A, 3))
    @inline function kernel(res, A, B)
        state = (blockIdx().x-1) * blockDim().x + threadIdx().x
        i,j,b = GPUArrays.gpu_ind2sub(res, state)
        i_A = (i-1) ÷ size(B,1) + 1
        j_A = (j-1) ÷ size(B,2) + 1
        i_B = (i-1) % size(B,1) + 1
        j_B = (j-1) % size(B,2) + 1
        state <= length(res) && (@inbounds res[state] = A[i_A, j_A, b]*B[i_B, j_B, b])
        return
    end

    X, Y = cudiv(length(res))
    @cuda threads=X blocks=Y kernel(res, A, B)
    res
end

function join(reg1::GPUReg{B}, reg2::GPUReg{B}) where {B}
    s1 = reg1 |> rank3
    s2 = reg2 |> rank3
    state = batched_kron(s1, s2)
    ArrayReg{B}(reshape(state, size(state, 1), :))
end

export insert_qubits!
function insert_qubits!(reg::GPUReg{B}, loc::Int; nqubits::Int=1) where B
    na = nactive(reg)
    focus!(reg, 1:loc-1)
    reg2 = join(zero_state(nqubits; nbatch=B) |> cu, reg) |> relax! |> focus!((1:na+nqubits)...)
    reg.state = reg2.state
    reg
end

for FUNC in [:measure!, :measure_collapseto!, :measure_remove!]
    @eval function $FUNC(rng::AbstractRNG, op::AbstractBlock, reg::GPUReg, al::AllLocs; kwargs...) where B
        E, V = eigen!(mat(op) |> Matrix)
        ei = Eigen(E|>cu, V|>cu)
        $FUNC(rng::AbstractRNG, ei, reg, al; kwargs...)
    end
end


text:
using LinearAlgebra
function softmax(a)
    if ndims(a) == 1
        return exp.(a) ./ sum(exp.(a))
    else
        return exp.(a) ./ sum(exp.(a), dims = 2)
    end
end
function fit(x, t; alpha = 0.01, tau_max = 1000)
    function CEE(w, x, t) #this is Cross entropy Error
        p = softmax(x * w)
        grad = -(x' * (t - p))
        return grad / length(x[:, 1])
    end
    if size(x)[1] != size(t)[1]#Processing when the matrix is organized by dependent variable
        x = x'
    end
    x = hcat(ones(size(x)[1], 1), x)
    w = ones(size(x)[2], size(t)[2])
    for tau in 1 : tau_max
        grad = CEE(w, x, t)
        w -= alpha * grad
    end
    return w
end

function forecast(x, w)
    x = hcat(ones(size(x)[1], 1), x)
    return softmax(x * w)
end

#Pass so that the row is each data sample and the column is each feature.
function predict(x, w)
    x = hcat(ones(size(x)[1], 1), x)
    s = softmax(x * w)
    p = [findfirst(s[i, :] .== maximum(s[i, :])) for i in 1:size(s)[1]]
end


text:
for N in [Float64, Rational{Int}, Float32]
    # singleton
    S = Singleton(N[1, -2, 3, -4])
    H = Hyperrectangle(zeros(N, 4), N[1, 2, 3, 4])
    @test symmetric_interval_hull(S) == H
    @test box_approximation_symmetric(S) == H

    # linear map of singleton
    M = Diagonal(ones(N, 4))
    @test symmetric_interval_hull(M * S) == H

    # Minkowski sum of singletons
    H2 = Hyperrectangle(zeros(N, 4), 2 * N[1, 2, 3, 4])
    @test symmetric_interval_hull(S + S) == H2

    # linear map of hyperrectangle
    M = Diagonal(N[1, -1, -2, 2])
    H1 = Hyperrectangle(5 * ones(N, 4), N[1, 2, 3, 4])
    H2 = Hyperrectangle(zeros(N, 4), N[6, 7, 16, 18])
    @test symmetric_interval_hull(M * H1) == H2
end

# tests that only work with Float64 and Float32
for N in [Float64, Float32]
    # exponential map of singleton
    M = sparse(diagm(N[1, 1, 1]))
    E = SparseMatrixExp(M) * Singleton(N[1, 2, 3])
    H = Hyperrectangle(zeros(N, 3), N[ℯ, 2ℯ, 3ℯ])
    @test symmetric_interval_hull(E) ≈ H

    # exponential map of hyperrectangle
    M = sparse(diagm(N[1, 1, 1]))
    E = SparseMatrixExp(M) * Hyperrectangle(N[1, 2, 3], N[1, 2, 3])
    H = symmetric_interval_hull(E)
    @test H ≈ Hyperrectangle(zeros(N, 3), N[2ℯ, 4ℯ, 6ℯ])
    # matrix with negative entries
    M = sparse(diagm(N[1, -1]))
    E = SparseMatrixExp(M) * Hyperrectangle(N[1, 1], N[1, 1])
    H = symmetric_interval_hull(E)
    @test H ≈ Hyperrectangle(zeros(N, 2), N[2ℯ, 2ℯ^(-1)])
end


text:
using Documenter, InfiniteOpt, Distributions, Literate, Random

if !@isdefined(EXAMPLE_DIR)
    const EXAMPLE_DIR = joinpath(@__DIR__, "src", "examples")
end

# Sort and collect all the files of a certain extension in a folder
sorted_files(dir, ext) = sort!(filter!(f -> endswith(f, ext), readdir(dir)))

# Remove any files generated by Literate that are still in a subdirectory
function remove_literate_files(dir)
    for file in append!(sorted_files(dir, ".md"), sorted_files(dir, ".png"))
        rm(joinpath(dir, file))
    end
    return
end

# Post process corrections
function format_fix(content)
    return replace(content, "nothing #hide" => "")
end

# Generate all the markdown files via Literate for a particular doc section
# Expect the directory to only have sub-directories that contain files
function generate_pages(dir)
    for subdir in sort!(readdir(dir))
        current_dir = joinpath(dir, subdir)
        remove_literate_files(current_dir)
        for file in sorted_files(current_dir, ".jl")
            jl_filename = joinpath(current_dir, file)
            Random.seed!(0)
            Literate.markdown(jl_filename, current_dir; documenter = true,
                              execute = true, postprocess = format_fix)
        end
    end
    return
end

# Generate the list of subsection pages with relative paths
function subsection_entry(dir, subdir)
    files = sorted_files(joinpath(dir, subdir), ".md")
    return map(file -> joinpath(basename(dir), subdir, file), files)
end

# Generate the list of subsections and their pages with relative paths
function section_entry(dir)
    return map(subdir -> subdir => subsection_entry(dir, subdir), readdir(dir))
end

# Run and convert the examples
generate_pages(EXAMPLE_DIR)

# Make the documentation via Documenter.jl
makedocs(;
    pages = [
        "Home" => "index.md",
        "Installation" => "install.md",
        "Tutorials" => ["Quick Start" => "tutorials/quick_start.md"], # TODO add more tutorials
        "Examples" =>  section_entry(EXAMPLE_DIR),
        # "Background" => [], # TODO add the background
        "User Guide" => [
            "Infinite Models" => "guide/model.md",
            "Infinite Domains" => "guide/domains.md",
            "Infinite Parameters" => "guide/parameter.md",
            "Finite Parameters" => "guide/finite_parameter.md",
            "Variables" => "guide/variable.md",
            "Derivatives" => "guide/derivative.md",
            "Expressions" => "guide/expression.md",
            "Measures" => "guide/measure.md",
            "Objectives" => "guide/objective.md",
            "Constraints" => "guide/constraint.md",
            "Model Transcription" => "guide/transcribe.md",
            "Optimization" => "guide/optimize.md",
            "Results" => "guide/result.md"
            ],
        "API Manual" => [
            "Infinite Models" => "manual/model.md",
            "Infinite Domains" => "manual/domains.md",
            "Infinite Parameters" => "manual/parameter.md",
            "Finite Parameters" => "manual/finite_parameter.md",
            "Variables" => "manual/variable.md",
            "Derivatives" => "manual/derivative.md",
            "Expressions" => "manual/expression.md",
            "Measures" => "manual/measure.md",
            "Objectives" => "manual/objective.md",
            "Constraints" => "manual/constraint.md",
            "Model Transcription" => "manual/transcribe.md",
            "Optimization" => "manual/optimize.md",
            "Results" => "manual/result.md"
            ],
        "Development" => [
            "Extensions" => "develop/extensions.md",
            "Getting Started" => "develop/start_guide.md",
            "Style Guide" => "develop/style.md"
            ],
        ],
    repo = "https://github.com/pulsipher/InfiniteOpt.jl/blob/{commit}{path}#L{line}",
    sitename = "InfiniteOpt.jl",
    authors = "Joshua Pulsipher and Weiqi Zhang",
    doctest = true,
    checkdocs = :exports,
    linkcheck = true,
    linkcheck_ignore = [r"https://www.youtube.com/.*"],
    strict = true,
    format = Documenter.HTML(
        # See https://github.com/JuliaDocs/Documenter.jl/issues/868
        prettyurls = get(ENV, "CI", nothing) == "true",
        analytics = "UA-178297470-1",
        collapselevel = 1,
        assets = ["assets/extra_styles.css"],

    )
)

# Push the built HTML files to the correct GitHub branch if building via CI
deploydocs(;
    repo = "github.com/pulsipher/InfiniteOpt.jl",
    push_preview = true
)


text:
using EmpiricalBayes
using Test
using Distributions
using QuadGK

# Move these to other parts eventually

# test if convolution matrix works for small n..
A = convolution_matrix(DiscretizedNormalConvolutionProblem, range(-1, stop=1, length=4), range(-7, stop=7, length=2000))
h = range(-7, stop=7, length=2000)[2] - range(-7, stop=7, length=2000)[1]
@test isapprox(sum(A,1)*h, ones(4)'; atol=1e-4)

# test discretization approach vs. properly integrating..
n_marginal_grid = 1001
marginal_grid =  range(-3, stop=3, length=n_marginal_grid)
marginal_h = marginal_grid[2]-marginal_grid[1]

n_prior_grid = 1001
prior_grid = range(-6, stop=6, length=n_prior_grid)

normal_normal = DiscretizedNormalConvolutionProblem(Normal(0,0.5),
                    prior_grid,
                    marginal_grid)

marginal_dist = Normal(0,sqrt(1+0.5^2))
marginal_l = marginal_grid .- marginal_h/2
marginal_r = marginal_grid .+ marginal_h/2

f_marginal = [QuadGK.quadgk(x->pdf.(marginal_dist,x),marginal_l[i],marginal_r[i])[1] for i=1:n_marginal_grid];
f_marginal ./= sum(f_marginal);

@test isapprox(f_marginal, normal_normal.marginal; atol= 1e-7)

@test sum(normal_normal.prior) ≈ 1
@test sum(normal_normal.marginal) ≈ 1

# check if posterior mean calculated correctly (normal-normal conjugate we know this ofc)
post_mean = posterior_stats(normal_normal, μ -> μ, 2)
@test isapprox(post_mean[3], 2/1/(1+1/0.5^2))

# now check for indicator μ >= 0
# actually tolerances here are rather bad...esp. for ground truth should do a better job.

true_num = quadgk(β ->  pdf(Normal(0,0.5), β)*pdf(Normal(), 2-β) , 0,100)[1]
true_denom = quadgk(β ->  pdf(Normal(0,0.5), β)*pdf(Normal(), 2-β) , -100,100)[1]
true_ratio = true_num/true_denom

posterior_positive = posterior_stats(normal_normal, μ ->  (μ >= 0).*1, 2)
@test isapprox(true_num, posterior_positive[1]; atol=1e-3 )
@test isapprox(true_denom, posterior_positive[2]; atol=1e-3 )
@test isapprox(true_ratio, posterior_positive[3]; atol=1e-2)

# again in this case we know exactly what the value is..

post_Normal = Normal(2/1/(1+1/0.5^2), sqrt(1*0.5^2/(1+0.5^2)))
post_positive_true = ccdf(post_Normal, 0 )

@test isapprox(post_positive_true, true_ratio)
@test isapprox(post_positive_true, posterior_positive[3]; atol=1e-2)

## OK NOW let us check with the more exact object

normal_normal_exact = NormalConvolutionProblem(Normal(0,0.5), marginal_grid)
@test isapprox(f_marginal, normal_normal_exact.marginal; atol= 1e-9)

non_disr_post_positive = posterior_stats(normal_normal_exact, μ ->  (μ >= 0).*1, 2)
@test isapprox(post_positive_true, non_disr_post_positive[3]; atol=1e-8)


text:
"""
    AliasStyle

A trait that determines the aliasing behavior of a constructor or function,
for example whether or not a function or constructor might return an alias
of one of the inputs (i.e. the output shares memory with one of the inputs,
such that modifying the output also modifies the input or vice versa).

See also [`AllowAlias`](@ref) and [`NeverAlias`](@ref).
"""
abstract type AliasStyle end

"""
    AllowAlias

Singleton type used in a constructor or function indicating
that the constructor or function may return an alias of the input data when
possible, i.e. the output may share data with the input. For a constructor
`T(AllowAlias(), args...)`, this would act like `Base.convert(T, args...)`.

See also [`AliasStyle`](@ref) and [`NeverAlias`](@ref).
"""
struct AllowAlias <: AliasStyle end

"""
    NeverAlias

Singleton type used in a constructor or function indicating
that the constructor or function will never return an alias of the input data,
i.e. the output will never share data with one of the inputs.

See also [`AliasStyle`](@ref) and [`AllowAlias`](@ref).
"""
struct NeverAlias <: AliasStyle end


text:
# We're going to create a struct that executes a bunch of parallel branches and then
# concatenates the results back together
struct InceptionBlock{T <: Tuple}
    branches::T
end

InceptionBlock(x, y, args...) = InceptionBlock((x, y, args...))

# Flux likes the "@treelike" macro for doing backpropogation.
Flux.@functor InceptionBlock

function (block::InceptionBlock)(x)
    # Call branches
    results = map(f -> f(x), block.branches)
    return cat(results...; dims = 3)
end

const A_SIZE = (35, 35, 384)
const B_SIZE = (17, 17, 1024)
const C_SIZE = (8, 8, 1536)

# We're gonna be REALLY mean to type inference. Sorry Julia!
function stem()
    #println("Stem")

    return Chain(
        # Initial Stem
        Conv((3, 3), 3 => 32, relu; pad = 0, stride = 2),
        Conv((3, 3), 32 => 32, relu; pad = 0),
        Conv((3, 3), 32 => 64, relu; pad = 1),

        # First Split
        InceptionBlock(
            MaxPool((3,3); pad = 0, stride = 2),
            Conv((3,3), 64 => 96, relu; pad = 0, stride = 2),
        ),

        # Second Split
        InceptionBlock(
            Chain(
                Conv((1,1), 160 => 64, relu; pad = 0, stride = 1),
                Conv((3,3), 64 => 96, relu; pad = 0),
            ),
            Chain(
                Conv((1,1), 160 => 64, relu; pad = 0),
                Conv((7,1), 64 => 64, relu; pad = (3, 3, 0, 0)),
                Conv((1,7), 64 => 64, relu; pad = (0, 0, 3, 3)),
                Conv((3,3), 64 => 96, relu; pad = 0),
            ),
        ),

        # Final Split
        InceptionBlock(
            Conv((3,3), 192 => 192, relu; pad = 0, stride = 2),
            MaxPool((3,3); pad = 0, stride = 2),
        )
    )
end

function inception_a()
    #println("A Block")
    return InceptionBlock(
        Chain(
            x -> meanpool(x, (3,3); pad = 1, stride = 1),
            Conv((1,1), 384 => 96, relu; pad = 0)
        ),
        Conv((1,1), 384 => 96, relu; pad = 0),
        Chain(
            Conv((1,1), 384 => 64, relu; pad = 0),
            Conv((3,3), 64 => 96, relu; pad = 1)
        ),
        Chain(
            Conv((1,1), 384 => 64, relu; pad = 0),
            Conv((3,3), 64 => 96, relu; pad = 1),
            Conv((3,3), 96 => 96, relu; pad = 1)
        )
    )
end

function inception_b()
    #println("B Block")
    S = B_SIZE[3]
    return InceptionBlock(
        Chain(
            x -> meanpool(x, (3,3); pad = 1, stride = 1),
            Conv((1,1), S => 128)
        ),
        Conv((1,1), S => 384, relu),
        Chain(
            Conv((1,1), S => 192, relu),
            Conv((1,7), 192 => 224, relu; pad = (0, 0, 3, 3)),
            Conv((7,1), 224 => 256, relu; pad = (3, 3, 0, 0)),
        ),
        Chain(
            Conv((1,1), S => 192, relu),
            Conv((1,7), 192 => 192, relu; pad = (0, 0, 3, 3)),
            Conv((7,1), 192 => 224, relu; pad = (3, 3, 0, 0)),
            Conv((1,7), 224 => 224, relu; pad = (0, 0, 3, 3)),
            Conv((7,1), 224 => 256, relu; pad = (3, 3, 0, 0)),
        )
    )
end

function inception_c()
    #println("C Block")
    S = C_SIZE[3]
    return InceptionBlock(
        Chain(
            x -> meanpool(x, (3,3); pad = 1, stride = 1),
            Conv((1,1), S => 256, relu)
        ),
        Conv((1,1), S => 256, relu),
        Chain(
            Conv((1, 1), S => 384, relu),
            InceptionBlock(
                Conv((1, 3), 384 => 256, relu; pad = (0, 0, 1, 1)),
                Conv((3, 1), 384 => 256, relu; pad = (1, 1, 0, 0)),
            ),
        ),
        Chain(
            Conv((1,1), S => 384, relu),
            Conv((1,3), 384 => 448, relu; pad = (0, 0, 1, 1)),
            Conv((3,1), 448 => 512, relu; pad = (1, 1, 0, 0)),
            InceptionBlock(
                Conv((1,3), 512 => 256, relu; pad = (0, 0, 1, 1)),
                Conv((3,1), 512 => 256, relu; pad = (1, 1, 0, 0)),
            ),
        )
    )
end

function inception_ra(k, l, m, n)
    #println("A Reduction")
    S = A_SIZE[3]

    return InceptionBlock(
        x -> maxpool(x, (3,3); pad = 0, stride = 2),
        Conv((3,3), S => n, relu; pad = 0, stride = 2),
        Chain(
            Conv((1,1), S => k, relu),
            Conv((3,3), k => l, relu; pad = 1),
            Conv((3,3), l => m, relu; pad = 0, stride = 2)
        ),
    )
end

function inception_rb()
    #println("B Reduction")
    S = B_SIZE[3]

    return InceptionBlock(
        x -> maxpool(x, (3,3); pad = 0, stride = 2),

        Chain(
            Conv((1,1), S => 192, relu),
            Conv((3,3), 192 => 192, relu; pad = 0, stride = 2)
        ),

        Chain(
            Conv((1,1), S => 256, relu; pad = 0),
            Conv((1,7), 256 => 256, relu; pad = (0, 0, 3, 3)),
            Conv((7,1), 256 => 320, relu; pad = (3, 3, 0, 0)),
            Conv((3,3), 320 => 320, relu; pad = 0, stride = 2)
        ),
    )
end

# Legacy from when I didn't know the sizes of the intermediate layers.
#
# Removed the automatic size tracking because it took FOREVER to run the functions
# in the basic NNlib implementations.
mutable struct SizeTracker
    layers::Vector{Any}
    array::Any
end

function Base.push!(S::SizeTracker, f, call = true)
    # This is kind of mindbending ...
    #
    # We call the provided function to get another function that we append onto the
    # layers. Then, we call that generated function to get the new array size.
    if call
        push!(S.layers, f())
    else
        push!(S.layers, f)
    end
end

function inception_v4(x)
    layers = SizeTracker([], x)
    push!(layers, stem)
    for _ in 1:4
        push!(layers, inception_a)
    end
    push!(layers, () -> inception_ra(192, 224, 256, 384))

    for _ in 1:7
        push!(layers, inception_b)
    end
    push!(layers, inception_rb)

    for _ in 1:3
        push!(layers, inception_c)
    end

    #kernel_size = size.(Ref(layers.array), (1, 2))
    kernel_size = (C_SIZE[1], C_SIZE[2])
    push!(layers, x -> meanpool(x, kernel_size; pad = 0, stride = 1), false)
    # dropout

    push!(layers, x -> reshape(x, :, size(x,4)), false)
    push!(layers, Dense(1536, 1000, relu), false)
    push!(layers, x -> log.(max.(x, Float32(1e-7))), false),
    push!(layers, softmax, false)

    return UnstableChain(layers.layers)
end

function inception_v4_training(batchsize; kw...)
    x = rand(Float32, 299, 299, 3, batchsize)
    X = (x .- mean(x)) ./ std(x)

    Y = zeros(Float32, 1000, batchsize)
    random_labels!(Y)

    f = ForwardLoss(inception_v4(x), Flux.crossentropy)
    return Actualizer(f, X, Y; optimizer = nGraph.SGD(Float32(0.05)))
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
# use SGD for linear regression

using SGDOptim

function linreg_sgd(θ_g::Vector{Float64}, n::Int, σ::Float64)

    # prepare experimental data
    d = length(θ_g) - 1
    X = randn(d, n)
    y = vec(θ_g[1:d]'X) + θ_g[d+1] + σ * randn(n)

    # initialize solution
    θ_0 = zeros(d + 1)

    # optimize
    rmodel = riskmodel(AffinePred(d), SqrLoss())

    θ = sgd(rmodel, θ_0,
        minibatch_seq(X, y, 10),          # configure the way data are supplied
        reg = SqrL2Reg(1.0e-4),           # regularization
        lrate = t->1.0 / (100.0 + t),     # learing rate policy
        cbinterval = 5,                   # how frequently callback is invoked
        callback = gtcompare_trace(θ_g))  # the callback function

    # compare solution with initial guess
    println()
    @printf("Initial:  deviation = %.4e | avg.risk = %.4e\n",
        vecnorm(θ_0 - θ_g), value(rmodel, θ_0, X, y) / n)
    @printf("Solution: deviation = %.4e | avg.risk = %.4e\n",
        vecnorm(θ - θ_g), value(rmodel, θ, X, y) / n)
end

linreg_sgd([3.0, 5.0, 2.0], 10000, 0.1)


text:
"""
    GMM(K; max=1e+8, th=1e-4)
Gauss Mixture Model. This is useful when the data follow a Gaussian distribution.

# Parameters:
- `K`: number of class
- `max`: maximum number of repitition
- `th`: converge threshold

# Example
```jldoctest
julia> model = GMM(3)
GMM{3}(Float64[], Matrix{Float64}(undef, 0, 0), Array{Float64, 3}(undef, 0, 0, 0), 100000000, 0.0001)

julia> using HorseML.Clustering: fit!

julia> fit!(model, x)

julia> model(x) |> size
(100, 3)
```
"""
mutable struct GMM{K}
    π::Vector{Float64}
    μ::Matrix{Float64}
    Σ::Array{Float64, 3}
    max::Int64
    th::Float64
    GMM(K; max = 1e+8, th = 1e-4) = new{K}(Array{Float64}(undef, 0), Array{Float64}(undef, 0, 0), Array{Float64}(undef, 0, 0, 0), max, th)
end

function initialize(x, K, N)
    d = x[:, 1]
    ma = maximum(d)
    mi = minimum(d)
    interval = (ma-mi) / K
    γ = zeros(N, K)
    for i in 1 : K-1
        γ[findall(x -> mi <= x < mi+interval, d), i] .= 1
        mi += interval
    end
    γ[findall(x -> mi < x <= mi+interval, d), K] .= 1
    Mstep(x, γ, size(x)..., K)
end

function gauss(x, μ, σ)
    N, D = size(x)
    c1 = 1 / (2 * π)^(D/2)
    c2 = 1 / det(σ)^0.5
    c3 = x .- μ'
    c4 = c3 * σ'
    c5 = zeros(N)
    for d in 1 : D
        c5 += c4[:, d] .* c3[:, d]
    end
    p = @. c1 * c2 * exp(-c5 / 2)
    return p
end

function Estep(x, π, μ, Σ, N, K)
    y = zeros(N, K)
    for k in 1 : K
        y[:, k] = gauss(x, μ[k, :], Σ[k, :, :])
    end
    γ = zeros(N, K)
    for n in 1 : N
        wk = [π[k] * y[n, k] for k in 1 : K]
        γ[n, :] = wk / sum(wk)
    end
    return γ
end

function Mstep(x, γ, N, D, K)
    π = sum(γ, dims = 1)
    μ = zeros(K, D)
    for k in 1 : K
        μ[k, :] = γ[:, k]' * (x / sum(γ[:, k]))
    end
    Σ = zeros(K, D, D)
    for k in 1 : K
        for n in 1 : N
            wk = x .- μ[k, :]'
            wk = wk[n, :, :]
            Σ[k, :, :] += γ[n, k] * (wk * wk')
        end
        Σ[k, :, :] /= sum(γ[:, k])
    end
    return π, μ, Σ
end

function fit!(model::GMM{K}, x) where {K}
    coverge(x, th) = @. abs(x) < th
    function coverge(π, μ, Σ, th)
        πb = coverge.(π, th) == trues(size(π)...)
        μb = coverge.(μ, th) == trues(size(μ)...)
        Σb = coverge.(Σ, th) == trues(size(Σ)...)
        πb && μb && Σb
    end
    N, D = size(x)
    π, μ, Σ = initialize(x, K, N)
    for _ in 1 : model.max
        γ = Estep(x, π, μ, Σ, N, K)
        π_new, μ_new, Σ_new = Mstep(x, γ, N, D, K)
        if coverge(π_new - π, μ_new - μ, Σ_new - Σ, model.th)
            model.π, model.μ, model.Σ = vec(π_new), μ_new, Σ_new
            return
        end
        π, μ, Σ = π_new, μ_new, Σ_new
    end
    @warn "Not Converged!"
    model.π, model.μ, model.Σ = vec(π), μ, Σ
end

function (model::GMM{K})(x) where {K}
    N = size(x, 1)
    Estep(x, model.π, model.μ, model.Σ, N, K)
end


text:
using BigO
using Test
using Random
using  RecipesBase

@testset "BigO.jl" begin

    Random.seed!(42069);

    report = RunReport([sort, sort!], rand, 10:10:50, samples=3)
    @test report.names == ["sort", "sort!"]
    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), report)
    @test rec[1].plotattributes[:title] == "Runtimes"
    @test rec[1].plotattributes[:label] == ["sort" "sort!"]

    report1 = RunReport(sort, rand, 10:10:50, samples=3)
    @test report1.names == ["sort"]

    @test typeof(bigos(report)) == Dict{String, String}

    testfuns = [
        (x -> 0.3 + 0.7 * x + rand() / 5, "O(n)"),
        (x -> 7 + 2 * log(x) + rand() / 10, "O(log n)"),
        (x -> 27 + 13 * x ^ 1.7 + rand(), "O(n^p)"),
        (x -> 1000 + 5 * x + 1.8 ^ x + rand(), "O(k^n)"),
        (x -> 100 + 2.5 * x * log(5, x) + rand() / 3, "O(n log n)"),

    ]

    for (truf, name) in testfuns
        x = 50:20:400;
        y = truf.(x);
        guess = guessbigo(x, y)
        @test guess[1].name == name
    end

    @test (guessbigo(x -> 10) |> best) == "O(1)"

end


text:
module PartonDensity

include("parametrisation.jl")
include("qcdnum_interface.jl")
include("cross_section.jl")
include("zeus.jl")

end


text:
using Optim
using DataFrames

require("preprocess.jl")

dataset = :adult

df = readtable("adult.data")
# df = readtable("census-income.data")

# complete_cases!(df)

X = featurize(df[1:end-1])
y = class(df[end])

X = X[1:10000,:]
y = y[1:10000,:]

data = cat(2, X, y)

num_features = size(X, 2)

data = data'

theta = randn(num_features)

lambda = 0.001

function func(x::Vector)
  return loss(data, x, lambda)
end

function gradient!(x::Vector, storage::Vector)
  copy!(storage, full_grad(grad!, data, x, lambda))
end

res = optimize(func, gradient!, theta, method = :l_bfgs, show_trace=true, grtol=1e-5)

dataset_info = ["name" => "UCI", "theta_star" => res.minimum]

max_stars_repo_name: pcmoritz/slbfgs

max_issues_repo_name: pcmoritz/slbfgs

max_forks_repo_name: pcmoritz/slbfgs


text:
Pkg.add("ImageMagick")
Pkg.build("ImageMagick")

# Pkg.clone("GR")
# Pkg.build("GR")

Pkg.clone("https://github.com/JuliaPlots/PlotReferenceImages.jl.git")

# Pkg.clone("https://github.com/JuliaStats/KernelDensity.jl.git")

Pkg.clone("StatPlots")
Pkg.checkout("PlotUtils")

# Pkg.clone("Blink")
# Pkg.build("Blink")
# import Blink
# Blink.AtomShell.install()
# Pkg.add("Rsvg")
# Pkg.add("PlotlyJS")

# Pkg.checkout("RecipesBase")
# Pkg.clone("VisualRegressionTests")

# need this to use Conda
ENV["PYTHON"] = ""
Pkg.add("PyPlot")
Pkg.build("PyPlot")

# Pkg.add("InspectDR")

Pkg.test("Plots"; coverage=false)


text:
export MelodyOneHotEncoding
import ..encode_event, ..decode_event

using ..NoteSequences: MIN_MIDI_PITCH, MAX_MIDI_PITCH

"""
    MelodyOneHotEncoding <: Any

One-Hot encoder for `Melody` events.

## Fields
* `minpitch::Int`: Minimum melody pitch to allow for encoding.
* `maxpitch::Int`: Maximum melody pitch to allow for encoding.
* `num_classes::Int`: Total number of event classes.
"""
struct MelodyOneHotEncoding
    minpitch::Int
    maxpitch::Int
    num_classes::Int

    function MelodyOneHotEncoding(minpitch::Int, maxpitch::Int)
        maxpitch <= minpitch &&
            throw(error("max pitch should be greater than min pitch"))
        minpitch < MIN_MIDI_PITCH &&
            throw(error("min pitch is lesser than 0 (minimum allowed midi pitch)"))
        maxpitch > MAX_MIDI_PITCH + 1 &&
            throw(error("max pitch is greater than 128"))

        num_classes = maxpitch - minpitch + NUM_SPECIAL_MELODY_EVENTS
        new(minpitch, maxpitch, num_classes)
    end
end

function Base.getproperty(encoder::MelodyOneHotEncoding, sym::Symbol)
    if sym === :labels
        return 1:encoder.num_classes
    elseif sym === :defaultevent
        return MELODY_NO_EVENT
    else
        getfield(encoder, sym)
    end
end

"""
    encode_event(event::Int, encoder::MelodyOneHotEncoding)

Encodes a `Melody` event to its corresponding one hot index.
"""
function encode_event(event::Int, encoder::MelodyOneHotEncoding)
    0 <= event < encoder.minpitch &&
        throw(error("Melody event is less than min melody pitch"))
    event >= encoder.maxpitch &&
        throw(error("Melody event is greater than max melody pitch"))
    event < -NUM_SPECIAL_MELODY_EVENTS &&
        throw(error("Invalid melody event"))

    if event < 0
        index = event + NUM_SPECIAL_MELODY_EVENTS
    else
        index = event - encoder.minpitch + NUM_SPECIAL_MELODY_EVENTS
    end

    # Add one to convert to 1-based indexing (starts at 1, ends at num_classes)
    index + 1
end

"""
    decode_event(index::Int, encoder::MelodyOneHotEncoding)

Decodes a one hot index to its corresponding `Melody` event.
"""
function decode_event(index::Int, encoder::MelodyOneHotEncoding)
    if index < NUM_SPECIAL_MELODY_EVENTS + 1
        event = index - NUM_SPECIAL_MELODY_EVENTS
    else
        event = index - NUM_SPECIAL_MELODY_EVENTS + encoder.minpitch
    end

    # Subtract one to convert from 1-based indexing
    event - 1
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "ASL"
version = v"0.1.2"

# Collection of sources required to build ThinASLBuilder
sources = [
    ArchiveSource("http://netlib.org/ampl/solvers.tgz",
                  "65bb7dc72ce072bd4fa49578164a8430c03b458969cd4abe5ff6d8e9c09ad2dc"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/solvers/
head -17 asl.h > LICENSE.txt
mkdir -p ${includedir}
mkdir -p ${libdir}

all_load="--whole-archive"
noall_load="--no-whole-archive"
makefile="makefile.u"
cflags=""

if [[ -f $WORKSPACE/srcdir/asl-extra/arith.h.$target ]]; then
    cp $WORKSPACE/srcdir/asl-extra/arith.h.$target ./arith.h
fi

if [[ "${target}" == "arm-linux-musleabihf" ]]; then
    cp $WORKSPACE/srcdir/asl-extra/arith.h.arm-linux-gnueabihf ./arith.h
elif [[ $target == "i686-linux-musl" ]]; then
    cp $WORKSPACE/srcdir/asl-extra/arith.h.i686-linux-gnu ./arith.h
elif [[ "${target}" == "aarch64-linux-musl" ]]; then
    cp $WORKSPACE/srcdir/asl-extra/arith.h.aarch64-linux-gnu ./arith.h
elif [[ "${target}" == *-freebsd* ]]; then
    cp $WORKSPACE/srcdir/asl-extra/arith.h.x86_64-linux-gnu ./arith.h
    cflags="-D__XSI_VISIBLE=1"
elif [[ "${target}" == *-mingw* ]]; then
    makefile="$WORKSPACE/srcdir/asl-extra/makefile.mingw"
elif [[ "${target}" == *-apple-* ]]; then
    all_load="-all_load"
    noall_load="-noall_load"
fi

make -f $makefile CC="$CC" CFLAGS="-O -fPIC $cflags"
c++ -fPIC -shared -I$WORKSPACE/srcdir/asl-extra -I. $WORKSPACE/srcdir/asl-extra/aslinterface.cc -Wl,${all_load} amplsolver.a -Wl,${noall_load} -o libasl.${dlext}
cp libasl.${dlext} ${libdir}
cp *.h ${includedir}
cp *.hd ${includedir}
install_license LICENSE.txt
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libasl", :libasl)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# ╔═╡ f5780e90-247d-11eb-04d0-75c1e456dc12
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add(["Images", "ImageIO", "ImageMagick", "TestImages"])
	Pkg.add("ColorTypes")
	using Images, ColorTypes
end

# ╔═╡ 3e6381c0-247e-11eb-0328-f37226f82f93
md"""Checkout Color_Lab.jl or Color Block jl file for more things, as how to download a image and details as such"""

# ╔═╡ 3e49b830-247e-11eb-3ea2-2573e0f9bdac
White1 = "D:\\MyGit\\Colors\\Images\\3by3_White.jpg"

# ╔═╡ cd12a1a0-2486-11eb-19eb-97b70a9160b8
Black1 = "D:\\MyGit\\Colors\\Images\\3by3_Black.jpg"

# ╔═╡ ccf8ff20-2486-11eb-039e-c9de29594cdf
Green1 = "D:\\MyGit\\Colors\\Images\\3by3_Green.jpg"

# ╔═╡ cce01ff0-2486-11eb-1798-79c44f5b9931
Red1 = "D:\\MyGit\\Colors\\Images\\3by3_Red.jpg"

# ╔═╡ 124906a0-2488-11eb-2fb2-d56a1f33c05c
Blacknw1 = "D:\\MyGit\\Colors\\Images\\3by3_Black&white.jpg"

# ╔═╡ 3e2fa080-247e-11eb-0ad9-5dd8644886f9
White01 = load(White1)

# ╔═╡ 44956b40-2487-11eb-0c5c-9fe206234318
Black01 = load(Black1)

# ╔═╡ 4458141e-2487-11eb-3f74-9d8aaf6d688f
Green01 = load(Green1)

# ╔═╡ 4408bba0-2487-11eb-3bd9-9f8f5cc2b5c7
Red01 = load(Red1)

# ╔═╡ 48f39fd2-2488-11eb-203b-470163c3153d
BnW01 = load(Blacknw1)

# ╔═╡ 9758b490-2487-11eb-1c56-3d65bbe388be
channelview(White01), [2,3,1]

# ╔═╡ 973b4180-2487-11eb-2007-eba73bee14d6
channelview(Black01), [2,3,1]

# ╔═╡ 971fca40-2487-11eb-168c-63c1fb2b6014
channelview(Green01), [2,3,1]

# ╔═╡ c7cc3eae-247f-11eb-3e1b-a31f6e1ce000
channelview(Red01), [2,3,1]

# ╔═╡ c79bb8d0-247f-11eb-3107-7b86ccb91b04
channelview(BnW01), [2,3,1]

# ╔═╡ c78327c0-247f-11eb-0020-09c5507187ba
BnW011 = load(Blacknw1)

# ╔═╡ c752a1e0-247f-11eb-246a-ffb8fb0dbae5
typeof(BnW011)

# ╔═╡ acaf1250-248f-11eb-347e-39cd376a6085
channelview(BnW01), [2,3,1]

# ╔═╡ 59fa9500-255e-11eb-0951-bbdbeb5c013d
md"""below will display a pixel of an image"""

# ╔═╡ 2705d7e0-255e-11eb-2f48-b7a48c280ddf
CD = BnW01[1,2]

# ╔═╡ 74ec9f70-255e-11eb-240f-cd14eca7c3d2
dump(c24)

# ╔═╡ 74d2aed0-255e-11eb-112d-256463b71abe
c24.color

# ╔═╡ 26d048f0-255e-11eb-308d-69e7a209938e
RGB24(RGB(1, 0, 0)).color

# ╔═╡ 32f91e80-255f-11eb-022f-d94070720a5f
s,d = size(BnW01)

# ╔═╡ 3233aba0-255f-11eb-2605-bfbe991e6337
begin
	 W01[1:s, 1:d] .= RGB24(BnW01[1:s, 1:d]).color
end

# ╔═╡ 6b882f20-255f-11eb-32a7-b7763bb7cb17
RGB24(BnW01).color

# ╔═╡ 6b6e3e80-255f-11eb-2cc9-cb60bafc3512

# ╔═╡ 6b5649b0-255f-11eb-10ba-c766f8e6c4cd

# ╔═╡ 26b51fd0-255e-11eb-0be7-39c12681a6bb
begin
	Dhwaj2 = copy(Tricolor)
	Dhwaj2[70:90, 70:110] .= RGB(0, 1, 1)
	Dhwaj2
end

# ╔═╡ 2698e540-255e-11eb-2306-6996fa3150f8

# ╔═╡ 26800612-255e-11eb-362e-c7b1e43fb2a6

# ╔═╡ ac79d180-248f-11eb-319b-0d04957afbf2

# ╔═╡ e4d06160-255d-11eb-3f63-490f5332f468

# ╔═╡ e4b5ad70-255d-11eb-1648-7730c93b2e33

# ╔═╡ e49a844e-255d-11eb-3135-79f1bcb5ad8a

# ╔═╡ e4817e10-255d-11eb-007e-41de1bb4d68c

# ╔═╡ e263f540-255d-11eb-15c1-438852517e77

# ╔═╡ e248cc20-255d-11eb-3cd7-dbfa1460d37f

# ╔═╡ e22edb80-255d-11eb-09fe-aff782cd12bf

# ╔═╡ e22eb470-255d-11eb-190b-bf34def0f996

# ╔═╡ e1ffdc40-255d-11eb-0c3b-7198c874ff3a

# ╔═╡ 9c361c80-248e-11eb-02fc-33e6ab8801bf
RGB(1, 0, 0)

# ╔═╡ 9c1b1a6e-248e-11eb-395c-db91bf6eee0a
[RGB(1, 0, 0) RGB(0, 1, 0);
 RGB(0, 0, 1) RGB(0, 0, 0)]

# ╔═╡ 9c00b4a0-248e-11eb-382f-e9b74cb77e12

# ╔═╡ 9be23020-248e-11eb-068a-3ba581ce11b6

# ╔═╡ 3de2e010-247e-11eb-38d8-6b8bffd99436
md"""http://blog.robindeits.com/2017/02/27/reading-raw-memory-with-julia/ """

# ╔═╡ 3dc8ef70-247e-11eb-376b-05ad5147ca85
md"""https://craftofcoding.wordpress.com/2017/03/25/plotting-image-histograms-with-julia/"""

# ╔═╡ 3db01042-247e-11eb-2f66-dfc3c4e9e7bc
md"""In julia, you can just say convert(Gray, rgb)"""

# ╔═╡ 3d950e30-247e-11eb-2fa2-f7f1c1bc5c7f
md"""https://docs.julialang.org/en/v1/base/numbers/"""

# ╔═╡ 3d7af680-247e-11eb-29ea-231b8b849389
md"""https://github.com/JuliaGraphics/ColorTypes.jl"""

# ╔═╡ 3d0c0810-247e-11eb-040c-77c7433ccedd

# ╔═╡ 3c636c00-247e-11eb-369c-cde034bcdd36

# ╔═╡ 75069010-255e-11eb-0bbd-f9f4a97cc8f6
c24 = RGB24(CD)

# ╔═╡ 584766c0-255e-11eb-03be-b17de3aacabc
c24 = RGB24(c); dump(c24)

# ╔═╡ Cell order:
# ╠═f5780e90-247d-11eb-04d0-75c1e456dc12
# ╟─3e6381c0-247e-11eb-0328-f37226f82f93
# ╠═3e49b830-247e-11eb-3ea2-2573e0f9bdac
# ╠═cd12a1a0-2486-11eb-19eb-97b70a9160b8
# ╠═ccf8ff20-2486-11eb-039e-c9de29594cdf
# ╠═cce01ff0-2486-11eb-1798-79c44f5b9931
# ╠═124906a0-2488-11eb-2fb2-d56a1f33c05c
# ╠═3e2fa080-247e-11eb-0ad9-5dd8644886f9
# ╠═44956b40-2487-11eb-0c5c-9fe206234318
# ╠═4458141e-2487-11eb-3f74-9d8aaf6d688f
# ╠═4408bba0-2487-11eb-3bd9-9f8f5cc2b5c7
# ╠═48f39fd2-2488-11eb-203b-470163c3153d
# ╠═9758b490-2487-11eb-1c56-3d65bbe388be
# ╠═973b4180-2487-11eb-2007-eba73bee14d6
# ╠═971fca40-2487-11eb-168c-63c1fb2b6014
# ╠═c7cc3eae-247f-11eb-3e1b-a31f6e1ce000
# ╠═c79bb8d0-247f-11eb-3107-7b86ccb91b04
# ╠═c78327c0-247f-11eb-0020-09c5507187ba
# ╠═c752a1e0-247f-11eb-246a-ffb8fb0dbae5
# ╠═acaf1250-248f-11eb-347e-39cd376a6085
# ╟─59fa9500-255e-11eb-0951-bbdbeb5c013d
# ╠═2705d7e0-255e-11eb-2f48-b7a48c280ddf
# ╠═75069010-255e-11eb-0bbd-f9f4a97cc8f6
# ╠═74ec9f70-255e-11eb-240f-cd14eca7c3d2
# ╠═74d2aed0-255e-11eb-112d-256463b71abe
# ╠═584766c0-255e-11eb-03be-b17de3aacabc
# ╠═26d048f0-255e-11eb-308d-69e7a209938e
# ╠═32f91e80-255f-11eb-022f-d94070720a5f
# ╠═3233aba0-255f-11eb-2605-bfbe991e6337
# ╠═6b882f20-255f-11eb-32a7-b7763bb7cb17
# ╠═6b6e3e80-255f-11eb-2cc9-cb60bafc3512
# ╠═6b5649b0-255f-11eb-10ba-c766f8e6c4cd
# ╠═26b51fd0-255e-11eb-0be7-39c12681a6bb
# ╠═2698e540-255e-11eb-2306-6996fa3150f8
# ╠═26800612-255e-11eb-362e-c7b1e43fb2a6
# ╠═ac79d180-248f-11eb-319b-0d04957afbf2
# ╠═e4d06160-255d-11eb-3f63-490f5332f468
# ╠═e4b5ad70-255d-11eb-1648-7730c93b2e33
# ╠═e49a844e-255d-11eb-3135-79f1bcb5ad8a
# ╠═e4817e10-255d-11eb-007e-41de1bb4d68c
# ╠═e263f540-255d-11eb-15c1-438852517e77
# ╠═e248cc20-255d-11eb-3cd7-dbfa1460d37f
# ╠═e22edb80-255d-11eb-09fe-aff782cd12bf
# ╠═e22eb470-255d-11eb-190b-bf34def0f996
# ╠═e1ffdc40-255d-11eb-0c3b-7198c874ff3a
# ╠═9c361c80-248e-11eb-02fc-33e6ab8801bf
# ╠═9c1b1a6e-248e-11eb-395c-db91bf6eee0a
# ╠═9c00b4a0-248e-11eb-382f-e9b74cb77e12
# ╠═9be23020-248e-11eb-068a-3ba581ce11b6
# ╟─3de2e010-247e-11eb-38d8-6b8bffd99436
# ╟─3dc8ef70-247e-11eb-376b-05ad5147ca85
# ╟─3db01042-247e-11eb-2f66-dfc3c4e9e7bc
# ╟─3d950e30-247e-11eb-2fa2-f7f1c1bc5c7f
# ╟─3d7af680-247e-11eb-29ea-231b8b849389
# ╠═3d0c0810-247e-11eb-040c-77c7433ccedd
# ╠═3c636c00-247e-11eb-369c-cde034bcdd36


text:
include("steps.jl")

config = Dict{String, Any}()
for d in readdir()
    if isdir(d) && isfile(joinpath(d, "config.jl"))
        include(joinpath(d, "config.jl"))
    end
end

function main()
    setup_threading()

    run_step_on_cases(firstrun, config)
    run_step_on_cases(optimal_damping, config)
    run_step_on_cases(scf, config)
end


text:
"""
Calculate distances, angles, speeds and points for aviation navigation purposes.

All angles are in radians, distances are standard in meters, and speeds are
standard in m/s.

These methods are not to be used for operational purposes.

Implemented Functions:
* distance
* bearing
* final_bearing
* midpoint
* intermediate_point
* destination_point
* intersection_point
* along_track_distance
* Vground
* head_wind
* cross_wind
* normalize

Implemented Types:
* Point(ϕ, λ)

Implemented constants:
* Rₑ    Radius Earth in [m]

Based on
source: https://www.movable-type.co.uk/scripts/latlong.html
source: http://edwilliams.org/avform.htm
"""
module Navigation

import Base.-, Base.Math.rad2deg, Base.Math.deg2rad

export Point, distance, bearing, final_bearing, midpoint,
intermediate_point, Vground, destination_point, intersection_point,
cross_track_distance, along_track_distance, head_wind, cross_wind

"""
Radius Earth [m]
https://en.wikipedia.org/wiki/Earth_radius
"""
const Rₑ = 6371008.7714

"Point type with latitude ϕ [rad] and longitude λ [rad]"
struct Point{T<:Float64}
    ϕ::T
    λ::T
end

(-)(x::Point, y::Point) = Point(x.ϕ - y.ϕ, x.λ - y.λ)
(rad2deg)(x::Point) = Point(rad2deg(x.ϕ), rad2deg(x.λ))
(deg2rad)(x::Point) = Point(deg2rad(x.ϕ), deg2rad(x.λ))

"""
"""
# function normalize_angle(angle::Float64)
#     return (angle + 2.0*π)
# end

"""
distance(pos₁::Point, pos₂::Point, radius::Float64=Rₑ)

Distance (=[Rₑ]) on a sphere calculated using the haversine formula
The haversine gives also good estimations at short distances

source: https://www.movable-type.co.uk/scripts/latlong.html
"""
function distance(pos₁::Point, pos₂::Point, radius::Float64=Rₑ)
    Δpos = pos₂ - pos₁
    a = sin(Δpos.ϕ/2.0)^2 + cos(pos₁.ϕ)*cos(pos₂.ϕ)*sin(Δpos.λ/2.0)^2
    c = 2.0 * atan(√a, √(1.0 - a))
    return radius * c
end

"""
bearing(pos₁::Point, pos₂::Point)

Initial bearing [rad] of the great circle line between the positions 1
(Point) and 2 (Point) on a sphere

source: https://www.movable-type.co.uk/scripts/latlong.html
"""
function bearing(pos₁::Point, pos₂::Point)
    Δpos = pos₂ - pos₁
    return mod2pi(atan(sin(Δpos.λ)*cos(pos₂.ϕ)
                ,cos(pos₁.ϕ)*sin(pos₂.ϕ) - sin(pos₁.ϕ)*cos(pos₂.ϕ)*cos(Δpos.λ)))
end

"""
final_bearing(pos₁::Point, pos₂::Point)

Final bearing [rad] of the great circle line between the positions 1
(Point) and 2 (Point) on a sphere

source: https://www.movable-type.co.uk/scripts/latlong.html
"""
function final_bearing(pos₁::Point, pos₂::Point)
    return mod2pi(bearing(pos₂, pos₁) + π)
end

"""
midpoint(pos₁::Point, pos₂::Point)

The half-way point (Point) on a great circle path between two points (Point)

source: https://www.movable-type.co.uk/scripts/latlong.html
"""
function midpoint(pos₁::Point, pos₂::Point)
    Δpos = pos₂ - pos₁
    Bx = cos(pos₂.ϕ) * cos(Δpos.λ)
    By = cos(pos₂.ϕ) * sin(Δpos.λ)
    ϕₘ = atan(sin(pos₁.ϕ) + sin(pos₂.ϕ), √((cos(pos₁.ϕ) + Bx)^2 + By^2))
    λₘ = pos₁.λ + atan(By, cos(pos₁.ϕ) + Bx)
    return Point(ϕₘ, normalize(λₘ, -π, Float64(π)))
end

"""
intermediate_point(pos₁::Point, pos₂::Point, fraction::Float64 = 0.5)

An intermediate point (Point) at any fraction along the great circle path
between two points 1 (Point) and 2 (Point). The fraction along the great circle
route is such that fraction = 0.0 is at point 1 and fraction 1.0 is at point 2.

source: https://www.movable-type.co.uk/scripts/latlong.html
"""
function intermediate_point(pos₁::Point, pos₂::Point, fraction::Float64 = 0.5)
    δ = distance(pos₁, pos₂) / Rₑ
    a = sin((1.0 - fraction) * δ) / sin(δ)
    b = sin(fraction * δ) / sin(δ)
    x = a * cos(pos₁.ϕ) * cos(pos₁.λ) + b * cos(pos₂.ϕ) * cos(pos₂.λ)
    y = a * cos(pos₁.ϕ) * sin(pos₁.λ) + b * cos(pos₂.ϕ) * sin(pos₂.λ)
    z = a * sin(pos₁.ϕ) + b * sin(pos₂.ϕ)
    ϕᵢ = atan(z, √(x*x + y*y))
    λᵢ = atan(y, x)
    return Point(ϕᵢ, λᵢ)
end

"""
destination_point(start_pos::Point, distance::Float64, bearing::Float64
, radius::Float64=Rₑ)

Given a start point (Point), initial bearing [rad] (clockwise from North),
distance (=[Rₑ]), and the earth radius (Rₑ) calculate the destination point
(Point) travelling along a (shortest distance) great circle arc.
(The distance must use the same radius as the radius of the earth.)

source: https://www.movable-type.co.uk/scripts/latlong.html
"""
function destination_point(start_pos::Point, distance::Float64, bearing::Float64,
    radius::Float64=Rₑ)
    δ = distance / radius  # Angular distance
    ϕ₂ = asin(sin(start_pos.ϕ) * cos(δ) + cos(start_pos.ϕ) * sin(δ) * cos(bearing))
    λ₂ = start_pos.λ + atan(sin(bearing) * sin(δ) * cos(start_pos.ϕ),
    cos(δ) - sin(start_pos.ϕ) * sin(ϕ₂))
    return Point(ϕ₂, λ₂)
end

"""
intersection_point(pos₁::Point, pos₂::Point, bearing₁₃::Float64, bearing₂₃::Float64)

The intersection point pos₃ (Point) of two great circle paths given two start
points (Point) pos₁ and pos₂ and two bearings [rad] from pos₁ to pos₃, and from
pos₂ to pos₃.

Under certain circumstances the results can be an infinite or ambiguous solution.

source: http://edwilliams.org/avform.htm
"""
function intersection_point(pos₁::Point, pos₂::Point, bearing₁₃::Float64,
    bearing₂₃::Float64)
    Δpos = pos₂ - pos₁
    δ₁₂ = 2asin(√(sin(Δpos.ϕ / 2)^2 + cos(pos₁.ϕ) * cos(pos₂.ϕ) * sin(Δpos.λ/2)^2))
    θₐ = acos((sin(pos₂.ϕ) - sin(pos₁.ϕ) * cos(δ₁₂)) / (sin(δ₁₂) * cos(pos₁.ϕ)))
    θᵦ = acos((sin(pos₁.ϕ) - sin(pos₂.ϕ) * cos(δ₁₂)) / (sin(δ₁₂) * cos(pos₂.ϕ)))
    local θ₁₂, θ₂₁
    if sin(Δpos.λ) > 0.0
        θ₁₂ = θₐ
        θ₂₁ = 2π - θᵦ
    else
        θ₁₂ = 2π - θₐ
        θ₂₁ = θᵦ
    end
    α₁ = (bearing₁₃ - θ₁₂ + π)%(2π) - π
    α₂ = (θ₂₁ - bearing₂₃ + π)%(2π) - π

    if sin(α₁) ≈ 0.0 && sin(α₂) ≈ 0.0
        return Point(Inf, Inf)  # infinity of intersections
    elseif sin(α₁) * sin(α₂) < 0.0
        return Point(Inf, Inf)  # intersection ambiguous
    else
        α₃ = acos(-cos(α₁) * cos(α₂) + sin(α₁) * sin(α₂) * cos(δ₁₂))
        δ₁₃ = atan(sin(δ₁₂) * sin(α₁) * sin(α₂), cos(α₂) + cos(α₁) * cos(α₃))
        ϕ₃ = asin(sin(pos₁.ϕ) * cos(δ₁₃) + cos(pos₁.ϕ) * sin(δ₁₃) * cos(bearing₁₃))
        Δλ₁₃ = atan(sin(bearing₁₃) * sin(δ₁₃) * cos(pos₁.ϕ), cos(δ₁₃)
        - sin(pos₁.ϕ) * sin(ϕ₃))
        λ₃ = pos₁.λ + Δλ₁₃
        return Point(ϕ₃, λ₃)
    end
end

"""
cross_track_distance(pos₁::Point, pos₂::Point, pos₃::Point, radius::Float64=Rₑ)

Cross-track distance from a point pos₃ (Point) to a great circle path defined
by the points pos₁ and pos₂ (Point).

cross_track_distance(pos₁::Point, bearing::Float64, pos₃::Point, radius::Float64=Rₑ)

Cross-track distance from a point pos₃ (Point) to a great circle path defined
by the point pos₁ and bearing [rad].

cross_track_distance(∠distance₁₃::Float64, bearing₁₂::Float64,
    bearing₁₃::Float64, radius::Float64=Rₑ)

Cross-track distance from a point pos₃ to a great circle path. The distance is
calculated using the angular distance between points pos₁ and pos₃ (Point),
the bearing [rad] between points pos₁ and pos₃ (Point), the bearing [rad]
between points pos₁ and pos₂ (Point), and the radius of the earth.

source: https://www.movable-type.co.uk/scripts/latlong.html
"""
function cross_track_distance(pos₁::Point, pos₂::Point, pos₃::Point,
    radius::Float64=Rₑ)
    ∠distance₁₃ = distance(pos₁, pos₃ , radius) / radius
    bearing₁₃ = bearing(pos₁, pos₃)
    bearing₁₂ = bearing(pos₁, pos₂)
    return cross_track_distance(∠distance₁₃, bearing₁₂, bearing₁₃, radius)
end

function cross_track_distance(pos₁::Point, bearing::Float64, pos₃::Point,
    radius::Float64=Rₑ)
    pos₂ = destination_point(pos₁, radius, bearing, radius)
    return cross_track_distance(pos₁, pos₂, pos₃, radius)
end

function cross_track_distance(∠distance₁₃::Float64, bearing₁₂::Float64,
    bearing₁₃::Float64, radius::Float64=Rₑ)
    return asin(sin(∠distance₁₃) * sin(bearing₁₃ - bearing₁₂)) * radius
end

"""
along_track_distance(pos₁::Point, pos₂::Point, pos₃::Point, radius::Float64=Rₑ)

The along-track distance from the start point pos₁ (Point) to the closest point
on the great circle path (defined by points pos₁ and pos₂ (Point)) to the point
pos₃ (Point). The radius of the earth can also be given.

along_track_distance(∠distance₁₃::Float64, cross_track_∠distance::Float64, radius::Float64=Rₑ)

source: https://www.movable-type.co.uk/scripts/latlong.html
"""
function along_track_distance(pos₁::Point, pos₂::Point, pos₃::Point,
    radius::Float64=Rₑ)
    ∠distance₁₃ = distance(pos₁, pos₃, radius)
    cross_track_∠distance = cross_track_distance(pos₁, pos₂, pos₃, radius) / radius
    return along_track_distance(∠distance₁₃, cross_track_∠distance, radius)
end

function along_track_distance(∠distance₁₃::Float64,
    cross_track_∠distance::Float64, radius::Float64=Rₑ)
    return acos(cos(∠distance₁₃) / cos(cross_track_∠distance)) * radius
end

"""
Vground(Vtas::Float64, Vwind::Float64, ∠wind::Float64, course::Float64)

Ground speed (=[Vtas]) given the course [rad], airspeed,
wind speed (=[Vtas]) and wind direction [rad]. The speeds need to use the same
speed unit.

source: http://edwilliams.org/avform.htm
"""
function Vground(Vtas::Float64, Vwind::Float64, ∠wind::Float64, course::Float64)
    swc = (Vwind / Vtas) * sin(∠wind - course)
    return Vtas * √(1-swc*swc) - Vwind * cos(∠wind - course)
end

"""
head_wind(Vwind::Float64, ∠wind::Float64, course::Float64)

Head-wind component. A tail-wind is a negative head-wind.
"""
function head_wind(Vwind::Float64, ∠wind::Float64, course::Float64)
    return Vwind*cos(∠wind - course)
end

"""
cross_wind(Vwind::Float64, ∠wind::Float64, course::Float64)

Cross-wind component. A positive cross-wind component indicates a wind from the
right.
"""
function cross_wind(Vwind::Float64, ∠wind::Float64, course::Float64)
    return Vwind*sin(∠wind - course)
end

"""
normalize(value::Float64, lower::Float64 = 0.0, upper::Float64 = 360.0)

Normalize a value to stay within the lower and upper limit:
normalize(370.0, lower=0.0, upper=360.0)
10.0
"""
function normalize(value::Float64, lower::Float64 = 0.0, upper::Float64 = 360.0)
    result = value
    if result > upper || value == lower
        value = lower + abs(value + upper) % (abs(lower) + abs(upper))
    end
    if result < lower || value == upper
        value = upper - abs(value - lower) % (abs(lower) + abs(upper))
    end
    value == upper ? result = lower : result = value
    return result
end

#TODO Closest point to the poles
#TODO Crossing parallels
#TODO Rhumb lines

end # module


text:
scene3 = AbstractPlotting.vbox(scene1, scene2)
boundingbox(scene2)

scene3
scene3.children[2].plots[2].transformation.model[]
scene = Scene(resolution = (500, 500))
vx = -1:0.1:1;
vy = -1:0.1:1;

f(x, y) = (sin(x*10) + cos(y*10)) / 4
psurf = surface(vx, vy, f)

pos = lift_node(psurf[:x], psurf[:y], psurf[:z]) do x, y, z
    vec(Point3f0.(x, y', z .+ 0.5))
end
pscat = scatter(pos)
plines = lines(view(pos, 1:2:length(pos)))
center!(scene)
@theme theme = begin
    markersize = to_markersize2d(0.01)
    strokecolor = to_color(:white)
    strokewidth = to_float(0.01)
end
# this pushes all the values from theme to the plot
push!(pscat, theme)
pscat[:glow_color] = to_node(RGBA(0, 0, 0, 0.4), x->to_color((), x))
# apply it to the scene
custom_theme(scene)
# From now everything will be plotted with new theme
psurf = surface(vx, 1:0.1:2, psurf[:z])
center!(scene)

#cell
scene = Scene(resolution = (500, 500))

x = map([:dot, :dash, :dashdot], [2, 3, 4]) do ls, lw
    linesegment(range(1, stop=5, length=100), rand(100), rand(100), linestyle = ls, linewidth = lw)
end
push!(x, scatter(range(1, stop=5, length=100), rand(100), rand(100)))
center!(scene)
l = Makie.legend(x, ["attribute $i" for i in 1:4])
l[:position] = (0, 1)
l[:backgroundcolor] = RGBA(0.95, 0.95, 0.95)
l[:strokecolor] = RGB(0.8, 0.8, 0.8)
l[:gap] = 30
l[:textsize] = 19
l[:linepattern] = Point2f0[(0,-0.2), (0.5, 0.2), (0.5, 0.2), (1.0, -0.2)]
l[:scatterpattern] = decompose(Point2f0, Circle(Point2f0(0.5, 0), 0.3f0), 9)
l[:markersize] = 2f0
scene

#cell
scene = Scene(resolution = (500, 500))
cmap = collect(range(to_color(:red), stop=to_color(:blue), length=20))
l = Makie.legend(cmap, 1:4)
l[:position] = (1.0,1.0)
l[:textcolor] = :blue
l[:strokecolor] = :black
l[:strokewidth] = 1
l[:textsize] = 15
l[:textgap] = 5
scene

#cell
using Makie, GeometryTypes, ColorTypes
scene = Scene();
scatter([Point2f0(1.0f0,1.0f0),Point2f0(1.0f0,0.0f0)])
center!(scene);
text_overlay!(scene, "test", position = Point2f0(1.0f0,1.0f0), textsize=200,color= RGBA(0.0f0,0.0f0,0.0f0,1.0f0))
text_overlay!(scene, "test", position = Point2f0(1.0f0,0.0f0), textsize=200,color= RGBA(0.0f0,0.0f0,0.0f0,1.0f0))

scene = Scene();
scatter([Point2f0(1.0f0,1.0f0),Point2f0(1.0f0,0.0f0)])
center!(scene);

text_overlay!(scene,:scatter, "test", "test", textsize=200,color= RGBA(0.0f0,0.0f0,0.0f0,1.0f0))

scene = Scene();
scatter([Point2f0(1.0f0,1.0f0),Point2f0(1.0f0,0.0f0)])
center!(scene);
text_overlay!(scene, :scatter, 1=>"test1", 2=>"test2", textsize=200,color= RGBA(0.0f0,0.0f0,0.0f0,1.0f0))

#cell

# needs to be in a function for ∇ˢf to be fast and inferable
function test(scene)
    n = 20
    f   = (x,y,z) -> x*exp(cos(y)*z)
    ∇f  = (x,y,z) -> Point3f0(exp(cos(y)*z), -sin(y)*z*x*exp(cos(y)*z), x*cos(y)*exp(cos(y)*z))
    ∇ˢf = (x,y,z) -> ∇f(x,y,z) - Point3f0(x,y,z)*dot(Point3f0(x,y,z), ∇f(x,y,z))
    θ = [0;(0.5:n-0.5)/n;1]
    φ = [(0:2n-2)*2/(2n-1);2]
    x = [cospi(φ)*sinpi(θ) for θ in θ, φ in φ]
    y = [sinpi(φ)*sinpi(θ) for θ in θ, φ in φ]
    z = [cospi(θ) for θ in θ, φ in φ]

    pts = vec(Point3f0.(x, y, z))
    lns = Makie.streamlines!(scene, pts, ∇ˢf)
    # those can be changed interactively:
    lns[:color] = :black
    lns[:h] = 0.06
    lns[:linewidth] = 1.0
    lns
end

using Makie
main = Scene()
cam3d!(main)
main
plots = [
    heatmap(0..1, 0..1, rand(100, 100)),
    meshscatter(rand(10), rand(10), rand(10)),
    scatter(rand(10), rand(10)),
    mesh(Makie.loadasset("cat.obj")),
    volume(0..1, 0..1, 0..1, rand(32, 32, 32)),
]

for p in plots
    push!(main, p)
    translate!(p, rand()*3, rand()*3, rand()*3)
end

using Makie, GeometryTypes
scene = Scene()
ui = Scene(scene, lift(x-> IRect(0, 0, widths(x)[1], 100), pixelarea(scene)))
plots = Scene(scene, lift(x-> IRect(0, 100, widths(x) .- Vec(0, 100)), pixelarea(scene)))
campixel!(ui)
s1 = slider!(ui, 1:10, raw = true)[end]
s2 = slider!(ui, 1:10, raw = true)[end]
s3 = slider!(ui, range(0, stop=1, length=100), raw = true)[end]
AbstractPlotting.vbox(s1, s2, s3)
heatmap!(plots, rand(100, 100))
scene

@inline function collide(p, bounds)
    mini = p .<= minimum(bounds)
    any(mini) && return true, normalize(Point3f0(mini))
    maxi = p .>= maximum(bounds)
    any(maxi) && return true, normalize(-Point3f0(maxi))
    false, p
end

@inline function particle_inner(posj, posi, vel)
    d = posj - posi
    distsq = dot(d, d) + 1f0
    vel .+ (reverse(d)/distsq)
end

function solve_particles!(
        positions::AbstractVector{P}, velocity::AbstractVector, bounds, dt::T = T(0.01)
    ) where P <: Point{N, T} where {N, T}
    @inbounds for i in eachindex(positions)
        vel = velocity[i]
        posi = positions[i]
        for j in eachindex(positions)
            posj = positions[j]
            d = posj .- posi
            distsq = dot(d, d) + T(1)
            vel = vel .+ (d ./ distsq)
            any(x-> abs(x) > T(0.8), vel) && break # restrict velocity
        end
        col, normal = collide(posi, bounds)
        if col
            vel = -2f0 * (dot(vel, normal) * normal + vel)
        end
        velocity[i] = vel
        positions[i] = posi .+ dt*vel
    end
    return
end

startpositions(N::Integer, radius::T, n) where T = startpositions(Val(N), radius, n)
function startpositions(::Val{N}, radius::T, n) where {N, T}
    sphere = HyperSphere(Point{N, T}(0), T(radius))
    n = N == 3 ? floor(Int, sqrt(n)) : n # n must be n^2 for 3D Sphere
    pos = decompose(Point{N, T}, sphere, n)
    map!(pos, pos) do p
        p .+ ((rand(Point{N, T}) .- 0.5) .* radius * T(0.1))
    end
end
using GeometryTypes
FRect3D = HyperRectangle
bounds = FRect3D(Vec3f0(-1), Vec3f0(2))
N = 3; T = Float32; n = 10^5
positions = startpositions(N, T(0.5), n)
velocities = rand(Point3f0, length(positions))

using Makie
scene = scatter(positions, markersize = 0.006, color = norm.(velocities))
particles = scene[end]
linesegments!(scene, bounds)
display(Makie.global_gl_screen(), scene)

@async while isopen(scene)
    @time solve_particles!(positions, velocities, bounds)
    particles[1] = positions
    sleep(1/60)
end
# keep_runnin[] = true
scene
velocities
solve_particles2!(positions, velocities, bounds, 0.1f0)

Profile.print()
any.(x-> abs(x) > T(0.8), velocities)

using Makie
import RDatasets
singers = RDatasets.dataset("lattice","singer")

keys = unique(singers[:VoicePart])
map(typeof, singers[:VoicePart])
findfirst
first(keys) == first(singers[:VoicePart])
findfirst(keys, first(singers[:VoicePart]))

x = map(x-> findfirst(keys, x), singers[:VoicePart])

box(x, singers[:Height])

scene = Scene()
cam2d!(scene)
axis2d!(scene, ((-3, 4), (3, 8)))
center!(scene)
scene
@which AbstractPlotting.plots_from_camera(scene)
@which boundingbox()