text:
"""
    getapeanalysis(file2d,file3d,outputfile,initial_timeindex,final_timeindex,partnumber,outputInterval,position,FloatType::Type=Float64)
"""
function getapeanalysis_with_condensate_loading(file2d,file3d,outputfile,initial_timeindex,final_timeindex,partnumber,outputInterval,position,FloatType::Type=Float64;tropopause_height = 15000)

    @info "Starting APE analysis routine"
    day           = 86400
    sst           = 300
    dt            = outputInterval

    dayLength     = 60*60*24÷outputInterval; #How many data points make one day
    #final_time    = floor(Int,final_time_days*dayLength)
    #initial_time  = floor(Int,initial_time_days*dayLength)+1

    Pref          = 1000*1e2                 #Pa
    Ts            = sst                      #Sea surface temperature
    qs            = 25.7*1e-3
    Tvs           = Ts*(1+epsilon*qs)
    c1            = (Dryair.R/Dryair.cp)

    #iterator_time_2d    = initial_timeindex*2-1:2:final_timeindex*2
    #iterator_time_3d    = initial_timeindex:1:final_timeindex

    smooth_x      = smooth_y = 15 #it was 11
    smooth_time   = floor(Int,dayLength*5)+1

    if position == 1
        iterator_time_2d    = initial_timeindex*2-1:2:(final_timeindex+div(smooth_time-1,2))*2
        iterator_time_3d    = initial_timeindex:1:final_timeindex+div(smooth_time-1,2)
    elseif position == 2
        iterator_time_2d    = (initial_timeindex-div(smooth_time-1,2))*2-1:2:(final_timeindex+div(smooth_time-1,2))*2
        iterator_time_3d    = initial_timeindex-div(smooth_time-1,2):1:final_timeindex+div(smooth_time-1,2)
    elseif position == 3
        iterator_time_2d    = (initial_timeindex-div(smooth_time-1,2))*2-1:2:final_timeindex*2
        iterator_time_3d    = initial_timeindex-div(smooth_time-1,2):1:final_timeindex
    elseif position == 4
    end
    @info " Reading files: " file2d file3d

    ds3d                = Dataset(file3d)
    ds2d                = Dataset(file2d)
    x                   = variable(ds3d,"x")[:]    :: Array{Float32,1}
    y                   = variable(ds3d,"y")[:]    :: Array{Float32,1}
    z                   = variable(ds3d,"z")[:]    :: Array{Float32,1}
    t                   = variable(ds3d,"time")[iterator_time_3d] :: Array{Float32,1}
    P0                  = variable(ds3d,"p")[:]    :: Array{Float32,1}
    U                   = variable(ds3d,"U")[:,:,:,iterator_time_3d] :: Array{Float32,4}
    @info "Read 10%"
    V                   = variable(ds3d,"V")[:,:,:,iterator_time_3d] :: Array{Float32,4}
    @info "Read 20%"
    W                   = variable(ds3d,"W")[:,:,:,iterator_time_3d] :: Array{Float32,4}
    @info "Read 30%"
    RAD                 = variable(ds3d,"QRAD")[:,:,:,iterator_time_3d] :: Array{Float32,4}
    @info "Read 40%"
    T                   = variable(ds3d,"TABS")[:,:,:,iterator_time_3d] :: Array{Float32,4}
    @info "Read 50%"
    Tv                  = variable(ds3d,"QV")[:,:,:,iterator_time_3d] :: Array{Float32,4}
    @info "Read 60%"
    Condensate          = variable(ds3d,"QP")[:,:,:,iterator_time_3d] :: Array{Float32,4}
    Condensate          .= Condensate .+ variable(ds3d,"QN")[:,:,:,iterator_time_3d] :: Array{Float32,4}
    @info "Read 60%"
    PP                  = variable(ds3d,"PP")[:,:,:,iterator_time_3d] :: Array{Float32,4}
    @info "Read 70%"
    SHF                 = variable(ds2d,"SHF")[:,:,iterator_time_2d] :: Array{Float32,3}
    @info "Read 80%"
    LHF                 = variable(ds2d,"LHF")[:,:,iterator_time_2d] :: Array{Float32,3}
    @info "Read 90%"

    close(ds3d)
    close(ds2d)

    #if 3 nodes:
    # node 1: workers 1, 2, 5,8
    # node 2: workers 3, 6, 9
    # node 3: workers 4, 7, 10
    @info "Smoothing batch 1"
    fU   = @spawnat 3 cutborders!(filter_array(U,smooth_x,smooth_time,1),smooth_time,position)
    fV   = @spawnat 4 cutborders!(filter_array(V,smooth_x,smooth_time,1),smooth_time,position)
    fW   = @spawnat 6 cutborders!(filter_array(W,smooth_x,smooth_time,1),smooth_time,position)
    fTv  = @spawnat 7 cutborders!(filter_array(Tv,smooth_x,smooth_time,1),smooth_time,position)
    U = fetch(fU)
    V = fetch(fV)
    W = fetch(fW)
    Tv = fetch(fTv)

    fetch(@spawnat 3 GC.gc())
    fetch(@spawnat 4 GC.gc())
    fetch(@spawnat 6 GC.gc())
    fetch(@spawnat 7 GC.gc())
    @info "smoothing batch 2"
    fT   = @spawnat 3 cutborders!(filter_array(T,smooth_x,smooth_time,1),smooth_time,position)
    fPP  = @spawnat 4 cutborders!(filter_array(PP,smooth_x,smooth_time,1),smooth_time,position)
    fRAD = @spawnat 6 cutborders!(filter_array(RAD,smooth_x,smooth_time,1),smooth_time,position)
    fCondensate = @spawnat 7 cutborders!(filter_array(Condensate,smooth_x,smooth_time,1),smooth_time,position)
    SHF = cutborders!(filter_array(SHF,smooth_x,smooth_time,1),smooth_time,position)
    LHF = cutborders!(filter_array(LHF,smooth_x,smooth_time,1),smooth_time,position)
    t = cutborders!(t,smooth_time,position)

    T = fetch(fT)
    PP = fetch(fPP)
    RAD = fetch(fRAD)
    Condensate = fetch(fCondensate)
    @everywhere GC.gc()
    @everywhere GC.gc()

    ###### Processing

    @. SHF    = g/(1*Dryair.cp*Ts)*(SHF)
    @. LHF    = g/(1*Dryair.cp*Ts)*(epsilon*Dryair.cp*Ts/Liquidwater.Lv*LHF)
    @. SHF    +=  LHF     # Now it is transformed

    ThetaV       = similar(T)
    xBar_Pt      = Array{eltype(T),4}(undef,1,1,size(PP,3),size(PP,4))
    xBar_Tv      = Array{eltype(T),4}(undef,1,1,size(Tv,3),size(Tv,4))
    xBar_ThetaV  = Array{eltype(T),4}(undef,1,1,size(ThetaV,3),size(ThetaV,4))
    #******
    dx           = x[2]-x[1]
    dy           = y[2]-y[1]                            # Grid size
    kz           = length(z)                            # vertical levels
    kx           = length(x)                            # # of horizonal grid points
    ky           = length(y)                            # # of horizonal grid points
    @.  RAD      = RAD/day                              # K/s #Heating rate per second
    @.  Tv       = (1 + 1e-3*epsilon*Tv)*T              # Virtual temperature
    @.  P0       = P0*1e2
    PP          .= PP .+ reshape(P0,(1,1,kz,1))
    ThetaV      .= Tv.*(P0[1]./PP).^c1 # Virtual potential temp
    mean!(xBar_Pt,PP)
    mean!(xBar_Tv,Tv)
    mean!(xBar_ThetaV,ThetaV)
    var_Tv       =  Tv     .- xBar_Tv
    var_ThetaV   =  ThetaV .- xBar_ThetaV
    rho0         = dropdims(xBar_Pt./R./xBar_Tv,dims=(1,2))
    B            = g .* (var_ThetaV./xBar_ThetaV .- 1e-3Condensate)
    @. RAD       = RAD*(g/xBar_Tv)                    # convert unit to buoyancy
    PP           = []

    N2           = compute_N2(xBar_Tv,z)

    LHF        = []
    T          = []
    ThetaV     = []
    Tv         = []
    var_ThetaV = []
    var_Tv     = []
    @info "Computing buoyancy budget"
    ### NOTE that SHF is now the sum, saving memory
    # Buoyancy budget
    dz          = 50
    @info size(B), size(RAD), size(SHF), size(U),size(V) ,size(W), size(N2), size(dx),size(dy), size(dz), size(dt), size(x),size(y), size(z), size(t)
    Diabatic_other = fetch(@spawnat 3  get_diabatic_as_residual_buoyancy(B, RAD, SHF, U,V ,W, N2, dx,dy, dz, dt))
    @everywhere GC.gc()
    @everywhere GC.gc()
    @info "Computing ape budget for the whole troposphere"

    # APE budget
    z_up        = tropopause_height
    z_BL        = 2000

     (int_mass,
     int_KE,
     int_APE,
     int_APE_rate,
     int_APE_Ub2,
     int_APE_Vb2,
     int_APE_WN2,
     int_APE_RAD,
     int_APE_DIA,
     xBar_APE_Fs,
     residual) = fetch( @spawnat 3 getapebudget(B, U,V, W, N2, RAD, SHF, Diabatic_other, rho0, x,y, z, t, dx,dy, dz, dt, z_up))
    @everywhere GC.gc()
    @everywhere GC.gc()

    Diabatic_other .= Diabatic_other .- mean(Diabatic_other,dims=(1,2)) #They are now perturbations
    RAD            .= RAD .- mean(RAD,dims=(1,2))

    dia_ape = Diabatic_other.*B
    rad_ape = RAD .*(B .- mean(B,dims=(1,2)))

    @info "WritingData"

    jldopen(string(outputfile,partnumber,".jld"), "w") do file
        write(file,"rho0",rho0)
        write(file,"B",B)
        write(file,"N2",N2)
        write(file,"int_APE",int_APE)
        write(file,"int_KE",int_KE)
        write(file,"int_RAD",int_APE_RAD)
        write(file,"int_DIA",int_APE_DIA)
        write(file,"int_WN2",int_APE_WN2)
        write(file,"int_Ub2",int_APE_Ub2)
        write(file,"int_Vb2",int_APE_Vb2)
        write(file,"int_APE_rate",int_APE_rate)
        write(file,"APE_Fs",xBar_APE_Fs)
        write(file,"convec_heating_anomaly",Diabatic_other)
        write(file,"rad_heating_anomaly",RAD)
        write(file,"buoyancy_anomaly",B)
        write(file,"radiative_ape_production",dia_ape)
        write(file,"convective_ape_production",rad_ape)
    end
    GC.gc()
    @everywhere GC.gc()
    return nothing

end

function computebudgets_with_condensate(exp_name)

    #file2d = "/global/cscratch1/sd/aramreye/sam3d/subsetsForLin/f5e-4_2km_1000km_control_2d.nc"
    #file3d = "/global/cscratch1/sd/aramreye/sam3d/subsetsForLin/f5e-4_2km_1000km_control_3d.nc"
    outputfile = string("/global/cscratch1/sd/aramreye/for_postprocessing/ApeBudgetOutputs2/",exp_name,"_with_condensate")
    input_folder = "/global/cscratch1/sd/aramreye/for_postprocessing/largencfiles/"
    file2d = string(input_folder,exp_name,"_2d.nc")
    file3d = string(input_folder,exp_name,"_3d.nc")
    total_days = 100
    output_timestep = 7200 #seconds
    indices_in_part = 50
    indices_in_day = 86400 ÷ output_timestep
    number_of_parts = total_days*indices_in_day÷indices_in_part

    for part_number in 1:number_of_parts
        if part_number == 1
            position = 1
        elseif part_number == number_of_parts
            position = 3
        else
            position = 2
        end
        initial_timestep = indices_in_part*(part_number-1) + 1
        final_timestep = indices_in_part*(part_number)
        @info "Starting part $part_number with indices from" initial_timestep final_timestep
        flush(stdout)
        flush(stderr)
        getapeanalysis_with_condensate_loading(file2d,file3d,outputfile,initial_timestep,final_timestep,part_number,output_timestep,position,Float32)

    end

end


text:
using PyPlot, PyCall, Dates
basemap = pyimport("mpl_toolkits.basemap")

# testing thing here
# Miller projection
map = basemap.Basemap(projection="mill", lon_0=180)

# Plot coastlines, draw label meridians and parallels
#map.pcolormesh(lons,lats,sst,shading="flat",cmap="RdYlBu_r",latlon=true)
map.drawcoastlines(linewidth=0.5)
map.drawcountries(linewidth=0.25)
#map.fillcontinents(color="coral", lake_color="aqua")
#map.drawmeridians(collect(map.lonmin:60:map.lonmax+30), labels=[0,0,0,1])
#map.drawparallels(collect(-90:30:90), labels=[1,0,0,0])
#map.drawmapboundary(fill_color="aqua")

# Shade the night areas using current time in UTC
#CS = map.nightshade(now(Dates.UTC))
#title("Day/Night Map for " * Libc.strftime(time()) * " (UTC)")

ext: jl

max_stars_repo_name: rctngle/OpenWeatherDecoder

max_issues_repo_name: rctngle/OpenWeatherDecoder

max_forks_repo_name: rctngle/OpenWeatherDecoder


text:
export
    TrajdataEditMode,

    remove_offroad_vehicles!,
    remove_unaligned_vehicles!,
    interpolate_to_timestep!

type TrajdataEditMode
    roadway::Roadway
    scenes::Vector{Scene}
    time::Vector{Float64}
end
function Base.convert(::Type{TrajdataEditMode}, trajdata::Trajdata)
    roadway = trajdata.roadway
    scenes = Scene[get!(Scene(carsinframe(trajdata, frame)), trajdata, frame) for frame in 1 : nframes(trajdata)]
    time = Float64[frame.t for frame in trajdata.frames]
    TrajdataEditMode(roadway, scenes, time)
end
function Base.convert(::Type{Trajdata}, tdem::TrajdataEditMode)
    roadway = tdem.roadway
    vehicles = Dict{Int, VehicleDef}()
    states = Array{TrajdataState}(sum(s->length(s), tdem.scenes))
    frames = Array{TrajdataFrame}(length(tdem.scenes))

    states_index = 0
    for (frame_index, scene) in enumerate(tdem.scenes)
        lo = states_index + 1

        for veh in scene
            if !haskey(vehicles, veh.def.id)
                vehicles[veh.def.id] = veh.def
            else
                @assert(vehicles[veh.def.id] == veh.def)
            end
            states[states_index+=1] = TrajdataState(veh.def.id, veh.state)
        end
        frames[frame_index] = TrajdataFrame(lo, states_index, tdem.time[frame_index])
    end

    Trajdata(roadway, vehicles, states, frames)
end

function remove_offroad_vehicles!(tdem::TrajdataEditMode, threshold_lane_lateral_offset::Float64=2.5)
    for scene in tdem.scenes
        i = 0
        while i < length(scene)
            i += 1
            veh = scene[i]
            if abs(veh.state.posF.t) > threshold_lane_lateral_offset
                deleteat!(scene, i)
                i -= 1
            end
        end
    end
    tdem
end
function remove_unaligned_vehicles!(tdem::TrajdataEditMode, threshold_lane_angle::Float64=deg2rad(45))
    for scene in tdem.scenes
        i = 0
        while i < length(scene)
            i += 1
            veh = scene[i]
            if abs(veh.state.posF.ϕ) > threshold_lane_angle
                deleteat!(scene, i)
                i -= 1
            end
        end
    end
    tdem
end

function interpolate_to_timestep!(tdem::TrajdataEditMode, timestep::Float64)
    time = collect(tdem.time[1]:timestep:tdem.time[end])
    scenes = Array{Scene}(length(time))
    scenes[1] = tdem.scenes[1]

    frame_old = 1
    for i in 2 : length(time)
        t = time[i]

        # move frame_old to one past the time, if it exists
        while tdem.time[frame_old] < t
            frame_old += 1
        end

        # interpolation value ∈ [0,1]
        γ = (t - tdem.time[frame_old-1]) / (tdem.time[frame_old] - tdem.time[frame_old-1])

        scene_lo = tdem.scenes[frame_old-1]
        scene_hi = tdem.scenes[frame_old]

        # interpolate all vehicles that are in both scenes
        # if a vehicle is not in both scenes, do not inlude it
        scenes[i] = Scene(min(length(scene_lo), length(scene_hi)))
        for veh_lo in scene_lo
            id = veh_lo.def.id
            veh_index = get_index_of_first_vehicle_with_id(scene_hi, id)
            if veh_index != 0
                veh_hi = scene_hi[veh_index]
                @assert(veh_lo.def == veh_hi.def)
                veh_interp = Vehicle(lerp(veh_lo.state, veh_hi.state, γ, tdem.roadway), veh_lo.def)
                push!(scenes[i], veh_interp)
            end
        end
    end

    tdem.time = time
    tdem.scenes = scenes
    tdem
end

# smooth trajectories

# remove ghost cars


text:
using FactCheck
using BioFeatures
using Compat
include(Pkg.dir("GenomicTiles", "src","Tiles.jl"))
include(Pkg.dir("GenomicTiles", "src","Cluster.jl"))


text:
# Support Data for color maps
include("maps_data.jl")

# color scale generation
# ----------------------

"""
    colors = distinguishable_colors(n, seed=RGB{N0f8}[];
                                    dropseed=false,
                                    transform=identity,
                                    lchoices=range(0, stop=100, length=15),
                                    cchoices=range(0, stop=100, length=15),
                                    hchoices=range(0, stop=342, length=20))

Generate n maximally distinguishable colors.

This uses a greedy brute-force approach to choose n colors that are maximally
distinguishable. Given seed color(s), and a set of possible hue, chroma, and
lightness values (in LCHab space), it repeatedly chooses the next color as the
one that maximizes the minimum pairwise distance to any of the colors already
in the palette.

# Arguments

- `n`: Number of colors to generate.
- `seed`: Initial color(s) included in the palette.

# Keyword arguments

- `dropseed`: if true, the `seed` values will be dropped. This provides an easy
  mechanism to ensure that the chosen colors are distinguishable from the seed value(s).
  When true, `n` does not include the seed color(s).
- `transform`: Transform applied to colors before measuring distance. Default is `identity`;
  other choices include `deuteranopic` to simulate color-blindness.
- `lchoices`: Possible lightness values
- `cchoices`: Possible chroma values
- `hchoices`: Possible hue values

Returns a `Vector` of colors of length `n`, of the type specified in `seed`.
"""
function distinguishable_colors(n::Integer,
        seed::AbstractVector{T};
        dropseed = false,
        transform::Function = identity,
        lchoices::AbstractVector = range(0, stop=100, length=15),
        cchoices::AbstractVector = range(0, stop=100, length=15),
        hchoices::AbstractVector = range(0, stop=342, length=20)) where T<:Color

    if n <= length(seed) && !dropseed
        return seed[1:n]
    end

    # Candidate colors
    N = length(lchoices)*length(cchoices)*length(hchoices)
    candidate = Vector{Lab{Float64}}(undef, N)
    j = 0
    for h in hchoices, c in cchoices, l in lchoices
        rgb = convert(RGB, LCHab(l, c, h))
        candidate[j+=1] = convert(LCHab, rgb)
    end

    # Transformed colors
    candidate_t = Vector{Lab{Float64}}(undef, N)
    for i = 1:N
        candidate_t[i] = transform(candidate[i])
    end

    # Start with the seed colors
    n += dropseed ? length(seed) : 0
    colors = Vector{T}(undef, n)
    copyto!(colors, seed)

    # Minimum distances of the current color to each previously selected color.
    ds = fill(Inf, N)
    for i = 1:length(seed)
        ts = convert(Lab{Float64}, transform(seed[i]))::Lab{Float64}
        for k = 1:N
            ds[k] = min(ds[k], colordiff(ts, candidate_t[k]))
        end
    end

    for i in length(seed)+1:n
        j = argmax(ds)
        colors[i] = candidate[j]
        tc = candidate_t[j]
        for k = 1:N
            d = colordiff(tc, candidate_t[k])
            ds[k] = min(ds[k], d)
        end
    end

    dropseed && deleteat!(colors, 1:length(seed))

    return colors
end

distinguishable_colors(n::Integer, seed::Color; kwargs...) = distinguishable_colors(n, [seed]; kwargs...)
distinguishable_colors(n::Integer; kwargs...) = distinguishable_colors(n, Vector{RGB{N0f8}}(); kwargs...)

@deprecate distinguishable_colors(n::Integer,
                                transform::Function,
                                seed::Color,
                                ls::Vector{Float64},
                                cs::Vector{Float64},
                                hs::Vector{Float64})    distinguishable_colors(n, [seed], transform = transform, lchoices = ls, cchoices = cs, hchoices = hs)

"""
    sequential_palette(h, N::Int=100; <keyword arguments>)

Implements the color palette creation technique by [Wijffelaars, M., et al. (2008)](http://magnaview.nl/documents/MagnaView-M_Wijffelaars-Generating_color_palettes_using_intuitive_parameters.pdf).

Colormaps are formed using Bézier curves in LCHuv colorspace
with some constant hue. In addition, start and end points can be given
that are then blended to the original hue smoothly.

# Arguments

- N        - number of colors
- h        - the main hue [0,360]
- c        - the overall lightness contrast [0,1]
- s        - saturation [0,1]
- b        - brightness [0,1]
- w        - cold/warm parameter, i.e. the strength of the starting color [0,1]
- d        - depth of the ending color [0,1]
- wcolor   - starting color (warmness)
- dcolor   - ending color (depth)
- logscale - true/false for toggling logspacing
"""
function sequential_palette(h,
                            N::Int=100;
                            c=0.88,
                            s=0.6,
                            b=0.75,
                            w=0.15,
                            d=0.0,
                            wcolor=RGB(1,1,0),
                            dcolor=RGB(0,0,1),
                            logscale=false)

    function MixHue(a,h0,h1)
        M=mod(180.0+h1-h0, 360)-180.0
        mod(h0+a*M, 360)
    end
    function mix_linearly(a::C, b::C, s) where C<:Color
        base_color_type(C)((1-s)*comp1(a)+s*comp1(b), (1-s)*comp2(a)+s*comp2(b), (1-s)*comp3(a)+s*comp3(b))
    end

    pstart=convert(LCHuv, wcolor)
    p1=MSC(h)
    #p0=LCHuv(0,0,h) #original end point
    pend=convert(LCHuv, dcolor)

    #multi-hue start point
    p2l = 100 * (1. - w) + w * pstart.l
    p2h=MixHue(w,h,pstart.h)
    p2c=min(MSC(p2h,p2l), w*s*pstart.c)
    p2=LCHuv(p2l,p2c,p2h)

    #multi-hue ending point
    p0l=20.0*d
    p0h=MixHue(d,h,pend.h)
    p0c=min(MSC(p0h,p0l), d*s*pend.c)
    p0=LCHuv(p0l,p0c,p0h)

    q0=mix_linearly(p0,p1,s)
    q2=mix_linearly(p2,p1,s)
    q1=mix_linearly(q0,q2,0.5)

    pal = RGB{Float64}[]

    if logscale
        absc = exp10.(range(-2., stop=0., length=N))
    else
        absc = range(0.,stop=1.,length=N)
    end

    for t in absc
        u=1.0-t

        #Change grid to favor light colors and to be uniform along the curve
        u = (125.0-125.0*0.2^((1.0-c)*b+u*c))
        u = invBezier(u,p0.l,p2.l,q0.l,q1.l,q2.l)

        #Get color components from Bezier curves
        ll = Bezier(u, p0.l, p2.l, q0.l, q1.l, q2.l)
        cc = Bezier(u, p0.c, p2.c, q0.c, q1.c, q2.c)
        hh = Bezier(u, p0.h, p2.h, q0.h, q1.h, q2.h)

        push!(pal, convert(RGB, LCHuv(ll,cc,hh)))
    end

    pal
end

"""
    diverging_palette(h1, h2, N::Int=100; <keyword arguments>)

Create diverging palettes by combining 2 sequential palettes

# Arguments

- N        - number of colors
- h1       - the main hue of the left side [0,360]
- h2       - the main hue of the right side [0,360]
- c        - the overall lightness contrast [0,1]
- s        - saturation [0,1]
- b        - brightness [0,1]
- w        - cold/warm parameter, i.e. the strength of the starting color [0,1]
- d1       - depth of the ending color in the left side [0,1]
- d2       - depth of the ending color in the right side [0,1]
- wcolor   - starting color (warmness)
- dcolor1  - ending color of the left side (depth)
- dcolor2  - ending color of the right side (depth)
- logscale - true/false for toggling logspacing
"""
function diverging_palette(h1,
                           h2,
                           N::Int=100;
                           mid=0.5,
                           c=0.88,
                           s=0.6,
                           b=0.75,
                           w=0.15,
                           d1=0.0,
                           d2=0.0,
                           wcolor=RGB(1,1,0),
                           dcolor1=RGB(1,0,0),
                           dcolor2=RGB(0,0,1),
                           logscale=false)

    if isodd(N)
        n=N-1
    else
        n=N
    end
    N1 = max(ceil(Int, mid*n), 1)
    N2 = Int(max(n-N1, 1))

    pal1 = sequential_palette(h1, N1+1, w=w, d=d1, c=c, s=s, b=b, wcolor=wcolor, dcolor=dcolor1, logscale=logscale)
    pal1 = reverse(pal1; dims=1)

    pal2 = sequential_palette(h2, N2+1, w=w, d=d2, c=c, s=s, b=b, wcolor=wcolor, dcolor=dcolor2, logscale=logscale)

    if isodd(N)
        midcol = weighted_color_mean(0.5, pal1[end], pal2[1])
        return [pal1[1:end-1]; midcol; pal2[2:end]]
    else
        return [pal1[1:end-1]; pal2[2:end]]
    end
end

# Colormap
# ----------------------
# Main function to handle different predefined colormaps
#
"""
    colormap(cname, N=100; mid=0.5, logscale=false, kvs...])

Returns a predefined sequential or diverging colormap computed using
the algorithm by Wijffelaars, M., et al. (2008).

Sequential colormaps `cname` choices are:

- `Blues`
- `Greens`
- `Grays`
- `Oranges`
- `Purples`,
- `Reds`

Diverging colormap choices are `RdBu`.

Optionally, you can specify the number of colors `N` (default 100).

Extra control is provided by keyword arguments.
- `mid` sets the position of the midpoint for diverging colormaps.
- `logscale=true` uses logarithmically-spaced steps in the colormap.
You can also use keyword argument names that match the argument names in
[`sequential_palette`](@ref) or [`diverging_palette`](@ref).
"""
function colormap(cname::AbstractString, N::Int=100; mid=0.5, logscale=false, kvs...)

    cname = lowercase(cname)
    if haskey(colormaps_sequential, cname)
        allowedkeys = [:h, :w, :d, :c, :s, :b, :wcolor, :dcolor]
        p = copy(colormaps_sequential[cname][1:8])

        for (k,v) in kvs
            ind = findfirst(e->e==k, allowedkeys)
            ind == nothing && throw(ArgumentError("Unknown keyword argument $k"))
            if ind > 0
                p[ind] = v
            end
        end

        return sequential_palette(p[1], N, w=p[2], d=p[3], c=p[4], s=p[5], b=p[6], wcolor=p[7],
                                  dcolor=p[8], logscale=logscale)

    elseif haskey(colormaps_diverging, cname)
        allowedkeys = [:h1, :h2, :w, :d1, :d2, :c, :s, :b, :wcolor, :dcolor1, :dcolor2]
        p = copy(colormaps_diverging[cname][1:11])

        for (k,v) in kvs
            ind = findfirst(e->e==k, allowedkeys)
            ind == nothing && throw(ArgumentError("Unknown keyword argument $k"))
            if ind > 0
                p[ind] = v
            end
        end

        return diverging_palette(p[1], p[2], N, w=p[3], d1=p[4], d2=p[5], c=p[6], s=p[7],
                                 b=p[8], wcolor=p[9], dcolor1=p[10], dcolor2=p[11], mid=mid,
                                 logscale=logscale)

    else
        throw(ArgumentError(string("Unknown colormap: ", cname)))
    end
end


text:
using POMDPs
using POMCPOW
using ProfileView
using POMDPModels

#=
using Gallium
breakpoint(Pkg.dir("POMCPOW", "src", "solver2.jl"), 100)
=#

solver = POMCPOWSolver(tree_queries=250_000,
                     eps=0.01,
                     criterion=MaxUCB(10.0),
                     enable_action_pw=true,
                     check_repeat_obs=false,
                     alpha_observation=1/8,
                     rng=MersenneTwister(2))

problem = LightDark1D()
policy = POMCPOWPlanner(solver, problem)
ib = initial_state_distribution(problem)

a = action(policy, ib)

@time a = action(policy, ib)

Profile.clear()
@profile a = action(policy, ib)
ProfileView.view()

#=
P = typeof(problem)
S = state_type(P)
A = action_type(P)
O = obs_type(P)
tree = POMCPOWTree{POWNodeBelief{S,A,O,P},A,O,typeof(ib)}(ib, 500_000)

@code_warntype POMCPOW.simulate(policy, POMCPOW.POWTreeObsNode(tree, 1), rand(Base.GLOBAL_RNG, ib), 10)
=#

#=
solver = POMCPOWSolver(tree_queries=30,
                     eps=0.01,
                     criterion=MaxUCB(10.0),
                     enable_action_pw=true,
                     check_repeat_obs=false,
                     alpha_observation=1/8,
                     rng=MersenneTwister(2))

policy = POMCPOWPlanner(solver, problem)
a = action(policy, ib)
blink(policy)
=#


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

using Markdown
using InteractiveUtils

# ╔═╡ 85412e62-9770-4a85-92cb-0b2366ba5a81
md"""
# 3월10일 강의영상

"""

# ╔═╡ 29b25571-c231-485a-9957-3dee5a9065e8
html"""

"""

# ╔═╡ 03bfb1bb-6de0-4836-92e6-5abafaad1784

# ╔═╡ bc068777-a366-4552-a8d7-d2fd4429f01c
md"# 3월10일 강의노트"

# ╔═╡ 7fd5c088-5ac0-432c-8aa5-251823a15dad
md"### 변수이름"

# ╔═╡ 8688fb1a-e5d1-498f-8fd5-2944191a832e
md"""
`-` 일반적인 변수이름
"""

# ╔═╡ 93015177-90dd-4a75-85d7-d7c4d22545dd
x=1.0

# ╔═╡ 7858119d-2300-4604-8a4a-e7f3d7a82fd2
md"""
`-` 좀 신기해보이는 변수이름도 가능함.
"""

# ╔═╡ 5f9cc659-444a-4aca-835d-4ebd3bc2db8a
δ=0.1 # \delta + tab

# ╔═╡ 037ea21a-bc67-4ede-9628-50d9ad29a6e7
π # 미리저장되어있음

# ╔═╡ 1674553a-8e42-4473-b574-6c814a8091fb
ℯ # 미리저장되어 있음 \euler + tab

# ╔═╡ fed45320-b7f4-478b-8a2a-bfb417f5e34d
x̄ = (1+2+3+3)/4 # x\bar + tab

# ╔═╡ 1f96e84b-4c3e-4fdc-8441-8fb4d3a49b3d
x̂ = 5.0 # x\hat + tab

# ╔═╡ 6100965a-ec16-4274-9411-64f9dedc67a0
β̂ = 4.0 # \beta + tab -> \hat + tab

# ╔═╡ 7558f29d-b09d-4e89-aa3e-0493cf735b49
xₜ = 2.0 # x\_t + tab

# ╔═╡ be0ed6cb-fbb1-427f-bd9b-ab6b2f96bf70
β̂₀ = 1.5 # \beta + tab -> \hat + tab -> \_0 + tab

# ╔═╡ 031ccf5a-fac5-4b03-97e1-abf1f88824e8
β̂ᵐˡᵉ = 2.2

# ╔═╡ a938bbfc-566d-4503-9543-24470a1568a5
♂ = "male"

# ╔═╡ 4e96a402-3eac-4a03-a3fb-050b67e5f165
최규빈 = "GuebinChoi"

# ╔═╡ 69ab7085-ffd6-41ea-9ac3-861cae5aa6cc
md"### 자료형"

# ╔═╡ 7a9ec94b-4bf5-45e9-ace2-9d57789af169
md"##### int, float"

# ╔═╡ aa3c4cfd-e260-47fc-aa96-103995d99216
123 # int

# ╔═╡ e5f215a7-cff9-4274-931a-36ff0c6bb63a
1.23 # float

# ╔═╡ c5a27cde-cccf-45bd-b278-6712744714a3
md"""
`-` 타입을 보려면?
"""

# ╔═╡ 3f7b6a98-8b81-4dda-85d8-01ed3bf498d6
typeof(123)

# ╔═╡ 7be5dda7-ec6a-4bac-b436-5efe82630bfe
typeof(1.23)

# ╔═╡ 1f09e243-6d8e-450c-a396-5781656ffaf2
md"""
##### Inf, NaN
"""

# ╔═╡ ecd47dbd-784a-448b-8bd9-bdf71b7c6192
md"""
`-` Inf
"""

# ╔═╡ 0de4283c-3416-4ed3-90b4-59785bac2add
Inf

# ╔═╡ 52ecb331-9838-4039-bc2d-d8baff6bd394
typeof(Inf)

# ╔═╡ a594053f-73ed-43cb-88ac-ee2306d1a32d
1/Inf

# ╔═╡ 3dc2979b-6f03-4337-af32-6b039475c198
1/0

# ╔═╡ 92ba7de9-87a5-4d2c-a6b9-da3fe695463c
Inf + 100

# ╔═╡ 1a34c0dd-3f49-45a5-968d-f8ac0ba70f96
Inf + Inf

# ╔═╡ 7d231a82-0126-4566-8b0b-308c3f45e065
500 - Inf

# ╔═╡ 10461e54-a79a-4c26-bfd0-f657b135de80
md"""
`-` NaN
"""

# ╔═╡ 213b178c-6c93-497c-ac33-5ffce4975a97
NaN

# ╔═╡ 03b6c227-b712-432e-9259-b5d349bed203
typeof(NaN)

# ╔═╡ 332839be-a047-483f-ae68-a453adee9998
0/0

# ╔═╡ e9b47540-4859-4d25-baff-b054d408afbe
typeof(0/0)

# ╔═╡ c0f74c89-01a5-4b2b-b8fe-ed0c50edfb38
Inf - Inf

# ╔═╡ a6b58372-9f41-45a2-87f7-ffe8da31cdda
Inf / Inf

# ╔═╡ fd9ebe5a-6c7f-4c97-8474-69cf100e4f0b
0 * Inf

# ╔═╡ aae82467-bf7d-4065-954a-0368aca73d8e
NaN + 2

# ╔═╡ d34f35b8-eaac-4097-82b0-c1018590bf0d
NaN - NaN

# ╔═╡ 34d1f92f-8861-4c1b-b1ec-b41a95ccb32a
NaN + Inf

# ╔═╡ e6350b7c-8943-4d10-be01-888112af125c
md"""
`-` Inf, Nan의 비교
- tip: 잘 모르겠으면 false
- 저는 이걸 이용해서 코딩하진 않아요.
- 그래도 성질은 한번씩 체크해보는 것은 좋음.
- 그리고 가끔 Inf를 쓰면 코딩이 편리할때가 있어요.
"""

# ╔═╡ 89daf982-a9c1-4583-a1bc-2e6366c6f125
Inf > 2

# ╔═╡ 3699cac6-046d-467d-b756-ce38dd42d953
Inf < 2

# ╔═╡ ef6c82fa-b164-4321-b539-1b3577ad431a
-Inf > 2

# ╔═╡ 364c7521-f34e-44fe-9aff-8fcece4f8f10
-Inf < 2

# ╔═╡ ff3a8440-164e-487e-aad7-d17f6e482e8c
-Inf < Inf

# ╔═╡ cb70b33d-f8ee-4e89-a5ee-90bc529460ea
Inf == Inf

# ╔═╡ 20df59d8-614c-4f10-914a-2e8b5c6ade45
Inf != Inf

# ╔═╡ 1f56d6f2-a59e-4a40-99e0-0abe9365bef1
NaN > 2

# ╔═╡ c27f0b69-796d-432c-8923-20c8868179ca
NaN == NaN

# ╔═╡ 07918808-d082-432a-8115-868ab95cba3f
NaN > NaN

# ╔═╡ 6947b7bb-bbd7-4b50-84f3-b3c760bdca93
NaN < NaN

# ╔═╡ 831328c6-fadf-4cff-9ab8-cc9d1839f6a8
NaN != NaN

# ╔═╡ 99bb3338-cc43-43aa-81ba-27d27afd61f6
NaN > Inf

# ╔═╡ 51a605f8-8e20-41d9-97c3-65b0c8aba622
NaN == Inf

# ╔═╡ d7f7edf5-e3c1-4ab2-aeee-d01c1fff1e11
NaN < Inf

# ╔═╡ cb50e1d2-f5fb-40f7-a97e-dbacfe60913d
NaN + Inf

# ╔═╡ 6215cb02-5b52-445b-bfa8-eeac87470537
md"""
##### 문자열
"""

# ╔═╡ 2ef84c02-f271-4bae-a50e-cfd6738a99ec
str1 = "asdf"

# ╔═╡ 1ab70342-5bd3-4c06-942c-599b06c80e98
str1[2] # 두번째 원소

# ╔═╡ 2922318a-9db5-4fe9-9f3b-d70e8b53eae6
str1[1:2]

# ╔═╡ 83ee9ee0-0eab-41c6-846d-9352338aca8b
length(str1)

# ╔═╡ c2909e6f-e0d2-4be6-9a36-6b53b622fa47
str1 * "efg" # 문자열 잇기 방법1

# ╔═╡ da88c5ab-f756-4c48-b32e-bc69482cbf69
string(str1 , "efg") # 문자열 잇기 방법2

# ╔═╡ 8b5c2f32-62f8-4973-ae0c-83a084868c58
"1더하기1은 $(1+1)" # 문자열 끼워넣기

# ╔═╡ c8497ded-75a5-481c-8a65-1b2bf3c190a0
begin
	word1 = "Julia"
	word2 = "data"
	word3 = "science"
end

# ╔═╡ 8ec75baf-85d1-43f6-8c6c-ef3e19536529
"$word1 is great for $word2 $word3"

# ╔═╡ 43f529e2-dd3f-472d-852b-bba8853ad9d8
md"""
##### 튜플
"""

# ╔═╡ ce71647d-fb14-40c3-b3e3-5d094f7e4187
(1,2,3)

# ╔═╡ 8c81c79a-4ce5-40be-9049-493067fccfec
1,2,3 # 괄호는 생략가능

# ╔═╡ 2ceac045-aac4-44f7-b969-d727595f30cb
a = (11,22,33) # a라는 변수에 튜플 (1,2,3)을 넣음

# ╔═╡ 2affadf8-d96c-467a-b725-d34888983acb
a[1]

# ╔═╡ 2d288956-7028-4356-a987-ee8db76ae472
a[1]=100

# ╔═╡ c5ad6523-07c2-444c-9cd3-8e1ba61cfef8
a₁, a₂, a₃ = 11,22,33
# (a₁, a₂, a₃) = (11,22,33)
# ==> a₁=11, a₂=22, a₃=33

# ╔═╡ e9df4f75-fec8-4fa0-8338-bd67c844bf4b
a₁

# ╔═╡ 0dc6ff92-b505-4ebc-97d9-a7761262f626
a₂

# ╔═╡ 9e83f7e9-9571-4195-b42b-62e12ea4d180
a₃

# ╔═╡ 5c347108-86f6-469c-9db3-a9e1a1e3c478
md"""
##### Array
"""

# ╔═╡ 556a13bc-2d9f-48ef-abe9-144ec33264bd
b = [11,22,33]

# ╔═╡ d738dd15-b1bb-4e55-b788-c9212a488ed1
b[1]

# ╔═╡ f0e05ff9-820d-4b95-ad0d-770c72589a3b
b[1] = 100

# ╔═╡ 3406a77f-9ebc-4bd2-9082-2dafda0a1f08
b

# ╔═╡ 72c6cdb7-608f-44a6-8f2e-6dfdf4af3bcc
md"""
`-` 브로드캐스팅
"""

# ╔═╡ 306a8f60-8cad-49ce-934f-60294c85069b
b+1

# ╔═╡ 28933892-c486-48f0-b420-e4ff03d26ede
broadcast(+,b,1)

# ╔═╡ 87c766e5-ba46-43f2-b840-02b8566056e1
b.+1

# ╔═╡ 9612c517-e2e4-4321-9045-e9cc70c8d650
md"""
`-` 아래는 다르다
"""

# ╔═╡ 5be10235-e7cc-4b9a-9638-0ce11ddb8e95
b1 = [1,2,3]

# ╔═╡ f76342cd-6bf4-41ad-b0ba-e3ea67c966b4
typeof(b1)

# ╔═╡ 1b8eb691-dd34-4b78-bb2f-29b563ed8892
b2 = [1 2 3]

# ╔═╡ c3b90a24-348a-40fe-bd9c-c2da508a11ea
typeof(b2)

# ╔═╡ 6e57bf96-a86c-41f3-96d6-82923ac359f3
md"""
`-` b1,b2를 브로드캐스팅 하면?
"""

# ╔═╡ e3f215a9-b173-4c9c-85cd-f0872cf9b0f5
b1.*b2

# ╔═╡ cc10ac69-667a-4a41-8439-0aea0406e140
md"""
`-` 브로드캐스팅의 또 다른 응용
"""

# ╔═╡ 546c0949-cbb7-4376-980b-cb5bae22b9ab
[1,2] == [1,2]

# ╔═╡ d588d8fa-2978-435b-8c1b-fad79c7d5810
[1,2] == [2,2]

# ╔═╡ cab2a579-7116-428a-bd45-4dff1eab3604
[1,2] == [1,1]

# ╔═╡ c7d4b4f3-74f2-4fa5-aaae-62dbfb3a4b3a
[1,2] == [2,1]

# ╔═╡ 493d28a8-4295-4eb5-a159-5f2b1fd7942a
[1,2] .== [1,2]

# ╔═╡ 98ec2369-df57-4668-a6ec-c738c52c458b
[1,2] .== [2,2]

# ╔═╡ 7155f100-1d83-4c04-a30c-d7ba45ee74e8
[1,2] .== [1,1]

# ╔═╡ 280937bb-9fa6-423c-874b-3cf814b6bd10
[1,2] .== [2,1]

# ╔═╡ 0de890bc-be78-487d-9001-190cc8664587
[1,2,3,4] .== [1,2,33,4]

# ╔═╡ e8e7d4c7-70f0-45ab-9469-26cd6f70a7f6
md"""
`-` 난수로 벡터선언
"""

# ╔═╡ d966ce95-faf6-484a-87a8-6ffd273f31ff
randn(5) # 정규분포

# ╔═╡ 7b951345-bfa5-464d-a71c-ff067d427010
rand(5) # rand의 첫번째 기능: 유니폼분포

# ╔═╡ 2399aaa0-6062-4517-ac38-be9d4ed8a25b
rand([1,2,3,4],2) # rand의 두번째 기능: sample

# ╔═╡ 9c2f38cd-3780-43de-b5d1-5e757613cfbf
md"""
`-` 리스트 컴프리헨션
"""

# ╔═╡ 67f92c8d-7167-4d28-8ab0-108b84e70c22
[√1, √2, √3, √4]

# ╔═╡ cf4cee44-2b0b-4219-837d-d3fc2309add0
[√i for i in [1,2,3,4]]

# ╔═╡ 94369be3-b932-4ae4-8c1e-fd306876a3b2

# ╔═╡ 4bca7dd3-5e5d-4eda-a98d-e9f2b0a26dc5
#과제 #

# ╔═╡ ec9be90d-4770-4e25-9703-551125bdaaca
[i for i in [2^0, 2^1, 2^2, 2^3]]

# ╔═╡ fd6ade5d-9c6f-47bd-85e0-13f83c00ac5b
md"""
- 이렇게 보면 편합니다. $\{\sqrt{1},\sqrt{2},\sqrt{3},\sqrt{4} \}= \{\sqrt{i}: i=1,2,3,4 \}$
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.2"
manifest_format = "2.0"

[deps]
"""

# ╔═╡ Cell order:
# ╠═85412e62-9770-4a85-92cb-0b2366ba5a81
# ╠═29b25571-c231-485a-9957-3dee5a9065e8
# ╠═03bfb1bb-6de0-4836-92e6-5abafaad1784
# ╠═bc068777-a366-4552-a8d7-d2fd4429f01c
# ╠═7fd5c088-5ac0-432c-8aa5-251823a15dad
# ╟─8688fb1a-e5d1-498f-8fd5-2944191a832e
# ╠═93015177-90dd-4a75-85d7-d7c4d22545dd
# ╟─7858119d-2300-4604-8a4a-e7f3d7a82fd2
# ╠═5f9cc659-444a-4aca-835d-4ebd3bc2db8a
# ╠═037ea21a-bc67-4ede-9628-50d9ad29a6e7
# ╠═1674553a-8e42-4473-b574-6c814a8091fb
# ╠═fed45320-b7f4-478b-8a2a-bfb417f5e34d
# ╠═1f96e84b-4c3e-4fdc-8441-8fb4d3a49b3d
# ╠═6100965a-ec16-4274-9411-64f9dedc67a0
# ╠═7558f29d-b09d-4e89-aa3e-0493cf735b49
# ╠═be0ed6cb-fbb1-427f-bd9b-ab6b2f96bf70
# ╠═031ccf5a-fac5-4b03-97e1-abf1f88824e8
# ╠═a938bbfc-566d-4503-9543-24470a1568a5
# ╠═4e96a402-3eac-4a03-a3fb-050b67e5f165
# ╟─69ab7085-ffd6-41ea-9ac3-861cae5aa6cc
# ╟─7a9ec94b-4bf5-45e9-ace2-9d57789af169
# ╠═aa3c4cfd-e260-47fc-aa96-103995d99216
# ╠═e5f215a7-cff9-4274-931a-36ff0c6bb63a
# ╟─c5a27cde-cccf-45bd-b278-6712744714a3
# ╠═3f7b6a98-8b81-4dda-85d8-01ed3bf498d6
# ╠═7be5dda7-ec6a-4bac-b436-5efe82630bfe
# ╟─1f09e243-6d8e-450c-a396-5781656ffaf2
# ╟─ecd47dbd-784a-448b-8bd9-bdf71b7c6192
# ╠═0de4283c-3416-4ed3-90b4-59785bac2add
# ╠═52ecb331-9838-4039-bc2d-d8baff6bd394
# ╠═a594053f-73ed-43cb-88ac-ee2306d1a32d
# ╠═3dc2979b-6f03-4337-af32-6b039475c198
# ╠═92ba7de9-87a5-4d2c-a6b9-da3fe695463c
# ╠═1a34c0dd-3f49-45a5-968d-f8ac0ba70f96
# ╠═7d231a82-0126-4566-8b0b-308c3f45e065
# ╟─10461e54-a79a-4c26-bfd0-f657b135de80
# ╠═213b178c-6c93-497c-ac33-5ffce4975a97
# ╠═03b6c227-b712-432e-9259-b5d349bed203
# ╠═332839be-a047-483f-ae68-a453adee9998
# ╠═e9b47540-4859-4d25-baff-b054d408afbe
# ╠═c0f74c89-01a5-4b2b-b8fe-ed0c50edfb38
# ╠═a6b58372-9f41-45a2-87f7-ffe8da31cdda
# ╠═fd9ebe5a-6c7f-4c97-8474-69cf100e4f0b
# ╠═aae82467-bf7d-4065-954a-0368aca73d8e
# ╠═d34f35b8-eaac-4097-82b0-c1018590bf0d
# ╠═34d1f92f-8861-4c1b-b1ec-b41a95ccb32a
# ╟─e6350b7c-8943-4d10-be01-888112af125c
# ╠═89daf982-a9c1-4583-a1bc-2e6366c6f125
# ╠═3699cac6-046d-467d-b756-ce38dd42d953
# ╠═ef6c82fa-b164-4321-b539-1b3577ad431a
# ╠═364c7521-f34e-44fe-9aff-8fcece4f8f10
# ╠═ff3a8440-164e-487e-aad7-d17f6e482e8c
# ╠═cb70b33d-f8ee-4e89-a5ee-90bc529460ea
# ╠═20df59d8-614c-4f10-914a-2e8b5c6ade45
# ╠═1f56d6f2-a59e-4a40-99e0-0abe9365bef1
# ╠═c27f0b69-796d-432c-8923-20c8868179ca
# ╠═07918808-d082-432a-8115-868ab95cba3f
# ╠═6947b7bb-bbd7-4b50-84f3-b3c760bdca93
# ╠═831328c6-fadf-4cff-9ab8-cc9d1839f6a8
# ╠═99bb3338-cc43-43aa-81ba-27d27afd61f6
# ╠═51a605f8-8e20-41d9-97c3-65b0c8aba622
# ╠═d7f7edf5-e3c1-4ab2-aeee-d01c1fff1e11
# ╠═cb50e1d2-f5fb-40f7-a97e-dbacfe60913d
# ╠═6215cb02-5b52-445b-bfa8-eeac87470537
# ╠═2ef84c02-f271-4bae-a50e-cfd6738a99ec
# ╠═1ab70342-5bd3-4c06-942c-599b06c80e98
# ╠═2922318a-9db5-4fe9-9f3b-d70e8b53eae6
# ╠═83ee9ee0-0eab-41c6-846d-9352338aca8b
# ╠═c2909e6f-e0d2-4be6-9a36-6b53b622fa47
# ╠═da88c5ab-f756-4c48-b32e-bc69482cbf69
# ╠═8b5c2f32-62f8-4973-ae0c-83a084868c58
# ╠═c8497ded-75a5-481c-8a65-1b2bf3c190a0
# ╠═8ec75baf-85d1-43f6-8c6c-ef3e19536529
# ╟─43f529e2-dd3f-472d-852b-bba8853ad9d8
# ╠═ce71647d-fb14-40c3-b3e3-5d094f7e4187
# ╠═8c81c79a-4ce5-40be-9049-493067fccfec
# ╠═2ceac045-aac4-44f7-b969-d727595f30cb
# ╠═2affadf8-d96c-467a-b725-d34888983acb
# ╠═2d288956-7028-4356-a987-ee8db76ae472
# ╠═c5ad6523-07c2-444c-9cd3-8e1ba61cfef8
# ╠═e9df4f75-fec8-4fa0-8338-bd67c844bf4b
# ╠═0dc6ff92-b505-4ebc-97d9-a7761262f626
# ╠═9e83f7e9-9571-4195-b42b-62e12ea4d180
# ╟─5c347108-86f6-469c-9db3-a9e1a1e3c478
# ╠═556a13bc-2d9f-48ef-abe9-144ec33264bd
# ╠═d738dd15-b1bb-4e55-b788-c9212a488ed1
# ╠═f0e05ff9-820d-4b95-ad0d-770c72589a3b
# ╠═3406a77f-9ebc-4bd2-9082-2dafda0a1f08
# ╟─72c6cdb7-608f-44a6-8f2e-6dfdf4af3bcc
# ╠═306a8f60-8cad-49ce-934f-60294c85069b
# ╠═28933892-c486-48f0-b420-e4ff03d26ede
# ╠═87c766e5-ba46-43f2-b840-02b8566056e1
# ╟─9612c517-e2e4-4321-9045-e9cc70c8d650
# ╠═5be10235-e7cc-4b9a-9638-0ce11ddb8e95
# ╠═f76342cd-6bf4-41ad-b0ba-e3ea67c966b4
# ╠═1b8eb691-dd34-4b78-bb2f-29b563ed8892
# ╠═c3b90a24-348a-40fe-bd9c-c2da508a11ea
# ╟─6e57bf96-a86c-41f3-96d6-82923ac359f3
# ╠═e3f215a9-b173-4c9c-85cd-f0872cf9b0f5
# ╟─cc10ac69-667a-4a41-8439-0aea0406e140
# ╠═546c0949-cbb7-4376-980b-cb5bae22b9ab
# ╠═d588d8fa-2978-435b-8c1b-fad79c7d5810
# ╠═cab2a579-7116-428a-bd45-4dff1eab3604
# ╠═c7d4b4f3-74f2-4fa5-aaae-62dbfb3a4b3a
# ╠═493d28a8-4295-4eb5-a159-5f2b1fd7942a
# ╠═98ec2369-df57-4668-a6ec-c738c52c458b
# ╠═7155f100-1d83-4c04-a30c-d7ba45ee74e8
# ╠═280937bb-9fa6-423c-874b-3cf814b6bd10
# ╠═0de890bc-be78-487d-9001-190cc8664587
# ╟─e8e7d4c7-70f0-45ab-9469-26cd6f70a7f6
# ╠═d966ce95-faf6-484a-87a8-6ffd273f31ff
# ╠═7b951345-bfa5-464d-a71c-ff067d427010
# ╠═2399aaa0-6062-4517-ac38-be9d4ed8a25b
# ╟─9c2f38cd-3780-43de-b5d1-5e757613cfbf
# ╠═67f92c8d-7167-4d28-8ab0-108b84e70c22
# ╠═cf4cee44-2b0b-4219-837d-d3fc2309add0
# ╠═94369be3-b932-4ae4-8c1e-fd306876a3b2
# ╠═4bca7dd3-5e5d-4eda-a98d-e9f2b0a26dc5
# ╠═ec9be90d-4770-4e25-9703-551125bdaaca
# ╠═fd6ade5d-9c6f-47bd-85e0-13f83c00ac5b
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)
include("fwi_util.jl")
include("fwi_util_op.jl")
np = pyimport("numpy")

# argsparse.jl
# ENV["CUDA_VISIBLE_DEVICES"] = 1
# ENV["PARAMDIR"] = "Src/params/"
# config = tf.ConfigProto(device_count = Dict("GPU"=>0))

nz = 400
nx = 400
dz = 20
dx = 20
nSteps = 2001
dt = 0.0025
f0 = 4.5
filter_para = [0, 0.1, 100.0, 200.0]
nPml = 32
isAc = false
nPad = 0
# x_src = collect(5:10:nx-2nPml-5)
# z_src = 2ones(Int64, size(x_src))
# x_rec = collect(5:100-nPml)
# z_rec = 2ones(Int64, size(x_rec))

x_src = [200-nPml]
z_src = [200-nPml]

z = (5:10:nz-2nPml-5)|>collect
x = (5:10:nx-2nPml-5)|>collect
x_rec, z_rec = np.meshgrid(x, z)
x_rec = x_rec[:]
z_rec = z_rec[:]

# x_src = 5
# z_src = [300-nPml]
# z_rec = collect(5:1:nz-2nPml-5)
# x_rec = (nx-2nPml-100) .* ones(Int64, size(z_rec))

para_fname = "./para_file.json"
survey_fname = "./survey_file.json"
data_dir_name = "./Data"
scratch_dir_name="./Scratch"
# paraGen(nz, nx, dz, dx, nSteps, dt, f0, nPml, nPad, filter_para, isAc, para_fname, survey_fname, data_dir_name, scratch_dir_name=scratch_dir_name)
# surveyGen(z_src, x_src, z_rec, x_rec, survey_fname)
paraGen(nz, nx, dz, dx, nSteps, dt, f0, nPml,
        nPad, para_fname, survey_fname, data_dir_name,
        scratch_dir_name=scratch_dir_name)
surveyGen(z_src, x_src, z_rec, x_rec, survey_fname)

cp = 3000ones(nz, nx)
# cp = (1. .+ 0.1*rand(nz, nx)) .* 3000.
cs = 3000.0/sqrt(3.0) .* ones(nz,nx)
# cs = zeros(nz,nx)
den = 2000.0 .* ones(nz, nx)
function vel2moduli(cp,cs,den)
    lambda = (cp.^2 - 2.0 .* cs.^2) .* den ./ 1e6
    mu = cs.^2 .* den ./ 1e6
    return lambda, mu
end
lambda, mu = vel2moduli(cp,cs,den)

tf_lambda = constant(lambda)
tf_mu = constant(mu)
tf_den = constant(den)

# # src = Matrix{Float64}(undef, 1, 2001)
# src[1,:] = Float64.(reinterpret(Float32, read("./Src/params/ricker_10Hz.bin")))
src = sourceGene(f0, nSteps, dt)
tf_stf = constant(repeat(src, outer=length(z_src)))
tf_para_fname = tf.strings.join([para_fname])
tf_gpu_id0 = constant(0, dtype=Int32)
tf_gpu_id1 = constant(1, dtype=Int32)
tf_shot_ids0 = constant(collect(Int32, 0:length(x_src)-1), dtype=Int32)
tf_shot_ids1 = constant(collect(Int32, 13:25), dtype=Int32)

res1 = fwi_obs_op(tf_lambda, tf_mu, tf_den, tf_stf, tf_gpu_id0, tf_shot_ids0, tf_para_fname)
# res2 = fwi_obs_op(tf_cp2, tf_cs2, tf_den2, tf_stf, tf_gpu_id1, tf_shot_ids0, tf_para_fname)

sess=Session();init(sess);
@time run(sess, res1)
# error("")

# gradient check -- v
function scalar_function(m)
    # return fwi_op(m, tf_mu, tf_den, tf_stf, tf_gpu_id0, tf_shot_ids0, tf_para_fname)
    return fwi_op(tf_lambda, m, tf_den, tf_stf, tf_gpu_id0, tf_shot_ids0, tf_para_fname)
    # return fwi_op(tf_lambda, tf_mu, m, tf_stf, tf_gpu_id0, tf_shot_ids0, tf_para_fname)
    # return fwi_op(tf_lambda, tf_mu, tf_den, m, tf_gpu_id0, tf_shot_ids0, tf_para_fname)
end

# lambda2,_ = vel2moduli(3200.0, 3200.0/sqrt(3.0), den)
# m_ = constant(lambda2)

_,mu2 = vel2moduli(3200.0, 3200.0/sqrt(3.0), den)
m_ = constant(mu2)

# m_ = constant(2200ones(nz,nx))

# # src2 = circshift(src, (0,30))
# src2 = sourceGene(f0,nSteps,dt) .*1.5
# m_ = constant(repeat(src2, outer=length(z_src)))

# for forward_backward wavefield comparison
# yy = scalar_function(m_)
# gradm = gradients(yy, m_)
# sess=Session();init(sess);
# @time G = run(sess, gradm)
# imshow(G);colorbar();
# A=read("SnapGPU.bin");A=reshape(reinterpret(Float32,A),(200,200));
# B=read("SnapGPU_back.bin");B=reshape(reinterpret(Float32,B),(200,200));
# imshow(A[33:end-32,33:end-32]-B[33:end-32,33:end-32]);colorbar();
# error("")

v0 = zeros(nz, nx)
# PLEASE!!!!!!!!!!!!!! Don't perturb in the CPML region!!!!!!!!!!!!!!!!!!!!!!!
v0[nPml+5:nz-nPml-5, nPml+5:nx-nPml-5] .= 1.0
# # perturb moduli
v_ = constant(Float64.((1. .+ 0.1*rand(nz, nx)) .* 1e3 .* v0))
# # perturb density
# v_ = constant(Float64.((1. .+ 0.1*rand(nz, nx)) .* 500 .* v0))
# perturb sft

# s0 = zeros(1, nSteps)
# s0[1, 20:end-20] .= 1.0
# src_perturb = rand(1,nSteps) .* s0 *0.1
# v_ = src_perturb

y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_ * v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session()
init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
# error("")
sval_ = [x[1] for x in sval_]
wval_ = [x[1] for x in wval_]
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="FWI gradient")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")

savefig("Convergence_test_mu.pdf", bbox_inches="tight",pad_inches = 0, dpi = 300);


text:
"""
Functions for applying empirical Bayes to a set of test statistics.
""";

using Discretizers
using Distributions
using Interpolations
using DataFrames
using GLM

"""
    discretize_test_statistics(test_statistics::AbstractArray{<:Real}, n)

Bin the values in test_statistics into n uniform width bins.
Returns the midpoints of the bins, the counts in each bin and the bin width.
"""
function discretize_test_statistics(test_statistics::AbstractArray{<:Real}, n)
    try
        n = Int(n) # Bad form, but allows values like 10.0 to be silently used as 10
    catch e
        error("n must be an Int or convertible to an Int")
    end

    # Calculate midpoints
    min_ts, max_ts = extrema(test_statistics)
    bin_width = (max_ts - min_ts) / n
    midpoints = [min_ts + i*bin_width - bin_width/2 for i in 1:n]

    # Get counts for each bin
    bin_edges = binedges(DiscretizeUniformWidth(n), test_statistics)
    lindisc = LinearDiscretizer(bin_edges)
    counts = get_discretization_counts(lindisc, test_statistics)

    return midpoints, counts, bin_width
end

"""
    keep_bottom_proportion_of_hist(midpoints, counts, proportion_to_keep, verbose)

Remove the highest proportion of values from a histogram. Returns the new
midpoints and counts for the truncated histogram.

# Arguments
- `midpoints` : list of midpoints of the histogram bins.
- `counts` : list of counts for the histogram bins.
- `proportion_to_keep` : proportion of lowest values in the histogram to keep.
"""
function keep_bottom_proportion_of_hist(midpoints, counts, proportion_to_keep, verbose)
    @assert 0.0 < proportion_to_keep <= 1.0

    # Number of values to keep
    num_vals_total = sum(counts)
    to_keep = Int(round(proportion_to_keep * num_vals_total))

    # Iterate through counts to find which bins to keep
    # Keep updating total until it is more than we want to keep
    cur_bin = 1
    cur_total = 0
    while (cur_total < to_keep) && (cur_bin <= length(counts))
        cur_total += counts[cur_bin]
        cur_bin += 1
    end

    # Only keep bins whose entire counts are within the kept interval
    final_bin = cur_total == to_keep ? cur_bin-1 : cur_bin-2
    counts = counts[1:final_bin]
    midpoints = midpoints[1:final_bin]

    if verbose
        kept_proportion = sum(counts) / num_vals_total
        println("Kept $kept_proportion of test statistics")
    end

    return midpoints, counts
end

"""
    estimate_fit_parameters(num_bins, bin_width, midpoints, counts, distr::Type{Gamma})

Estimates the poisson parameters required to fit a mode matched null distribution
to the input histogram. Solves the regression problem outlined in equations (6)
and (8) in Schwarzman 2009 (https://projecteuclid.org/euclid.aoas/1231424213).
Returns the parameters C, eta1, eta2
"""
function estimate_fit_parameters(num_bins, bin_width, midpoints, counts, distr::Type{Gamma})
    # Get the log of the counts
    log_counts = log.(counts)

    # Filter out non finite
    log_counts = [isfinite(x) ? x : eps() for x in log_counts]

    # Minus difference h to create regression target
    h = log(num_bins * bin_width)
    y = log_counts .- h

    # Construct and solve the regression problem
    data = DataFrame(Y = y, X1 = midpoints, X2 = log.(midpoints))
    C, eta1, eta2 = coef(glm(@formula(Y ~ X1 + X2), data, Normal(), minstepfac = eps()))

    return C, eta1, eta2
end
"""
    estimate_fit_parameters(num_bins, bin_width, midpoints, counts, distr::Type{Normal})

Estimates the parameters required to fit a mode matched null distribution to the
input histogram. Solves the regression problem outlined in equation (4.5) and (4.6)
in Efron 2007 (https://projecteuclid.org/euclid.aos/1188405614).
Returns the parameters beta0, beta1, beta2
"""
function estimate_fit_parameters(num_bins, bin_width, midpoints, counts, distr::Type{Normal})
    # Get the log of the counts
    log_counts = log.(counts)

    # Filter out non finite
    log_counts = [isfinite(x) ? x : eps() for x in log_counts]

    # Construct and solve the regression problem
    data = DataFrame(Y = log_counts, X1 = midpoints, X2 = midpoints.^2)
    beta0, beta1, beta2 = coef(glm(@formula(Y ~ X1 + X2), data, Normal(), minstepfac = eps()))

    return beta0, beta1, beta2
end

"""
    get_distr_parameters(C, eta1, eta2, distr::Type{Gamma})

Gets the parameters which describe a gamma distribution based on the poisson
parameters input. Converts the poisson parameter inputs into a chi squared
distribution, then converts that into a gamma distribution.
"""
function get_distr_parameters(C, eta1, eta2, distr::Type{Gamma})

    # Convert eta to chi squared parameters
    a = -(1 / (2*eta1))
    v = 2 * (eta2 + 1)

    # Convert chi squared parameters to gamma shape and scale parameters
    k = v/2
    theta = 2*a

    # TODO: Decide what to do about p0
    # Find p0 by solving the following:
    # log(p0) = C + log(gammafunc(eta2 + 1) / (-eta1)^(eta2 + 1))
    # p0 = exp(C + log(gamma(eta2 + 1) / (-eta1)^(eta2 + 1)))
    p0 = 0

    params_valid = k > zero(k) && theta > zero(theta)
    if !params_valid
        error_string = "Error: resulting gamma parameters are invalid. Got k = ", k, " and Θ = ", theta, ".
        Input data might not be gamma distributed, or number of bins needs to change.
        Exiting"
        error(error_string)
    end

    return k, theta, p0
end
"""
    get_distr_parameters(beta0, beta1, beta2, distr::Type{Normal})

Gets the parameters which describe a normal distribution based on the
parameters input. Converts the quadratic parameter inputs into a normal
distribution.
"""
function get_distr_parameters(beta0, beta1, beta2, distr::Type{Normal})

    sigma = (-2 * beta2)^-(1/2)

    mu = beta1 * sigma^2

    # TODO: Decide what to do about p0
    p0 = 0

    return mu, sigma, p0
end

"""
    fit_null_distribution(midpoints, counts, num_bins, bin_width, proportion_to_keep, distr, verbose=true)

Fit a distribution to the discretized inputs. Uses mode-matching, since
the input test statistics may have been truncated.
"""
function fit_null_distribution(midpoints, counts, num_bins, bin_width, proportion_to_keep, distr; verbose = true)

    # Remove highest test statistics, that dont fit the zero assumption
    midpoints, counts = keep_bottom_proportion_of_hist(midpoints, counts, proportion_to_keep, verbose)

    # Estimate parameters using regression
    param0, param1, param2 = estimate_fit_parameters(num_bins, bin_width, midpoints, counts, distr)

    # TODO: Decide what to do about p0
    # Calculate distribution parameters
    k, theta, p0 = get_distr_parameters(param0, param1, param2, distr)

    # Create distribution
    f0 = distr(k, theta)

    return f0
end

"""
    fit_mixture_distribution(midpoints, counts, bin_width)

Fit a cubic spline distribution to the discretized inputs.
"""
function fit_mixture_distribution(midpoints, counts, bin_width)
    # Add some zero bins on the beginning to prevent curve going negative
    number_of_zero_bins = 10 #3
    zero_bins = zeros(number_of_zero_bins)
    pushfirst!(counts, zero_bins...)
    for i in 1 : number_of_zero_bins
        pushfirst!(midpoints, midpoints[1] - bin_width)
    end

    # Normalize counts
    num_statsitics = sum(counts)
    counts = counts / (num_statsitics * bin_width)

    # Fit spline to histogram points
    interpolation = interpolate(counts, BSpline(Cubic(Line(OnCell()))))
    r = range(midpoints[1], stop=midpoints[end], length=length(midpoints))
    scaled_interpolation = Interpolations.scale(interpolation, r)
    scaled_extrapolation = extrapolate(scaled_interpolation, Line())

    fhat(x) = max(scaled_extrapolation(x), 0)
    fhat(x::AbstractArray) = [fhat(i) for i in x]

    return fhat
end

"""
    calculate_posterior(test_statistics, priors, null_distr, mixture_pdf, tail, w0=0.0)

Calculate the empirical Bayes posterior using the priors, null distribution and
mixture distribution.

# Arguments
- `test_statistics` : list of test statistics.
- `priors` : one dimensional list of prior values that correspond with the
   test_statistics, such that `priors[i]` is the prior for `test_statistics[i]`.
- `null_distr` : the null distribution of test statistics
- `mixture_pdf` : function representing the mixture distribution of test statistics
- `tail` : Whether the test is two-tailed (:two) or one-tailed (:lower or :upper)
- `w0` : Default constant for the prior calculation
- `null_value = -Inf` : what to set the unwanted tail values to
"""
function calculate_posterior(test_statistics, priors, null_distr, mixture_pdf, tail; w0=0.0, null_value = -Inf)
    @assert length(test_statistics) == size(priors, 1)

    null_pdf(x) = pdf(null_distr, x)

    # Prior function
    prior_fn(x) = exp(w0) / ( exp(w0) + exp(x) )

    num_test_statistics = length(test_statistics)

    # Posterior array
    posterior = Array{Float64}(undef, num_test_statistics)

    for i in 1:num_test_statistics
        ts = test_statistics[i]
        prior_val = prior_fn(priors[i])
        null_val = null_pdf(ts)
        mix_val = mixture_pdf(ts)

        # if mixture distr equals 0, then just return a 0 posterior
        if mix_val == zero(mix_val)
            posterior[i] = 0.0
            continue
        end

        fdr = null_val / mix_val
        p1 = 1 - prior_val * fdr
        posterior[i] = p1
    end

    # Get rid of the lower tail for an upper-tailed test, or vice-versa
    if tail == :upper || tail == :lower
        null_mode = mode(null_distr)
        compare = tail == :upper ? (<=) : (>=)
        zero_indices = findall(t -> compare(t, null_mode), test_statistics)
        posterior[zero_indices] .= null_value
    end

    return posterior
end

"""
    calculate_posterior(test_statistics, null_distr, mixture_pdf, tail, w0=0.0)

Calculate the empirical Bayes posterior using no priors, null distribution and
mixture distribution.

# Arguments
- `test_statistics` : list of test statistics.
- `null_distr` : the null distribution of test statistics
- `mixture_pdf` : function representing the mixture distribution of test statistics
- `tail` : Whether the test is two-tailed (:two) or one-tailed (:lower or :upper)
- `w0` : Default constant for the prior calculation
- `null_value = -Inf` : what to set the unwanted tail values to
"""
function calculate_posterior(test_statistics, null_distr, mixture_pdf, tail; w0=0.0, null_value = -Inf)
    priors = [0 for _ in test_statistics]
    return calculate_posterior(test_statistics, priors, null_distr, mixture_pdf, tail, w0=w0)
end

"""
    empirical_bayes(test_statistics, priors, num_bins, distr, proportion_to_keep=1.0, tail=:two, w0=0.0)

Calculate the empirical Bayes posteriors of the input statistics using the priors.

# Arguments
- `test_statistics` : list of test statistics.
- `priors` : one dimensional list of prior values that correspond with the
   test_statistics, such that `priors[i]` is the prior for `test_statistics[i]`.
- `num_bins::Integer` : number of uniform width bins to discretize into.
- `distr` : form of the null distribution to be fitted.
- `proportion_to_keep=1.0` : Proportion of lowest test statistics to
   keep when calculating null distribution.
- `tail=:two` : Whether the test is two-tailed (:two) or one-tailed (:lower or :upper)
- `w0` : Default constant for the prior calculation
- `null_value = -Inf` : what to set the unwanted tail values to
"""
function empirical_bayes(test_statistics, priors, num_bins, distr; proportion_to_keep=1.0, tail=:two, w0=0.0, null_value = -Inf)
     midpoints, counts, bin_width = discretize_test_statistics(test_statistics, num_bins)
     null_distr = fit_null_distribution(midpoints, counts, num_bins, bin_width, proportion_to_keep, distr)
     mixture_pdf = fit_mixture_distribution(midpoints, counts, bin_width)
     posteriors = calculate_posterior(test_statistics, priors, null_distr, mixture_pdf, tail, w0=w0, null_value = null_value)
     return posteriors
end
function empirical_bayes(test_statistics, priors, num_bins, distr::Symbol; proportion_to_keep=1.0, tail=:two, w0=0.0, null_value = -Inf)
    return empirical_bayes(test_statistics, priors, num_bins, get_distr(distr); proportion_to_keep=proportion_to_keep, tail=tail, w0=w0, null_value = null_value)
end

"""
    empirical_bayes(test_statistics, num_bins, distr, proportion_to_keep=1.0, tail=:two, w0=0.0)

Calculate the empirical Bayes posteriors of the input statistics with null priors.

# Arguments
- `test_statistics` : list of test statistics.
- `num_bins::Integer` : number of uniform width bins to discretize into.
- `distr` : form of the null distribution to be fitted.
- `proportion_to_keep=1.0` : Proportion of lowest test statistics to
   keep when calculating null distribution.
- `tail=:two` : Whether the test is two-tailed (:two) or one-tailed (:lower or :upper)
- `w0` : Default constant for the prior calculation
- `null_value = -Inf` : what to set the unwanted tail values to
"""
function empirical_bayes(test_statistics, num_bins, distr; proportion_to_keep=1.0, tail=:two, w0=0.0, null_value = -Inf)
    priors = [0 for _ in test_statistics]
    return empirical_bayes(test_statistics, priors, num_bins, distr, proportion_to_keep=proportion_to_keep, tail=tail, w0=w0, null_value = null_value)
end
function empirical_bayes(test_statistics, num_bins, distr::Symbol; proportion_to_keep=1.0, tail=:two, w0=0.0, null_value = -Inf)
    return empirical_bayes(test_statistics, num_bins, get_distr(distr); proportion_to_keep=proportion_to_keep, tail=tail, w0=w0, null_value = null_value)
end

# Utility function for getting a distribution from a symbol, so that the caller needn't
# depend on Distributions.jl
function get_distr(d::Symbol)
    distrs = Dict(
        :Gamma => Gamma,
        :Normal => Normal
    )
    return distrs[d]
end


text:
"""
Module for loading Cartesian ephemeris data
from ASCII files, e.g. plaintext HORIZONS
outputs.
"""
module Ephemeris

export loadascii, interpolator

using Reexport
@reexport using ..AstrodynamicsCore

using DelimitedFiles, Interpolations

include("../Misc/DocStringExtensions.jl")
include("../Misc/UnitfulAliases.jl")

include("LoadASCII.jl")

end


text:
function (F::CNRfForcing)(s::AbstractSimulation)

    cPlaneNumber = 1
    dt = get_dt(s)
    Tf = F.Tf
    ts = Tf

    it = s.iteration[]

    omega = (2*π/(Tf*ts))
    alpha = F.α
    eps = Base.eps()
    _kf = F.Kf
    maxdk = F.maxDk
    avgWaveNumInShell2d = F.avgK
    nShells2d = length(avgWaveNumInShell2d)

    Ef = F.Ef
    R = F.R
    Em = F.Em
    Zf = F.Zf
    factor = F.factor
    N = F.N
    Nm1 = F.Nm1
    dtm2 = s.timestep.x.dt2[]

    s1 = fill!(F.forcex,0)
    s2 = fill!(F.forcey,0)

    u1 = s.u.c.x
    u2 = s.u.c.y

    # // get current spectrum for k=1 plane
    calculate_u1u2_spectrum!(Ef,s.u, cPlaneNumber)

    copyto!(Nm1,N)
    calculate_nn_spectrum!(N,s.u,s.rhs,cPlaneNumber)

    @inbounds for i=2:nShells2d
        if avgWaveNumInShell2d[i] <= _kf
            R[i] += dt*(-2*alpha*omega*(R[i] + N[i]) - (N[i]-Nm1[i])/dtm2 - omega*omega*(Ef[i] - Em[i]))
            R[i]=max(0.0, R[i])
            factor[i] = fsqrt(R[i]/max(Ef[i],eps))*Zf[i]*dt
        end
    end

    umag = 0.0
  #  # -------------------- Horizontal Forcing --------------------
    ff = 0.0
    fv = 0.0

    @inbounds for j=YRANGE
        conjFactX=1.0
        for i=XRANGE
            k = fsqrt(muladd(KX[i],KX[i],KY[j]^2))
            n = round(Int,k/maxdk) + 1
            if (1 < n <= nShells2d)
                if avgWaveNumInShell2d[n] <= _kf
                    s1[i,j,1] = u1[i,j,1]*factor[n]
                    s2[i,j,1] = u2[i,j,1]*factor[n]
                    ff += (abs2(s1[i,j,1]) + abs2(s2[i,j,1]))*conjFactX
                    fv += (proj(u1[i,j,1],s1[i,j,1]) + proj(u2[i,j,1],s2[i,j,1]))*conjFactX
                end
            end
            conjFactX=2.0
        end
    end
    ff = ff*0.5/dt
    fv = fv/dt
    horizontalPower = ff+fv
    F.hp[] = horizontalPower
    powerFraction = 0.99

    # # -------------------- Vertical Forcing --------------------
    # /*
    #   Force the points with kx=ky=0 and kz=1,2,3,4,5 so that the power
    #   added is (1-powerFraction) of that added to the horizontal */
    ff = 0.0
    fv = 0.0
    umag = 0.0
    #     // get the sum of the existing velocites
    @inbounds for k = 2:6
        umag += fsqrt(abs2(u1[1,1,k]) + abs2(u2[1,1,k])) + fsqrt(abs2(u1[1,1,NZ-k+2]) + abs2(u2[1,1,NZ-k+2]))
    end

    f = dt*sqrt(2.0*horizontalPower*(1.0-powerFraction))/umag

    @inbounds for k=2:6
        s1[1,1,k] = u1[1,1,k]*f
        s2[1,1,k] = u2[1,1,k]*f
        s1[1,1,NZ-k+2] = u1[1,1,NZ-k+2]*f
        s2[1,1,NZ-k+2] = u2[1,1,NZ-k+2]*f
        ff += abs2(s1[1,1,k]) + abs2(s2[1,1,k]) + abs2(s1[1,1,NZ-k+2]) + abs2(s2[1,1,NZ-k+2])
        fv += proj(u1[1,1,k],s1[1,1,k]) + proj(u2[1,1,k],s2[1,1,k]) + proj(u1[1,1,NZ-k+2],s1[1,1,NZ-k+2]) + proj(u2[1,1,NZ-k+2],s2[1,1,NZ-k+2])
    end
    ff = ff*0.5/dt
    fv = fv/dt
    verticalPower = ff+fv
    F.vp[] = verticalPower
    return 1/dt
end

function calculate_nn_spectrum!(N,u,rhs,cplane)
    ux = u.c.x
    uy = u.c.y
    rhsx = rhs.c.x
    rhsy = rhs.c.y
    # Initialize the shells to zeros
    fill!(N,0)
    maxdk2d = max(KX[2],KY[2])
    nshells = min(NX,NY÷2)
    @inbounds for j in YRANGE
        conjFactX=1.0
        for i in XRANGE
            kxy2 = KX[i]^2 + KY[j]^2
            k2 = kxy2 + KZ[cplane]^2
            k = fsqrt(kxy2)
            n = round(Int, k/maxdk2d) + 1
            if n <= nshells
                magsq = proj(ux[i,j,cplane],rhsx[i,j,cplane] - ν*k2*ux[i,j,cplane]) + proj(uy[i,j,cplane],rhsy[i,j,cplane] - ν*k2*uy[i,j,cplane])
                ee = conjFactX * magsq;
                N[n]+=ee;
            end
            conjFactX=2.0;
        end
    end
end


text:
"""
    batchrunner(agent_step, model_step, model::AbstractModel, nsteps::Integer, properties::Array{Symbol}, aggregators::AbstractVector{T}, steps_to_collect_data::Array{Int64}, replicates::Integer) where T<:Function

Runs `replicates` number of simulation replicates and returns a `DataFrame`.
"""
function batchrunner(agent_step, model_step, model::AbstractModel, nsteps::Integer, properties::Array{Symbol}, aggregators::AbstractVector{T}, steps_to_collect_data::Array{Int64}, replicates::Integer) where T<:Function
  dataall = step!(agent_step, model_step, model, nsteps, properties, aggregators, steps_to_collect_data)
  for i in 2:replicates
    data = step!(agent_step, model_step, model, nsteps, properties, aggregators, steps_to_collect_data)
    dataall = join(dataall, data, on=:step, kind=:outer, makeunique=true)
  end
  return dataall
end

"""
    batchrunner(agent_step, model_step, model::AbstractModel, nsteps::Integer, properties::Array{Symbol}, steps_to_collect_data::Array{Int64}, replicates::Integer)

Runs `replicates` number of simulation replicates and returns a `DataFrame`.
"""
function batchrunner(agent_step, model_step, model::AbstractModel, nsteps::Integer, properties::Array{Symbol}, steps_to_collect_data::Array{Int64}, replicates::Integer)
  dataall = step!(agent_step, model_step, model, nsteps, properties, steps_to_collect_data)
  for i in 2:replicates
    data = step!(agent_step, model_step, model, nsteps, properties, steps_to_collect_data)
    dataall = join(dataall, data, on=:id, kind=:outer, makeunique=true)
  end
  return dataall
end

"""
    batchrunner(agent_step, model::AbstractModel, nsteps::Integer, properties::Array{Symbol}, aggregators::AbstractVector{T}, steps_to_collect_data::Array{Int64}, replicates::Integer) where T<:Function

Runs `replicates` number of simulation replicates and returns a `DataFrame`.
"""
function batchrunner(agent_step, model::AbstractModel, nsteps::Integer, properties::Array{Symbol}, aggregators::AbstractVector{T}, steps_to_collect_data::Array{Int64}, replicates::Integer) where T<:Function
  dataall = step!(agent_step, model, nsteps, properties, aggregators, steps_to_collect_data)
  for i in 2:replicates
    data = step!(agent_step, model, nsteps, properties, aggregators, steps_to_collect_data)
    dataall = join(dataall, data, on=:step, kind=:outer, makeunique=true)
  end
  return dataall
end

"""
    batchrunner(agent_step, model::AbstractModel, nsteps::Integer, properties::Array{Symbol}, steps_to_collect_data::Array{Int64}, replicates::Integer)

Runs `replicates` number of simulation replicates and returns a `DataFrame`.
"""
function batchrunner(agent_step, model::AbstractModel, nsteps::Integer, properties::Array{Symbol}, steps_to_collect_data::Array{Int64}, replicates::Integer)
  dataall = step!(agent_step, model, nsteps, properties, steps_to_collect_data)
  for i in 2:replicates
    data = step!(agent_step, model, nsteps, properties, steps_to_collect_data)
    dataall = join(dataall, data, on=:id, kind=:outer, makeunique=true)
  end
  return dataall
end

"""
A function to be used in `pmap` in `batchrunner_parallel`. It runs the `step!` function, but has a `dummyvar` parameter that does nothing, but is required for the `pmap` function.
"""
function parallel_step_dummy!(agent_step!::T, model::AbstractModel, nsteps::Integer, agent_properties::Array{Symbol}, steps_to_collect_data::AbstractArray{V}, dummyvar) where {T<:Function, V<:Integer}
  data = step!(agent_step!, model, nsteps, agent_properties, steps_to_collect_data);
  return data
end

function parallel_step_dummy!(agent_step!::T, model::AbstractModel, nsteps::Integer, agent_properties::Array{Symbol}, aggregators::AbstractVector{W}, steps_to_collect_data::AbstractArray{V}, dummyvar) where {T<:Function, V<:Integer, W<:Function}
  data = step!(agent_step!, model, nsteps, agent_properties, aggregators, steps_to_collect_data);
  return data
end

function parallel_step_dummy!(agent_step!::T, model_step!::U, model::AbstractModel, nsteps::Integer, properties::Array{Symbol}, steps_to_collect_data::AbstractArray{V}, dummyvar) where {T<:Function, U<:Function, V<:Integer}
  data = step!(agent_step!, model_step!, deepcopy(model), nsteps, properties, steps_to_collect_data);
  return data
end

function parallel_step_dummy!(agent_step!::T, model_step!::U, model::AbstractModel, nsteps::Integer, properties::Array{Symbol}, aggregators::AbstractVector{W}, steps_to_collect_data::AbstractArray{V}, dummyvar) where {T<:Function, U<:Function, V<:Integer, W<:Function}
  data = step!(agent_step!, model_step!, deepcopy(model), nsteps, properties, aggregators, steps_to_collect_data);
  return data
end

"""
    batchrunner_parallel(agent_step!, model::AbstractModel, nsteps::Integer, agent_properties::Array{Symbol}, steps_to_collect_data::AbstractArray{Integer}, nreplicates::Integer)

Runs `nreplicates` number of simulations in parallel and returns a `DataFrame`.
"""
function batchrunner_parallel(agent_step!::V, model::AbstractModel, nsteps::T, agent_properties::Array{Symbol}, steps_to_collect_data::AbstractArray{T}, nreplicates::T) where {V<:Function, T<:Integer}
  dd = step!(agent_step!, deepcopy(model), nsteps, agent_properties, steps_to_collect_data);
  all_data = pmap(j-> parallel_step_dummy!(agent_step!, model, nsteps, agent_properties, steps_to_collect_data, j), 1:nreplicates)
  for d in all_data
    dd = join(dd, d, on=:id, kind=:outer, makeunique=true)
  end
  return dd
end

"""
    batchrunner_parallel(agent_step!::T, model::AbstractModel, nsteps::V, properties::Array{Symbol}, aggregators::AbstractVector{T}, steps_to_collect_data::AbstractArray{V}, replicates::V) where {T<:Function, V<:Integer}

"""
function batchrunner_parallel(agent_step!::T, model::AbstractModel, nsteps::V, properties::Array{Symbol}, aggregators::AbstractVector{T}, steps_to_collect_data::AbstractArray{V}, replicates::V) where {T<:Function, V<:Integer}
  dd = step!(agent_step!, deepcopy(model), nsteps, properties, aggregators, steps_to_collect_data)
  all_data = pmap(j-> parallel_step_dummy!(agent_step!, model, nsteps, properties, aggregators, steps_to_collect_data, j), 1:nreplicates)
  for d in all_data
    dd = join(dd, d, on=:step, kind=:outer, makeunique=true)
  end
  return dd
end

"""
    batchrunner_parallel(agent_step!::T, model_step!::U, model::AbstractModel, nsteps::V, properties::Array{Symbol}, steps_to_collect_data::AbstractArray, replicates::V) where {T<:Function, U<:Function, V<:Integer}

"""
function batchrunner_parallel(agent_step!::T, model_step!::U, model::AbstractModel, nsteps::V, properties::Array{Symbol}, steps_to_collect_data::AbstractArray, replicates::V) where {T<:Function, U<:Function, V<:Integer}
  dd = step!(agent_step!, model_step!, deepcopy(model), nsteps, properties, steps_to_collect_data)
  all_data = pmap(j-> parallel_step_dummy!(agent_step!, model_step!, model, nsteps, properties, steps_to_collect_data, j), 1:nreplicates)
  for d in all_data
    dd = join(dd, d, on=:id, kind=:outer, makeunique=true)
  end
  return dd
end

"""
    batchrunner(agent_step!::U, model_step!::V, model::AbstractModel, nsteps::Integer, properties::Array{Symbol}, aggregators::AbstractVector{T}, steps_to_collect_data::AbstractArray{X}, replicates::Integer) where {T<:Function,U<:Function,V<:Function,X<:Integer}

"""
function batchrunner(agent_step!::U, model_step!::V, model::AbstractModel, nsteps::Integer, properties::Array{Symbol}, aggregators::AbstractVector{T}, steps_to_collect_data::AbstractArray{X}, replicates::Integer) where {T<:Function,U<:Function,V<:Function,X<:Integer}
  dd = step!(agent_step!, model_step!, deepcopy(model), nsteps, properties, aggregators, steps_to_collect_data)
  all_data = pmap(j-> parallel_step_dummy!(agent_step!, model_step!, model, nsteps, properties, aggregators, steps_to_collect_data, j), 1:nreplicates)
  for d in all_data
    dd = join(dd, d, on=:step, kind=:outer, makeunique=true)
  end
  return dd
end


text:
## Direct calls

const ADJOINT_PARAMETER_COMPATABILITY_MESSAGE =
"""
Adjoint sensitivity analysis functionality requires being able to solve
a differential equation defined by the parameter struct `p`. Thus while
DifferentialEquations.jl can support any parameter struct type, usage
with adjoint sensitivity analysis requires that `p` could be a valid
type for being the initial condition `u0` of an array. This means that
many simple types, such as `Tuple`s and `NamedTuple`s, will work as
parameters in normal contexts but will fail during adjoint differentiation.
To work around this issue for complicated cases like nested structs, look
into defining `p` using `AbstractArray` libraries such as RecursiveArrayTools.jl
or ComponentArrays.jl so that `p` is an `AbstractArray` with a concrete element type.
"""

struct AdjointSensitivityParameterCompatibilityError <: Exception end

function Base.showerror(io::IO, e::AdjointSensitivityParameterCompatibilityError)
  print(io, ADJOINT_PARAMETER_COMPATABILITY_MESSAGE)
end

@doc doc"""
adjoint_sensitivities(sol,alg,g,t=nothing,dg=nothing;
                            abstol=1e-6,reltol=1e-3,
                            checkpoints=sol.t,
                            corfunc_analytical=nothing,
                            callback = nothing,
                            sensealg=InterpolatingAdjoint(),
                            kwargs...)

Adjoint sensitivity analysis is used to find the gradient of the solution
with respect to some functional of the solution. In many cases this is used
in an optimization problem to return the gradient with respect to some cost
function. It is equivalent to "backpropagation" or reverse-mode automatic
differentiation of a differential equation.

Using `adjoint_sensitivities` directly let's you do three things. One it can
allow you to be more efficient, since the sensitivity calculation can be done
directly on a cost function, avoiding the overhead of building the derivative
of the full concretized solution. It can also allow you to be more efficient
by directly controlling the forward solve that is then reversed over. Lastly,
it allows one to define a continuous cost function on the continuous solution,
instead of just at discrete data points.

!!! warning

      Adjoint sensitivity analysis functionality requires being able to solve
      a differential equation defined by the parameter struct `p`. Thus while
      DifferentialEquations.jl can support any parameter struct type, usage
      with adjoint sensitivity analysis requires that `p` could be a valid
      type for being the initial condition `u0` of an array. This means that
      many simple types, such as `Tuple`s and `NamedTuple`s, will work as
      parameters in normal contexts but will fail during adjoint differentiation.
      To work around this issue for complicated cases like nested structs, look
      into defining `p` using `AbstractArray` libraries such as RecursiveArrayTools.jl
      or ComponentArrays.jl so that `p` is an `AbstractArray` with a concrete element type.

!!! warning

      Non-checkpointed InterpolatingAdjoint and QuadratureAdjoint sensealgs
      require that the forward solution `sol(t)` has an accurate dense
      solution unless checkpointing is used. This means that you should
      not use `solve(prob,alg,saveat=ts)` unless checkpointing. If specific
      saving is required, one should solve dense `solve(prob,alg)`, use the
      solution in the adjoint, and then `sol(ts)` interpolate.

### Syntax

There are two forms. For discrete adjoints, the form is:

```julia
du0,dp = adjoint_sensitivities(sol,alg,dg,ts;sensealg=InterpolatingAdjoint(),
                               checkpoints=sol.t,kwargs...)
```

where `alg` is the ODE algorithm to solve the adjoint problem, `dg` is the jump
function, `sensealg` is the sensitivity algorithm, and `ts` is the time points
for data. `dg` is given by:

```julia
dg(out,u,p,t,i)
```

which is the in-place gradient of the cost functional `g` at time point `ts[i]`
with `u=u(t)`.

For continuous functionals, the form is:

```julia
du0,dp = adjoint_sensitivities(sol,alg,g,nothing,(dgdu,dgdp);sensealg=InterpolatingAdjoint(),
                               checkpoints=sol.t,,kwargs...)
```

for the cost functional

```julia
g(u,p,t)
```

with in-place gradient

```julia
dgdu(out,u,p,t)
dgdp(out,u,p,t)
```

If the gradient is omitted, i.e.

```julia
du0,dp = adjoint_sensitivities(sol,alg,g,nothing;kwargs...)
```

then we assume `dgdp` is zero and `dgdu` will be computed automatically using ForwardDiff or finite
differencing, depending on the `autodiff` setting in the `AbstractSensitivityAlgorithm`.
Note that the keyword arguments are passed to the internal ODE solver for
solving the adjoint problem.

### Example discrete adjoints on a cost function

In this example we will show solving for the adjoint sensitivities of a discrete
cost functional. First let's solve the ODE and get a high quality continuous
solution:

```julia
function f(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -p[3]*u[2] + u[1]*u[2]
end

p = [1.5,1.0,3.0]
prob = ODEProblem(f,[1.0;1.0],(0.0,10.0),p)
sol = solve(prob,Vern9(),abstol=1e-10,reltol=1e-10)
```

Now let's calculate the sensitivity of the ``\ell_2`` error against 1 at evenly spaced
points in time, that is:

```math
L(u,p,t)=\sum_{i=1}^{n}\frac{\Vert1-u(t_{i},p)\Vert^{2}}{2}
```

for ``t_i = 0.5i``. This is the assumption that the data is `data[i]=1.0`.
For this function, notice we have that:

```math
\begin{aligned}
dg_{1}&=1-u_{1} \\
dg_{2}&=1-u_{2} \\
& \quad \vdots
\end{aligned}
```

and thus:

```julia
dg(out,u,p,t,i) = (out.=1.0.-u)
```

Also, we can omit `dgdp`, because the cost function doesn't dependent on `p`. If we had data, we'd just replace `1.0` with `data[i]`. To get the adjoint
sensitivities, call:

```julia
ts = 0:0.5:10
res = adjoint_sensitivities(sol,Vern9(),dg,ts,abstol=1e-14,
                            reltol=1e-14)
```

This is super high accuracy. As always, there's a tradeoff between accuracy
and computation time. We can check this almost exactly matches the
autodifferentiation and numerical differentiation results:

```julia
using ForwardDiff,Calculus,Tracker
function G(p)
  tmp_prob = remake(prob,u0=convert.(eltype(p),prob.u0),p=p)
  sol = solve(tmp_prob,Vern9(),abstol=1e-14,reltol=1e-14,saveat=ts,
              sensealg=SensitivityADPassThrough())
  A = convert(Array,sol)
  sum(((1 .- A).^2)./2)
end
G([1.5,1.0,3.0])
res2 = ForwardDiff.gradient(G,[1.5,1.0,3.0])
res3 = Calculus.gradient(G,[1.5,1.0,3.0])
res4 = Tracker.gradient(G,[1.5,1.0,3.0])
res5 = ReverseDiff.gradient(G,[1.5,1.0,3.0])
```

and see this gives the same values.

### Example controlling adjoint method choices and checkpointing

In the previous examples, all calculations were done using the interpolating
method. This maximizes speed but at a cost of requiring a dense `sol`. If it
is not possible to hold a dense forward solution in memory, then one can use
checkpointing. For example:

```julia
ts = [0.0,0.2,0.5,0.7]
sol = solve(prob,Vern9(),saveat=ts)
```

Creates a non-dense solution with checkpoints at `[0.0,0.2,0.5,0.7]`. Now we
can do:

```julia
res = adjoint_sensitivities(sol,Vern9(),dg,ts,
                            sensealg=InterpolatingAdjoint(checkpointing=true))
```

When grabbing a Jacobian value during the backwards solution, it will no longer
interpolate to get the value. Instead, it will start a forward solution at the
nearest checkpoint to build local interpolants in a way that conserves memory.
By default the checkpoints are at `sol.t`, but we can override this:

```julia
res = adjoint_sensitivities(sol,Vern9(),dg,ts,
                            sensealg=InterpolatingAdjoint(checkpointing=true),
                            checkpoints = [0.0,0.5])
```

### Example continuous adjoints on an energy functional

In this case we'd like to calculate the adjoint sensitivity of the scalar energy
functional:

```math
G(u,p)=\int_{0}^{T}\frac{\sum_{i=1}^{n}u_{i}^{2}(t)}{2}dt
```

which is:

```julia
g(u,p,t) = (sum(u).^2) ./ 2
```

Notice that the gradient of this function with respect to the state `u` is:

```julia
function dg(out,u,p,t)
  out[1]= u[1] + u[2]
  out[2]= u[1] + u[2]
end
```

To get the adjoint sensitivities, we call:

```julia
res = adjoint_sensitivities(sol,Vern9(),g,nothing,dg,abstol=1e-8,
                                 reltol=1e-8,iabstol=1e-8,ireltol=1e-8)
```

Notice that we can check this against autodifferentiation and numerical
differentiation as follows:

```julia
using QuadGK
function G(p)
  tmp_prob = remake(prob,p=p)
  sol = solve(tmp_prob,Vern9(),abstol=1e-14,reltol=1e-14)
  res,err = quadgk((t)-> (sum(sol(t)).^2)./2,0.0,10.0,atol=1e-14,rtol=1e-10)
  res
end
res2 = ForwardDiff.gradient(G,[1.5,1.0,3.0])
res3 = Calculus.gradient(G,[1.5,1.0,3.0])
```
"""
function adjoint_sensitivities(sol,args...;
                                  sensealg=InterpolatingAdjoint(),
                                  kwargs...)
  _adjoint_sensitivities(sol,sensealg,args...;kwargs...)
end

function _adjoint_sensitivities(sol,sensealg,alg,g,t=nothing,dg=nothing;
                                   abstol=1e-6,reltol=1e-3,
                                   checkpoints=sol.t,
                                   corfunc_analytical=nothing,
                                   callback = nothing,
                                   kwargs...)

  if !(typeof(sol.prob.p) <: Union{Nothing,SciMLBase.NullParameters,AbstractArray}) || (sol.prob.p isa AbstractArray && !Base.isconcretetype(eltype(sol.prob.p)))
    throw(AdjointSensitivityParameterCompatibilityError())
  end

  if sol.prob isa ODEProblem
    adj_prob = ODEAdjointProblem(sol,sensealg,g,t,dg; checkpoints=checkpoints,
                                 callback = callback,
                                 abstol=abstol,reltol=reltol, kwargs...)

  elseif sol.prob isa SDEProblem
    adj_prob = SDEAdjointProblem(sol,sensealg,g,t,dg,checkpoints=checkpoints,
                                 callback = callback,
                                 abstol=abstol,reltol=reltol,
                                 corfunc_analytical=corfunc_analytical)
  elseif sol.prob isa RODEProblem
    adj_prob = RODEAdjointProblem(sol,sensealg,g,t,dg,checkpoints=checkpoints,
                                callback = callback,
                                abstol=abstol,reltol=reltol,
                                corfunc_analytical=corfunc_analytical)
  else
    error("Continuous adjoint sensitivities are only supported for ODE/SDE/RODE problems.")
  end

  tstops = ischeckpointing(sensealg, sol) ? checkpoints : similar(sol.t, 0)
  adj_sol = solve(adj_prob,alg;
                  save_everystep=false,save_start=false,saveat=eltype(sol[1])[],
                  tstops=tstops,abstol=abstol,reltol=reltol,kwargs...)

  p = sol.prob.p
  l = p === nothing || p === DiffEqBase.NullParameters() ? 0 : length(sol.prob.p)
  du0 = -adj_sol[end][1:length(sol.prob.u0)]

  if eltype(sol.prob.p) <: real(eltype(adj_sol[end]))
    dp = real.(adj_sol[end][(1:l) .+ length(sol.prob.u0)])'
  elseif p === nothing || p === DiffEqBase.NullParameters()
    dp = nothing
  else
    dp = adj_sol[end][(1:l) .+ length(sol.prob.u0)]'
  end

  du0,dp
end

function _adjoint_sensitivities(sol,sensealg::SteadyStateAdjoint,alg,g,dg=nothing;
                                   abstol=1e-6,reltol=1e-3,
                                   kwargs...)
  SteadyStateAdjointProblem(sol,sensealg,g,dg;kwargs...)
end

function _adjoint_sensitivities(sol,sensealg::SteadyStateAdjoint,alg;
                                   g=nothing,dg=nothing,
                                   abstol=1e-6,reltol=1e-3,
                                   kwargs...)
  SteadyStateAdjointProblem(sol,sensealg,g,dg;kwargs...)
end

@doc doc"""
H = second_order_sensitivities(loss,prob,alg,args...;
                               sensealg=ForwardDiffOverAdjoint(InterpolatingAdjoint(autojacvec=ReverseDiffVJP())),
                               kwargs...)

Second order sensitivity analysis is used for the fast calculation of Hessian
matrices.

!!! warning

      Adjoint sensitivity analysis functionality requires being able to solve
      a differential equation defined by the parameter struct `p`. Thus while
      DifferentialEquations.jl can support any parameter struct type, usage
      with adjoint sensitivity analysis requires that `p` could be a valid
      type for being the initial condition `u0` of an array. This means that
      many simple types, such as `Tuple`s and `NamedTuple`s, will work as
      parameters in normal contexts but will fail during adjoint differentiation.
      To work around this issue for complicated cases like nested structs, look
      into defining `p` using `AbstractArray` libraries such as RecursiveArrayTools.jl
      or ComponentArrays.jl so that `p` is an `AbstractArray` with a concrete element type.

### Example second order sensitivity analysis calculation

```julia
using DiffEqSensitivity, OrdinaryDiffEq, ForwardDiff
using Test

function lotka!(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -p[3]*u[2] + p[4]*u[1]*u[2]
end

p = [1.5,1.0,3.0,1.0]; u0 = [1.0;1.0]
prob = ODEProblem(lotka!,u0,(0.0,10.0),p)
loss(sol) = sum(sol)
v = ones(4)

H  = second_order_sensitivities(loss,prob,Vern9(),saveat=0.1,abstol=1e-12,reltol=1e-12)
```

## Arguments

The arguments for this function match `adjoint_sensitivities`. The only notable difference
is `sensealg` which requires a second order sensitivity algorithm, of which currently the
only choice is `ForwardDiffOverAdjoint` which uses forward-over-reverse to mix a forward-mode
sensitivity analysis with an adjoint sensitivity analysis for a faster computation than either
double forward or double reverse. `ForwardDiffOverAdjoint`'s positional argument just accepts
a first order sensitivity algorithm.
"""
function second_order_sensitivities(loss,prob,alg,args...;
                                    sensealg=ForwardDiffOverAdjoint(InterpolatingAdjoint(autojacvec=ReverseDiffVJP())),
                                    kwargs...)
  _second_order_sensitivities(loss,prob,alg,sensealg,args...;kwargs...)
end

@doc doc"""
Hv = second_order_sensitivity_product(loss,v,prob,alg,args...;
                               sensealg=ForwardDiffOverAdjoint(InterpolatingAdjoint(autojacvec=ReverseDiffVJP())),
                               kwargs...)

Second order sensitivity analysis product is used for the fast calculation of
Hessian-vector products ``Hv`` without requiring the construction of the Hessian
matrix.

!!! warning

      Adjoint sensitivity analysis functionality requires being able to solve
      a differential equation defined by the parameter struct `p`. Thus while
      DifferentialEquations.jl can support any parameter struct type, usage
      with adjoint sensitivity analysis requires that `p` could be a valid
      type for being the initial condition `u0` of an array. This means that
      many simple types, such as `Tuple`s and `NamedTuple`s, will work as
      parameters in normal contexts but will fail during adjoint differentiation.
      To work around this issue for complicated cases like nested structs, look
      into defining `p` using `AbstractArray` libraries such as RecursiveArrayTools.jl
      or ComponentArrays.jl so that `p` is an `AbstractArray` with a concrete element type.

### Example second order sensitivity analysis calculation

```julia
using DiffEqSensitivity, OrdinaryDiffEq, ForwardDiff
using Test

function lotka!(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -p[3]*u[2] + p[4]*u[1]*u[2]
end

p = [1.5,1.0,3.0,1.0]; u0 = [1.0;1.0]
prob = ODEProblem(lotka!,u0,(0.0,10.0),p)
loss(sol) = sum(sol)
v = ones(4)

Hv = second_order_sensitivity_product(loss,v,prob,Vern9(),saveat=0.1,abstol=1e-12,reltol=1e-12)
```

## Arguments

The arguments for this function match `adjoint_sensitivities`. The only notable difference
is `sensealg` which requires a second order sensitivity algorithm, of which currently the
only choice is `ForwardDiffOverAdjoint` which uses forward-over-reverse to mix a forward-mode
sensitivity analysis with an adjoint sensitivity analysis for a faster computation than either
double forward or double reverse. `ForwardDiffOverAdjoint`'s positional argument just accepts
a first order sensitivity algorithm.
"""
function second_order_sensitivity_product(loss,v,prob,alg,args...;
                                          sensealg=ForwardDiffOverAdjoint(InterpolatingAdjoint(autojacvec=ReverseDiffVJP())),
                                          kwargs...)
  _second_order_sensitivity_product(loss,v,prob,alg,sensealg,args...;kwargs...)
end


text:
# initialise evotree
function init_evotree(params::EvoTypes{T,U,S},
    X::AbstractMatrix, Y::AbstractVector; verbosity=1) where {T,U,S}

    K = 1
    levels = ""
    X = convert(Matrix{T}, X)
    if typeof(params.loss) == Logistic
        Y = T.(Y)
        μ = fill(logit(mean(Y)), 1)
    elseif typeof(params.loss) == Poisson
        Y = T.(Y)
        μ = fill(log(mean(Y)), 1)
    elseif typeof(params.loss) == Softmax
        if typeof(Y) <: AbstractCategoricalVector
            levels = CategoricalArray(CategoricalArrays.levels(Y))
            K = length(levels)
            μ = zeros(T, K)
            Y = MLJModelInterface.int.(Y)
        else
            levels = CategoricalArray(sort(unique(Y)))
            K = length(levels)
            μ = zeros(T, K)
            Y = UInt32.(Y)
        end
    elseif typeof(params.loss) == Gaussian
        K = 2
        Y = T.(Y)
        μ = SVector{2}([mean(Y), log(std(Y))])
    else
        Y = T.(Y)
        μ = fill(mean(Y), 1)
    end

    # initialize preds
    pred = zeros(SVector{K,T}, size(X,1))
    for i in eachindex(pred)
        pred[i] += μ
    end

    bias = Tree([TreeNode(SVector{K,T}(μ))])
    evotree = GBTree([bias], params, Metric(), K, levels)

    X_size = size(X)
    𝑖_ = collect(1:X_size[1])
    𝑗_ = collect(1:X_size[2])

    # initialize gradients and weights
    δ, δ² = zeros(SVector{evotree.K, T}, X_size[1]), zeros(SVector{evotree.K, T}, X_size[1])
    𝑤 = zeros(SVector{1, T}, X_size[1])
    𝑤_ini = SVector{1, T}(1)
    for i in 1:length(𝑤)
        𝑤[i] += 𝑤_ini
    end

    # binarize data into quantiles
    edges = get_edges(X, params.nbins)
    X_bin = binarize(X, edges)

    # initializde histograms
    hist_δ = Vector{Matrix{SVector{evotree.K, T}}}(undef, 2^params.max_depth-1)
    hist_δ² = Vector{Matrix{SVector{evotree.K, T}}}(undef, 2^params.max_depth-1)
    hist_𝑤 = Vector{Matrix{SVector{1, T}}}(undef, 2^params.max_depth-1)

    # initialize train nodes
    train_nodes = Vector{TrainNode{evotree.K, T, Int64}}(undef, 2^params.max_depth-1)

    for node in 1:2^params.max_depth-1
        train_nodes[node] = TrainNode(0, 0, SVector{evotree.K, T}(fill(T(-Inf), evotree.K)), SVector{evotree.K, T}(fill(T(-Inf), evotree.K)), SVector{1, T}(fill(T(-Inf), 1)), T(-Inf), [0], [0])

        hist_δ[node] = zeros(SVector{evotree.K, T}, params.nbins, X_size[2])
        hist_δ²[node] = zeros(SVector{evotree.K, T}, params.nbins, X_size[2])
        hist_𝑤[node] = zeros(SVector{1, T}, params.nbins, X_size[2])
    end

    splits = Vector{SplitInfo{evotree.K, T, Int64}}(undef, X_size[2])
    for feat in 𝑗_
        splits[feat] = SplitInfo{evotree.K, T, Int}(T(-Inf), SVector{evotree.K, T}(zeros(evotree.K)), SVector{evotree.K, T}(zeros(evotree.K)), SVector{1, T}(zeros(1)), SVector{evotree.K, T}(zeros(evotree.K)), SVector{evotree.K, T}(zeros(evotree.K)), SVector{1, T}(zeros(1)), T(-Inf), T(-Inf), 0, feat, 0.0)
    end

    cache = (params=deepcopy(params),
        X=X, Y=Y, pred=pred,
        𝑖_=𝑖_, 𝑗_=𝑗_, δ=δ, δ²=δ², 𝑤=𝑤,
        edges=edges, X_bin=X_bin,
        train_nodes=train_nodes, splits=splits,
        hist_δ=hist_δ, hist_δ²=hist_δ², hist_𝑤=hist_𝑤)

    cache.params.nrounds = 0

    return evotree, cache
end

function grow_evotree!(evotree::GBTree{L,T,S}, cache; verbosity=1) where {L,T,S}

    # initialize from cache
    params = evotree.params
    train_nodes = cache.train_nodes
    splits = cache.splits
    X_size = size(cache.X_bin)
    δnrounds = params.nrounds - cache.params.nrounds

    # loop over nrounds
    for i in 1:δnrounds

        # select random rows and cols
        𝑖 = cache.𝑖_[sample(params.rng, cache.𝑖_, ceil(Int, params.rowsample * X_size[1]), replace=false, ordered=true)]
        𝑗 = cache.𝑗_[sample(params.rng, cache.𝑗_, ceil(Int, params.colsample * X_size[2]), replace=false, ordered=true)]
        # reset gain to -Inf
        for feat in cache.𝑗_
            splits[feat].gain = T(-Inf)
        end

        # build a new tree
        update_grads!(params.loss, params.α, cache.pred, cache.Y, cache.δ, cache.δ², cache.𝑤)
        ∑δ, ∑δ², ∑𝑤 = sum(cache.δ[𝑖]), sum(cache.δ²[𝑖]), sum(cache.𝑤[𝑖])
        gain = get_gain(params.loss, ∑δ, ∑δ², ∑𝑤, params.λ)
        # assign a root and grow tree
        train_nodes[1] = TrainNode(0, 1, ∑δ, ∑δ², ∑𝑤, gain, 𝑖, 𝑗)
        tree = grow_tree(cache.δ, cache.δ², cache.𝑤, cache.hist_δ, cache.hist_δ², cache.hist_𝑤, params, train_nodes, splits, cache.edges, cache.X_bin)
        push!(evotree.trees, tree)
        predict!(cache.pred, tree, cache.X)

    end #end of nrounds

    cache.params.nrounds = params.nrounds
    # cache = (deepcopy(params), X, Y, pred, 𝑖_, 𝑗_, δ, δ², 𝑤, edges, X_bin, train_nodes, splits, hist_δ, hist_δ², hist_𝑤)
    # return model, cache
    return evotree
end

# grow a single tree
function grow_tree(δ, δ², 𝑤,
    hist_δ, hist_δ², hist_𝑤,
    params::EvoTypes{T,U,S},
    train_nodes::Vector{TrainNode{L,T,S}},
    splits::Vector{SplitInfo{L,T,Int}},
    edges, X_bin) where {T<:AbstractFloat, U, S, L}

    active_id = ones(Int, 1)
    leaf_count = one(Int)
    tree_depth = one(Int)
    tree = Tree(Vector{TreeNode{L, T, Int, Bool}}())

    # grow while there are remaining active nodes
    while size(active_id, 1) > 0 && tree_depth <= params.max_depth
        next_active_id = ones(Int, 0)
        # grow nodes
        for id in active_id
            node = train_nodes[id]
            if tree_depth == params.max_depth || node.∑𝑤[1] <= params.min_weight + 1e-8
                push!(tree.nodes, TreeNode(pred_leaf(params.loss, node, params, δ²)))
            else
                if id > 1 && id == tree.nodes[node.parent].right
                    # println("id is right:", id)
                    hist_δ[id] .= hist_δ[node.parent] .- hist_δ[id-1]
                    hist_δ²[id] .= hist_δ²[node.parent] .- hist_δ²[id-1]
                    hist_𝑤[id] .= hist_𝑤[node.parent] .- hist_𝑤[id-1]
                else
                    # println("id is left:", id)
                    update_hist!(hist_δ[id], hist_δ²[id], hist_𝑤[id], δ, δ², 𝑤, X_bin, node)
                end
                for j in node.𝑗
                    splits[j].gain = node.gain
                    find_split!(view(hist_δ[id],:,j), view(hist_δ²[id],:,j), view(hist_𝑤[id],:,j), params, node, splits[j], edges[j])
                end

                best = get_max_gain(splits)
                # grow node if best split improves gain
                if best.gain > node.gain + params.γ
                    left, right = update_set(node.𝑖, best.𝑖, view(X_bin,:,best.feat))
                    # println("id: ∑𝑤/length(node/left/right) / ", id, " : ", node.∑𝑤, " / ", length(node.𝑖), " / ", length(left), " / ", length(right), " / ", best.𝑖)
                    train_nodes[leaf_count + 1] = TrainNode(id, node.depth + 1, best.∑δL, best.∑δ²L, best.∑𝑤L, best.gainL, left, node.𝑗)
                    train_nodes[leaf_count + 2] = TrainNode(id, node.depth + 1, best.∑δR, best.∑δ²R, best.∑𝑤R, best.gainR, right, node.𝑗)
                    push!(tree.nodes, TreeNode(leaf_count + 1, leaf_count + 2, best.feat, best.cond, best.gain-node.gain, L))
                    push!(next_active_id, leaf_count + 1)
                    push!(next_active_id, leaf_count + 2)
                    leaf_count += 2
                else
                    push!(tree.nodes, TreeNode(pred_leaf(params.loss, node, params, δ²)))
                end # end of single node split search
            end
        end # end of loop over active ids for a given depth
        active_id = next_active_id
        tree_depth += 1
    end # end of tree growth
    return tree
end

# extract the gain value from the vector of best splits and return the split info associated with best split
function get_max_gain(splits::Vector{SplitInfo{L,T,S}}) where {L,T,S}
    gains = (x -> x.gain).(splits)
    feat = findmax(gains)[2]
    best = splits[feat]
    return best
end

function fit_evotree(params, X_train, Y_train;
    X_eval=nothing, Y_eval=nothing,
    early_stopping_rounds=9999,
    eval_every_n=1,
    print_every_n=9999,
    verbosity=1)

    # initialize metric
    iter_since_best = 0
    if params.metric != :none
        metric_track = Metric()
        metric_best = Metric()
    end

    nrounds_max = params.nrounds
    params.nrounds = 0
    model, cache = init_evotree(params, X_train, Y_train)
    iter = 1

    if params.metric != :none && X_eval !== nothing
        pred_eval = predict(model.trees[1], X_eval, model.K)
        Y_eval = convert.(eltype(cache.Y), Y_eval)
    end

    while model.params.nrounds < nrounds_max && iter_since_best < early_stopping_rounds
        model.params.nrounds += 1
        grow_evotree!(model, cache)
        # callback function
        if params.metric != :none
            if X_eval !== nothing
                predict!(pred_eval, model.trees[model.params.nrounds+1], X_eval)
                metric_track.metric = eval_metric(Val{params.metric}(), pred_eval, Y_eval, params.α)
            else
                metric_track.metric = eval_metric(Val{params.metric}(), cache.pred, cache.Y, params.α)
            end
            if metric_track.metric < metric_best.metric
                metric_best.metric = metric_track.metric
                metric_best.iter =  model.params.nrounds
                iter_since_best = 0
            else
                iter_since_best += 1
            end
            if mod(model.params.nrounds, print_every_n) == 0 && verbosity > 0
                display(string("iter:", model.params.nrounds, ", eval: ", metric_track.metric))
            end
        end # end of callback
    end
    if params.metric != :none
        model.metric.iter = metric_best.iter
        model.metric.metric = metric_best.metric
    end
    params.nrounds = nrounds_max
    return model
end


text:
module Eval

using Ekztazy

c = Client()
g = ENV["TESTGUILD"]

command!(c, g, "eval", "evaluates a string as julia code", legacy=false, options=Options(
    [String, "str", "the string to evaluate"]
)) do ctx, str::String
    @info "$str"
    reply(c, ctx, content="```julia\n$(eval(Meta.parse(str)))\n```")
end

start(c)

end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
using StatsBase

sample([2, 3, 5], ProbabilityWeights([0.2, 0.3, 0.5]))

# %%
using StatsBase

a = sample([2, 3, 5], ProbabilityWeights([0.2, 0.3, 0.5]), 10^6)
[(x, count(==(x), a)/length(a)) for x in sort(unique(a))]

# %%
?sample

# %%
using Distributions

rand(DiscreteNonParametric([2, 3, 5], [0.2, 0.3, 0.5]))

# %%
using Distributions

a = rand(DiscreteNonParametric([2, 3, 5], [0.2, 0.3, 0.5]), 10^6)
[(x, count(==(x), a)/length(a)) for x in sort(unique(a))]

# %%
?DiscreteNonParametric

# %%

max_stars_repo_path: 0021/sampling from a set of values.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0021/sampling from a set of values.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0021/sampling from a set of values.jl

max_forks_repo_name: genkuroki/public