text:
using PyPlot
using JLD
using Statistics

# load linear and logarithmic quantisation
D1 = load("/Users/milan/cams/error/linlogroundzfp_all.jld")

Elin = D1["E"][:,1,:]
Elog = D1["E"][:,2,:]
# make compression factors relative to Float64
Elin[:,1] .*= 2
Elog[:,1] .*= 2

D2 = load("/Users/milan/cams/error/round_zstd_all_gridded.jld")
Eround = D2["E"][:,3,:]  # make compression factors relative to unstructured grid
Eround[:,1] .= Eround[:,1]/(900*451)*(348528)*2

@load "/Users/milan/cams/error/roundzfp_all_opt_gridded2.jld"
Ezfp = E[:,5,:]
Ezfp[[33,34,54,55],:] = E[[33,34,54,55],6,:]

Ezfp[:,1] .= Ezfp[:,1]/(900*451)*(348528)*2

## sort into groups
aero = Array(1:15)
ozone = [35,44,45]
methane = [25,26,41]
dynamics = [33,34,54,55,53]
clouds = [22,27,28,31,32,51]
hydro = vcat(Array(36:40),46)
nitro = [42,43,52]
oopp = Array(47:50)
ces = vcat(Array(17:21),[23,24,16])
co12 = [29,30]
chydro = vcat(ces,hydro)

grouped = vcat(aero,co12,clouds,methane,chydro,dynamics,nitro,ozone,oopp)
groups = [aero,co12,clouds,methane,chydro,dynamics,nitro,ozone,oopp]

## geometric mean
cflin = 1/mean(inv,Elin[Elin[:,1].>0,1])
cflog = 1/mean(inv,Elog[Elog[:,1].>0,1])
cfround = 1/mean(inv,Eround[Eround[:,1].>0,1])
# cfround_un = 1/mean(inv,Eround_un[Eround_un[:,1].>0,1])
cfzfp = 1/mean(inv,Ezfp[Ezfp[:,1].>0,1])

## plotting
ioff()
fig,(ax1,ax2) = subplots(1,2,figsize=(10,5),sharey=true)

colours = ["grey","k","C1","C2"]
labls = ["LinQuant24","LogQuant16","Round+lossless","Zfp"]
labls2 = ["Aerosols","Carbon oxides","Clouds & water","Methane","Alkanes,\nalcohols",
            "Dynamics &\ntemperature","N&S oxides","Ozone","Others"]
symblist = ["v","^",">","<","d","s","D","p","h"]

alfa = 0.7

abserri = 4     # statistic to plot, 3 = median, 4 = 90%, 5 = max
decerri = 7     # statistic to plot, 6 = median, 7 = 90%, 8 = max

# ABSOLUTE ERROR
for (i,Es) in enumerate([Elin,Elog,Eround,Ezfp])
    for (ig,g) in enumerate(groups)
        ax1.scatter(Es[g,abserri],log2.(Es[g,1]),30,colours[i],alpha=alfa,
                marker=symblist[ig],edgecolor="k")
    end
end

# DECIMAL ERROR
for (i,Es) in enumerate([Elin,Elog,Eround,Ezfp])
    for (ig,g) in enumerate(groups)
        ax2.scatter(Es[g,decerri],log2.(Es[g,1]),30,colours[i],alpha=alfa,
                    marker=symblist[ig],edgecolor="k")
    end
end

# for varname in ["d"]
#     idx = Array(1:length(varnames))[varnames .== varname]
#     for Es in [Elin,Elog,Eround,Ezfp]
#         ax1.text(Es[idx,abserri],log2.(Es[idx,1]),varname)
#         ax2.text(Es[idx,decerri],log2.(Es[idx,1]),varname)
#     end
# end

ax1.set_yticks([0,1,2,3,4,5,6,7,8,9])
ax1.set_yticks(log2.(Array(4:4:256)),minor=true)
ax1.set_yticklabels([1,2,4,8,16,32,64,128,256,512])
ax1.set_xscale("log")
ax2.set_xscale("log")

ax1.set_xlim(2e-8,9)
ax2.set_xlim(2e-7,99)

ax1.set_ylim(1,log2(256))
ax2.set_ylim(1,log2(256))

ax1.set_title("Absolute error",loc="left")
ax2.set_title("Decimal error",loc="left")

ax1.set_xlabel("norm. absolute error")
ax2.set_xlabel("decimal error")
ax1.set_ylabel("Compression factor")

ax1.set_title("a",loc="right",fontweight="bold")
ax2.set_title("b",loc="right",fontweight="bold")

for i in 1:4
    ax1.fill_betweenx([0,0],[0,0],color=colours[i],label=labls[i],alpha=0.7)
end

for i in 1:9
    ax1.scatter(0,0,30,"white",alpha=alfa,marker=symblist[i],edgecolor="k",
                    label=labls2[i])
end

ax1.legend(loc=2,scatterpoints=3,fontsize=8)

ax1.axhline(log2(cfround),color="C1",zorder=-1)
ax1.axhline(log2(cfzfp),color="C2",zorder=-1)
ax2.axhline(log2(cfround),color="C1",zorder=-1)
ax2.axhline(log2(cfzfp),color="C2",zorder=-1)

ax2.text(1e-6,log2(cfround)+0.1,"$(Int(round(cfround)))x",color="C1",fontweight="bold")
ax2.text(1e-6,log2(cfzfp)+0.1,"$(Int(round(cfzfp)))x",color="C2",fontweight="bold")

tight_layout()
savefig("/Users/milan/git/Elefridge.jl/plots/linlogroundzfp_all.png")
close(fig)


text:
module GenomeFragments

## This is a module for handling high throughput sequencing reads

using GenomicFeatures, DataStructures, BioAlignments, Mmap, UnalignedVectors, XAM, IterTools, ProgressMeter

export FragMatrix, FragMatrixPair, FragMatrixSingle, get_frags, chrom_lt, load_frag_matrix, totalfrags, streambam, qualityfilt, filterpairfirst, filterpairfirstpp,
       barcode, ishighquality, isunique,
       inc_meta_frag!, inc_meta_cut!, inc_meta_mid!,
       inc_heat_frag!, inc_heat_cut!, inc_heat_mid!,
       inc_heat_atac_cut!, inc_heat_atac_frag!, inc_meta_atac_frag!,
       inc_meta_mid_width, inc_heat_mid_width,
       inc_heat_atac_cut_width, inc_meta_atac_cut_width,
       readsintersecting, fragregion,
       fragheatmap, fragheatmap!,
       streambam_singlecell

include("fragmatrix.jl")
include("buildfm.jl")
include("processing.jl")

end # module

ext: jl


text:
module Resonance

export airtable_metadata,
       brain_ingest,
       findprevstool,
       count_set,
       upset_dots!,
       load_metabolites,
       pull_row,
       loadings,
       varexplained,
       mdsaxis,
       permanovas,
       plot_permanovas,
       mantel,
       commonname,
       countmap,
       codebreastfeeding!,
       stp_overlap,
       startup

using Reexport
using Airtable
using CairoMakie
using AlgebraOfGraphics
using MultivariateStats
using Dictionaries
using CodecZlib
using Statistics
using Arrow
using Tables
using SparseArrays
using CodecZlib
using FilePaths
using ProgressLogging
using CategoricalArrays
using PERMANOVA
using LinearAlgebra
using Random

@reexport using BiobakeryUtils
@reexport using DataFrames
@reexport using CSV
@reexport using XLSX
@reexport using DataFramesMeta

const transform = DataFrames.transform

include("airtable.jl")
include("wrangle.jl")
include("metabolites.jl")
include("genefamilies.jl")
include("plotting.jl")
include("brain.jl")
include("startup.jl")
include("kneaddata.jl")
include("mantel.jl")

end

max_stars_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_issues_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_forks_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis


text:
using Schrodinger, PyPlot
using Schrodinger: gaussianpulse

function opt_pihalfx(n)
    Hd = qzero(2)
    Hc = [π*σx, π*σy]
    t = 1
    #ui = rand(n,1) .- 0.5
    ui = zeros(n,2)
    Ut = exp(-1im*π*σx/4)
    # Create objective function type
    O = NormPSU(Ut,Hd,Hc,t,n)
    grape(O,ui)
end

function opt_hadamard(n)
    Hd = σz
    Hc = [σx]
    t = 10.0
    #ui = rand(n,1) .- 0.5
    ui = zeros(n,1)
    Ut = Operator([1/√2 1/√2; 1/√2 -1/√2])
    # Create objective function type
    O = NormPSU(Ut,Hd,Hc,t,n)
    grape(O,ui)
end

function opt_RF_3lvlNOT(n)
    Δ = 2π*(-0.4) # anharmonicity (GHz)
    Hd = Δ*Operator(basis(3,2)) # drift Hamiltonian
    Hc = [create(3)/2+destroy(3)/2, im*create(3)/2-im*destroy(3)/2]
    t = 3 # (ns)
    tvec = LinRange(-t/2,t/2,n)
    ui = [-Δ*gaussianpulse.(tvec,[[t/2,t,0,0,pi]]) LinRange(-Δ/2,Δ/2,n)]
    Ut = Operator([0 1 0; 1 0 0; 0 0 1]) # 3lvl NOT gate
    # Create objective function type
    O = CoherentSubspaces(Ut,1:2,Hd,Hc,t,n) # care only about computational subspace
    grape(O,ui)
end

function opt_3lvlNOT(n)
    ω = 2π*7.0 # 0-1 transition freq (GHz)
    Δ = 2π*(-0.4) # anharmonicity
    Hd = ω*Operator(basis(3,1))+(2ω+Δ)*Operator(basis(3,2)) # drift Hamiltonian
    Hc = [create(3)+destroy(3)]
    t = 3 # (ns)
    tvec = LinRange(-t/2,t/2,n)
    ui = gaussianpulse.(tvec,[[t/2,t,ω,0,pi]])
    #ui = 2 .* cos.(ω.*tvec)
    #Ut = Operator([0 1 0; 1 0 0; 0 0 1]) # 3lvl NOT gate
    Ut = (Operator([0 1 0; 1 0 0; 0 0 0]),Operator([0 0 0; 0 0 0; 0 0 -1]))
    # Create objective function type
    O = CoherentSubspaces(Ut,[1],Hd,Hc,t,n) # care only about computational subspace
    grape(O,ui)
end

function opt_2Q_QFT(n)
    Hd = 0.5 * (σx⊗σx + σy⊗σy + σz⊗σz)
    Hc = [0.5*σx⊗σ0, 0.5*σy⊗σ0, 0.5*σ0⊗σx, 0.5*σ0⊗σy]
    t = 6.0
    #ui = zeros(n,4)
    #ui = rand(n,4) .- 0.5
    ui = Schrodinger.randomsmooth(n,4)
    Ut = Operator([0.5  0.5    0.5  0.5
                   0.5  0.5im -0.5 -0.5im
                   0.5 -0.5    0.5 -0.5
                   0.5 -0.5im -0.5  0.5im],
                  (2,2))
    # Create objective function type
    O = NormPSU(Ut,Hd,Hc,t,n)
    grape(O,ui)
end

benchgrape(O,u,n=100) = for i=1:n; grape(O,u); end

function run_examples()
    r1 = opt_pihalfx(20)
    r2 = opt_hadamard(200)
    r3 = opt_RF_3lvlNOT(300)
    r4 = opt_3lvlNOT(1000)
    r5 = opt_2Q_QFT(500)
    rr = [r1, r2, r3, r4, r5]
    [plotgrape(r) for r in rr]
    return rr
end

function plotgrape(res::Schrodinger.GrapeResult)
    ui = res.ui
    uf = res.uf
    t = 0 : res.t/size(uf,1) : res.t
    figure()
    step(t,[uf[1:1,:]; uf])
    gca()[:set_prop_cycle](nothing)
    step(t,[ui[1:1,:]; ui],linestyle="dashed")
    legend(["Control $i" for i in 1:size(uf,2)])
    tight_layout(true)
    grid(true)
end


text:
include("C:\\Users\\warner\\Desktop\\My Stuff\\Personal Learning\\Monte Carlo\\NautilusMC\\NautilusMC\\Materials.jl")
include("C:\\Users\\warner\\Desktop\\My Stuff\\Personal Learning\\Monte Carlo\\NautilusMC\\NautilusMC\\geometry.jl")

water = define_Material(
        density = 1.0,
        nuclides = ["O-nat", "H-nat"],
        atomic_fracs = [1.0, 2.0],
        color = [0, 200, 0])

B4C = define_Material(
        density = 2.52,
        nuclides = ["B-nat", "C-nat"],
        atomic_fracs = [4.0, 1.0],
        color = [100, 0, 100])

s1 = box([-1.0, -1.0, -1.0], [1.0, 1.0, 1.0])
s2 = box([-2.0, -2.0, -2.0], [2.0, 2.0, 2.0])
s3 = box([-3.0, -3.0, -3.0], [3.0, 3.0, 3.0])

C1 = Cell([s1], [], water, [])
C2 = Cell([s2], [s1], B4C, [])
outside = Cell([], [s2], nothing, [])

Root = Universe([C1, C2, outside], [])

max_stars_repo_name: wamcghee/NautilusMC

max_issues_repo_name: wamcghee/NautilusMC

max_forks_repo_name: wamcghee/NautilusMC


text:
module InfDecomp_Base

export My_Eval, condEntropy, create_My_Eval, Set_Data, create_stuff

using MathProgBase
# http://mathprogbasejl.readthedocs.io/en/latest/nlp.html

type My_Eval <: MathProgBase.AbstractNLPEvaluator
    n_x     :: Int64
    n_y     :: Int64
    n_z     :: Int64

    n     :: Int64         # number of variables
    m     :: Int64         # number of marginal equations

    varidx  :: Array{Int64,3}                   # 0 if variable not present; otherwise idx of var in {1,...,n}
    xyz     :: Vector{Tuple{Int64,Int64,Int64}} # xyz-triple of varidx

    eqidx   :: Dict{ Tuple{String,Int64,Int64},  Int64} # first idx is "xy" or "xz". Index of eqn in {1,...,m}
    mr_eq   :: Vector{ Tuple{String,Int64,Int64} }      # ("xy",x,y) / ("yz",y,z) of an eqn
    prb_xyz :: Array{Float64,3}
    marg_xy :: Array{Float64,2}
    marg_xz :: Array{Float64,2}

    count_eval_f      :: Int64
    count_eval_grad_f :: Int64
    count_eval_g      :: Int64
    count_eval_jac_g  :: Int64
    count_hesseval    :: Int64

    rhs   :: Vector{Float64}                # m-vector
    Gt    :: SparseMatrixCSC{Float64,Int64} # G^T n x m; transpose of constraint matrix
    Gt_K  :: Vector{Int64}                  # findn()-info of G^T: rows (vars)
    Gt_L  :: Vector{Int64}                  # findn()-info of G^T: columns (equations)

    TmpFloat       :: DataType
    bigfloat_nbits :: Int64
end

function condEntropy{TFloat_2,TFloat}(e::My_Eval, p::Vector{TFloat_2}, dummy::TFloat)   :: TFloat
    # m_yz = marg_yz(e,p,dummy)
    s::TFloat = TFloat(0.)
    for y = 1:e.n_y
        for z = 1:e.n_z
            P_yz = TFloat(0.)
            for x = 1:e.n_x
                i = e.varidx[x,y,z]
                if i>0
                    P_yz += p[i]
                end
            end
            for x = 1:e.n_x
                i = e.varidx[x,y,z]
                if i>0
                    p_xyz = p[i]
                    s  +=  (  (p_xyz ≤ 0 || P_yz  ≤ 0)  ?   TFloat(0.)   :   - p_xyz*log( p_xyz / P_yz )   )
                end
            end
        end
    end
    return s
    ;
end

function ∇f{TFloat,TFloat_2}(e::My_Eval, grad::Vector{TFloat_2}, p::Vector{TFloat_2}, dummy::TFloat) :: Void
    for y = 1:e.n_y
        for z = 1:e.n_z
            # make marginal P(*yz)
            P_yz::TFloat = TFloat(0.)
            for x = 1:e.n_x
                i = e.varidx[x,y,z]
                if i>0
                    P_yz += p[i]
                end
            end
            # make log-expressions  log( P(xyz) / P(*yz) )
            for x = 1:e.n_x
                i = e.varidx[x,y,z]
                if i>0
                    P_xyz::TFloat = TFloat( p[i] )
                    grad[i] = TFloat_2(   (P_xyz ≤ 0 || P_yz ≤ 0) ?  -log(TFloat(e.n_x))  : log( P_xyz / P_yz )  )
                end
            end
        end# for y
    end# for x
    ;
end # ∇f()

function create_My_Eval(q::Array{Float64,3}, tmpFloat::DataType, bigfloat_precision=256)
    const n_x::Int64 = size(q,1);
    const n_y::Int64 = size(q,2);
    const n_z::Int64 = size(q,3);

    # Create marginals
    prb_xyz::Array{Float64,3}  = zeros(n_x,n_y,n_z)
    marg_xy::Array{Float64,2}  = zeros(n_x,n_y)
    marg_xz::Array{Float64,2}  = zeros(n_x,n_z)

    count_eval_f      :: Int64 = 0
    count_eval_grad_f :: Int64 = 0
    count_eval_g      :: Int64 = 0
    count_eval_jac_g  :: Int64 = 0
    count_hesseval    :: Int64 = 0
    for x in 1:n_x
        for y in 1:n_y
            for z in 1:n_z
                marg_xy[x,y] += q[x,y,z]
                marg_xz[x,z] += q[x,y,z]
            end
        end
    end
    # Find the variables
    varidx::Array{Int64,3}                   = zeros(Bool,size(q));
    xyz   ::Vector{Tuple{Int64,Int64,Int64}} = [ (0,0,0) for i in 1:n_x*n_y*n_z ]
    n::Int64 = 0
    for x in 1:n_x
        for y in 1:n_y
            for z in 1:n_z
                if marg_xy[x,y] > 0  &&  marg_xz[x,z] > 0
                    n += 1
                    varidx[x,y,z] = n
                    xyz[n]        = (x,y,z)
                    prb_xyz[x,y,z]= q[x,y,z]
                else
                    varidx[x,y,z] = 0
                end#if
            end
        end
    end

    # Find the equations
    eqidx = Dict{ Tuple{String,Int64,Int64},Int64}() # first idx is "xy" or "xz"
    mr_eq ::Vector{ Tuple{String,Int64,Int64} }   = [ ("",0,0)   for i in 1:n_x*(n_y+n_z) ]
    m::Int64 = 0
    for x in 1:n_x
        for y in 1:n_y
            if marg_xy[x,y] > 0
                m += 1
                eqidx["xy",x,y] = m
                mr_eq[m]        = ("xy",x,y)
            else
                eqidx["xy",x,y] = 0
            end#if
        end
        for z in 1:n_z
            if marg_xz[x,z] > 0
                m += 1
                eqidx["xz",x,z] = m
                mr_eq[m]        = ("xz",x,z)
            else
                eqidx["xz",x,z] = 0
            end#if
        end
    end #for x

    rhs::Vector{Float64} = zeros(m)
    for k in 1:m
        mr = mr_eq[k]
        if mr[1] == "xy"
            rhs[k] = marg_xy[ mr[2], mr[3] ]
        elseif mr[1]=="xz"
            rhs[k] = marg_xz[ mr[2], mr[3] ]
        else
            print("Fuck! Bug!")
            return;
        end
    end #for all marg equations

    denseGt :: Array{Float64,2} = zeros(n,m)
    for l in 1:n
        (x,y,z) = xyz[l]
        for k in 1:m
            mr = mr_eq[k]
            if mr[1] == "xy"
                xy = mr[2:3]
                if xy[1]==x && xy[2]==y
                    denseGt[l,k] = 1.
                end
            elseif mr[1]=="xz"
                xz = mr[2:3]
                if xz[1]==x && xz[2]==z
                    denseGt[l,k] = 1.
                end
            else
                print("Fuck! Bug!")
                return;
            end
        end
    end

    Gt::SparseMatrixCSC{Float64,Int64} = sparse(denseGt)
    local Gt_K::Array{Int64,1}
    local Gt_L::Array{Int64,1}
    (Gt_K,Gt_L) = findn(Gt)

    TmpFloat       :: DataType  = tmpFloat
    bigfloat_nbits :: Int64     = bigfloat_precision

    return My_Eval(n_x,n_y,n_z, n,m, varidx,xyz, eqidx,mr_eq, prb_xyz, marg_xy,marg_xz, count_eval_f, count_eval_grad_f, count_eval_g, count_eval_jac_g, count_hesseval, rhs, Gt,Gt_K,Gt_L,  TmpFloat,bigfloat_nbits)
    ;
end #^ create_My_eval()

type Set_Data{T1,T2,T3}
    Xidx :: Dict{T1,Int64}
    Yidx :: Dict{T2,Int64}
    Zidx :: Dict{T3,Int64}
    X    :: Vector{T1}
    Y    :: Vector{T2}
    Z    :: Vector{T3}
end

function create_Set_Data{T1,T2,T3}(pdf::Dict{Tuple{T1,T2,T3},Float64}) :: Set_Data{T1,T2,T3}
    Xidx = Dict{T1,Int64}()
    Yidx = Dict{T2,Int64}()
    Zidx = Dict{T3,Int64}()

    anX=aY=aZ=nothing
    for (xyz,val) in pdf
        x,y,z = xyz
        if x ∉ keys(Xidx)
            Xidx[x] = length(Xidx)+1
            anX = x
        end
        if y ∉ keys(Yidx)
            Yidx[y] = length(Yidx)+1
            aY = y
        end
        if z ∉ keys(Zidx)
            Zidx[z] = length(Zidx)+1
            aZ = z
        end
    end

    X ::Vector{T1} = [ anX  for i = 1:length(Xidx) ]
    Y ::Vector{T2} = [ aY   for i = 1:length(Yidx) ]
    Z ::Vector{T3} = [ aZ   for i = 1:length(Zidx) ]

    for (x,i) in Xidx
        X[i] = x
    end
    for (y,i) in Yidx
        Y[i] = y
    end
    for (z,i) in Zidx
        Z[i] = z
    end

    @assert ( length(X) ≥ 2 ) "|Range(X)| ≥ 2 needed"
    @assert ( length(Y) ≥ 2 ) "|Range(Y)| ≥ 2 needed"
    @assert ( length(Z) ≥ 2 ) "|Range(Z)| ≥ 2 needed"

    return Set_Data{T1,T2,T3}(Xidx,Yidx,Zidx,X,Y,Z)
    ;
end #^ create_Set_Data()

function create_stuff{T1,T2,T3}(pdf::Dict{Tuple{T1,T2,T3},Float64}, tmpFloat::DataType) :: Tuple{ Set_Data{T1,T2,T3}, My_Eval }
    sd = create_Set_Data(pdf)

    Q = zeros(Float64,length(sd.X),length(sd.Y),length(sd.Z))
    for (xyz,val) in pdf
        Q[ sd.Xidx[xyz[1]], sd.Yidx[xyz[2]], sd.Zidx[xyz[3]] ] = val
    end

    e = create_My_Eval(Q,tmpFloat)

    return (sd,e)
    ;
end #^ create_stuff

end #^ module InfDecomp_base
;#EOF

max_stars_repo_name: dot-at/newpid

max_issues_repo_name: dot-at/newpid

max_forks_repo_name: dot-at/newpid


text:
####
#### Defines list of how-to-guides
####

how_to_guides = Any[
    "Balance Laws" => Any["How to make a balance law" => "HowToGuides/BalanceLaws/how_to_make_a_balance_law.md",],
    "Atmos" => Any[
        "Reference profiles" => "HowToGuides/Atmos/AtmosReferenceState.md",
        "Moisture model" => "HowToGuides/Atmos/MoistureModelChoices.md",
        "Precipitation model" => "HowToGuides/Atmos/PrecipitationModelChoices.md",
        "How to create an experiment with moisture and precipitation" => "HowToGuides/Atmos/MoistureAndPrecip.md",
    ],
    "Ocean" => Any[
    # "Home" => "HowToGuides/Ocean/index.md"
    ],
    "Land" => Any[
    # "Home" => "HowToGuides/Land/index.md"
    ],
    "Numerics" => Any[
        "Meshes" => Any[
        # "Home" => "HowToGuides/Numerics/Meshes/index.md",
        ],
        "ODE Solvers" => Any["Time-integration" => "HowToGuides/Numerics/ODESolvers/Timestepping.md",],
        "System Solvers" => Any["Iterative Solvers" => "HowToGuides/Numerics/SystemSolvers/IterativeSolvers.md",],
    ],
    "Diagnostics" => Any["Using Diagnostics" => "HowToGuides/Diagnostics/UsingDiagnostics.md",],
]


text:
include("lib/utils.jl")
include("lib/json.jl")

struct Token
  lineno::Int
  kind::String
  value::String
end

tokens = Token[]
pos = 1

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

function peek(offset = 0)
  tokens[pos + offset]
end

function consume(str)
  if peek().value == str
    global pos += 1
  else
    error("unexpected token ($( peek() ))")
  end
end

function inc_pos()
  global pos += 1
end

function is_end()
  length(tokens) < pos
end

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

function parse_arg()
  t = peek()
  inc_pos()

  if t.kind == "ident"
    t.value
  elseif t.kind == "int"
    parse(Int, t.value)
  else
    error("unexpected kind $(t)")
  end
end

function parse_args()
  args = []

  if peek().value == ")"
    return args
  end

  push!(args, parse_arg())

  while peek().value == ","
    consume(",")
    push!(args, parse_arg())
  end

  args
end

function parse_expr_right()
  t = peek()

  if t.value == "+"
    consume("+")
    expr_r = parse_expr()
    ["+", expr_r]
  elseif t.value == "*"
    consume("*")
    expr_r = parse_expr()
    ["*", expr_r]
  elseif t.value == "=="
    consume("==")
    expr_r = parse_expr()
    ["eq", expr_r]
  elseif t.value == "!="
    consume("!=")
    expr_r = parse_expr()
    ["neq", expr_r]
  else
    error("unexpected token")
  end
end

function next_binop()
  val = peek().value

  val == "+" ||
  val == "*" ||
  val == "==" ||
  val == "!="
end

function parse_expr()
  t_left = peek()

  if t_left.value == "("
    consume("(")
    expr_l = parse_expr()
    consume(")")
  else
    inc_pos()
    if t_left.kind == "int"
      expr_l = parse(Int, t_left.value)
    elseif t_left.kind == "ident"
      expr_l = t_left.value
    else
      error("unexpected token")
    end
  end

  if next_binop()
    op, expr_r = parse_expr_right()
    [op, expr_l, expr_r]
  else
    expr_l
  end
end

function parse_set()
  consume("set")

  var_name = peek().value
  inc_pos()

  consume("=")

  expr = parse_expr()

  consume(";")

  ["set", var_name, expr]
end

function _parse_funcall()
  fn_name = peek().value
  inc_pos()

  consume("(")
  args = parse_args()
  consume(")")

  vcat([fn_name], args)
end

function parse_call()
  consume("call")

  funcall = _parse_funcall()

  consume(";")

  vcat(["call"], funcall)
end

function parse_call_set()
  consume("call_set")

  var_name = peek().value
  inc_pos()

  consume("=")
  funcall = _parse_funcall()
  consume(";")

  ["call_set", var_name, funcall]
end

function parse_return()
  consume("return")

  expr = parse_expr()

  consume(";")

  ["return", expr]
end

function parse_while()
  consume("while")

  consume("(")
  expr = parse_expr()
  consume(")")

  consume("{")
  stmts = parse_stmts()
  consume("}")

  ["while", expr, stmts]
end

function _parse_when_clause()
  consume("(")
  expr = parse_expr()
  consume(")")

  consume("{")
  stmts = parse_stmts()
  consume("}")

  vcat([expr], stmts)
end

function parse_case()
  consume("case")
  consume("{")

  when_clauses = Any[]

  while peek().value != "}"
    push!(when_clauses, _parse_when_clause())
  end

  consume("}")

  vcat(["case"], when_clauses)
end

function parse_vm_comment()
  consume("_cmt")
  consume("(")

  vm_comment = peek().value
  inc_pos()

  consume(")")
  consume(";")

  ["_cmt", vm_comment]
end

function parse_stmt()
  t = peek()

  if t.value == "set"
    parse_set()
  elseif t.value == "call"
    parse_call()
  elseif t.value == "call_set"
    parse_call_set()
  elseif t.value == "return"
    parse_return()
  elseif t.value == "while"
    parse_while()
  elseif t.value == "case"
    parse_case()
  elseif t.value == "_cmt"
    parse_vm_comment()
  else
    error("unexpected token ($(t))")
  end
end

function parse_stmts()
  stmts = []

  while peek().value != "}"
    push!(stmts, parse_stmt())
  end

  stmts
end

function parse_var_declare()
  var_name = peek().value
  inc_pos()

  consume(";")

  ["var", var_name]
end

function parse_var_init()
  var_name = peek().value
  inc_pos()

  consume("=")

  expr = parse_expr()

  consume(";")

  ["var", var_name, expr]
end

function parse_var()
  consume("var")

  t1 = peek(1)
  if t1.value == ";"
    parse_var_declare()
  elseif t1.value == "="
    parse_var_init()
  else
    error("unexpected token")
  end
end

function parse_func()
  consume("func")

  fn_name = peek().value
  inc_pos()

  consume("(")
  args = parse_args()
  consume(")")

  consume("{")

  stmts = Any[]

  while peek().value != "}"
    if peek().value == "var"
      push!(stmts, parse_var())
    else
      push!(stmts, parse_stmt())
    end
  end

  consume("}")

  ["func", fn_name, args, stmts]
end

function parse_top_stmt()
  parse_func()
end

function parse_top_stmts()
  top_stmts = Any["top_stmts"]

  while ! is_end()
    push!(top_stmts, parse_top_stmt())
  end

  top_stmts
end

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

for line in readlines()
  list = json_parse(line)
  t = Token(
    list[1],
    list[2],
    list[3]
  )
  push!(tokens, t)
end

ast = parse_top_stmts()
json_print(ast, true)


text:
@doc raw"""
    SpecialOrthogonal{n} <: GroupManifold{ℝ,Rotations{n},MultiplicationOperation}

Special orthogonal group $\mathrm{SO}(n)$ represented by rotation matrices.

# Constructor
    SpecialOrthogonal(n)
"""
const SpecialOrthogonal{n} = GroupManifold{ℝ,Rotations{n},MultiplicationOperation}

invariant_metric_dispatch(::SpecialOrthogonal, ::ActionDirection) = Val(true)

function default_metric_dispatch(
    ::MetricManifold{𝔽,<:SpecialOrthogonal,EuclideanMetric},
) where {𝔽}
    return Val(true)
end
default_metric_dispatch(::SpecialOrthogonal, ::EuclideanMetric) = Val(true)

SpecialOrthogonal(n) = SpecialOrthogonal{n}(Rotations(n), MultiplicationOperation())

Base.show(io::IO, ::SpecialOrthogonal{n}) where {n} = print(io, "SpecialOrthogonal($(n))")

Base.inv(::SpecialOrthogonal, p) = transpose(p)
Base.inv(::SpecialOrthogonal, e::Identity{MultiplicationOperation}) = e

inverse_translate(G::SpecialOrthogonal, p, q, ::LeftAction) = inv(G, p) * q
inverse_translate(G::SpecialOrthogonal, p, q, ::RightAction) = q * inv(G, p)

translate_diff(::SpecialOrthogonal, p, q, X, ::LeftAction) = X
translate_diff(G::SpecialOrthogonal, p, q, X, ::RightAction) = inv(G, p) * X * p

function translate_diff!(G::SpecialOrthogonal, Y, p, q, X, conv::ActionDirection)
    return copyto!(Y, translate_diff(G, p, q, X, conv))
end

function inverse_translate_diff(G::SpecialOrthogonal, p, q, X, conv::ActionDirection)
    return translate_diff(G, inv(G, p), q, X, conv)
end

function inverse_translate_diff!(G::SpecialOrthogonal, Y, p, q, X, conv::ActionDirection)
    return copyto!(Y, inverse_translate_diff(G, p, q, X, conv))
end

function allocate_result(
    ::GT,
    ::typeof(exp),
    ::Identity,
    X,
) where {n,GT<:SpecialOrthogonal{n}}
    return allocate(X)
end
function allocate_result(
    ::GT,
    ::typeof(log),
    ::Identity,
    q,
) where {n,GT<:SpecialOrthogonal{n}}
    return allocate(q)
end


text:
using Distributed
addprocs(2)
@everywhere using DiffEqGPU, CUDA, OrdinaryDiffEq, Test, Random

@everywhere begin
    function lorenz_distributed(du,u,p,t)
        du[1] = p[1]*(u[2]-u[1])
        du[2] = u[1]*(p[2]-u[3]) - u[2]
        du[3] = u[1]*u[2] - p[3]*u[3]
    end
    CUDA.allowscalar(false)
    u0 = Float32[1.0;0.0;0.0]
    tspan = (0.0f0,100.0f0)
    p = (10.0f0,28.0f0,8/3f0)
    Random.seed!(1)
    pre_p_distributed = [rand(Float32,3) for i in 1:10]
    function prob_func_distributed(prob,i,repeat)
        remake(prob,p=pre_p_distributed[i].*p)
    end
end

prob = ODEProblem(lorenz_distributed,u0,tspan,p)
monteprob = EnsembleProblem(prob, prob_func = prob_func_distributed)

#Performance check with nvvp
# CUDAnative.CUDAdrv.@profile
@time sol = solve(monteprob,Tsit5(),EnsembleGPUArray(),trajectories=10,saveat=1.0f0)
@test length(filter(x -> x.u != sol.u[1].u, sol.u)) != 0 # 0 element array
@time sol = solve(monteprob,ROCK4(),EnsembleGPUArray(),trajectories=10,saveat=1.0f0)
@time sol2 = solve(monteprob,Tsit5(),EnsembleGPUArray(),trajectories=10,
                                                 batch_size=5,saveat=1.0f0)

@test length(filter(x -> x.u != sol.u[1].u, sol.u)) != 0 # 0 element array
@test length(filter(x -> x.u != sol2.u[6].u, sol.u)) != 0 # 0 element array
@test all(all(sol[i].prob.p .== pre_p_distributed[i].*p) for i in 1:10)
@test all(all(sol2[i].prob.p .== pre_p_distributed[i].*p) for i in 1:10)

#To set 1 GPU per device:
#=
using Distributed
addprocs(numgpus)
import CUDAdrv, CUDAnative

let gpuworkers = asyncmap(collect(zip(workers(), CUDAdrv.devices()))) do (p, d)
  remotecall_wait(CUDAnative.device!, p, d)
  p
end
=#


text:
## includes.jl.  Included by both ROC.jl and nomodule.jl
## (c) 2015 David A. van Leeuwen.
##
## Licensed under the MIT software license, see LICENSE.md

include("main.jl")
include("eer.jl")
include("cllr.jl")
include("auc.jl")
include("pav.jl")
include("dcf.jl")

include("plotrecipes.jl")

include("dataframes.jl")

include("delong.jl")


text:
module TestLBV

########################### submodule Rbd ###################

module Airframe
using Flight.LBV
export XAirframe

@define_node XAirframe (att = LBVLeaf{4}, vel = LBVLeaf{3}, pos = LBVLeaf{4})

end #submodule

#################### submodule Ldg ########################

module Ldg
using Flight.LBV
export XLdg

@define_node XLdg (nlg = LBVLeaf{2}, mlg = LBVLeaf{2})

end #submodule

#################### submodule Aircraft ########################

module Aircraft
using Flight.LBV
using ..Airframe #needed to access XAirframe
using ..Ldg #needed to access XLdg
export XAircraft

@define_node XAircraft (rbd = XAirframe, ldg = XLdg, pwp = LBVLeaf{4})

end #submodule

#################### Tests ########################

using Flight.LBV
using Reexport
@reexport using .Airframe
@reexport using .Ldg
@reexport using .Aircraft

export test_lbv

function test_lbv()

    #test that when broadcasting, the result is of the same type and numerically
    #correct
    #test that when broadcasting, numbers are promoted to the appropriate type
    #test that we can construct from a vector or from a view
    #when we extract a field, it is a view and it is equal to the appropriate
    #indices
    #when we change the field, it changes the original
    #when we extract a LBVLeaf, it is a Leaf
    #reject inappropriate length
    #test copy
    #test similar
    x = XAircraft()
    x_rbd = x.rbd
    x_rbd[4] = 0
    @show x
    v = XAircraft(view(x,:))
    v.pwp .= 3
    @show x
    #broadcasting with LBVLeafs of mixed parametric subtypes
    xv_rbd = x.pwp .+ v.pwp
    #broadcasting with LBVNodes of mixed parametric subtypes
    xv = x .+ v
    y = similar(x)
    @. y = x + sin(x) + 2x
    @show y

    z1 = XAircraft(zeros(Int32, length(XAircraft)))
    z2 = XAircraft(ones(Int64, length(XAircraft)))
    z = exp.(z1) + z2
    @show z

    # println(Aircraft.descriptor(z))

end

end #module


text:
# image_geom.jl

using MIRT: ImageGeom, image_geom, cbct
using MIRT: image_geom_ellipse
using MIRT: jim
using FillArrays: Trues

using Test: @test, @testset, @test_throws, @inferred

function image_geom_test2(ig::ImageGeom)
	# test 2D functions provided by the constructor
	ig.dim
	ig.dims
	ig.deltas
	ig.offsets
	ig.x
	ig.y
	ig.wx
	ig.wy
	ig.xg
	ig.yg
	ig.fovs
	ig.np
	ig.mask_outline
	ig.ones
	ig.zeros
	ig.u
	ig.v
	ig.ug
	ig.vg
	ig.fg
	@test ig.shape(vec(ig.ones)) == ig.ones
	@test ig.embed(ig.ones[ig.mask]) == Float32.(ig.mask)
	@test ig.maskit(ig.ones) == Float32.(ones(ig.np))
#	@inferred # todo
	ig.unitv()
	ig.unitv(j=4)
	ig.unitv(i=ntuple(i->1, length(ig.dim)))
	ig.unitv(c=ntuple(i->0, length(ig.dim)))
#= todo-i: why do these fail?
	@inferred image_geom_ellipse(8, 10, 1, 2)
	@inferred ig.circ()
	@inferred ig.plot()
	@inferred ig.unitv()
=#
	ig.circ()
	ig.plot(jim)
#	isinteractive() && gui()
	ig.unitv()
	@inferred ig.down(2)
	@inferred ig.over(2)
	image_geom_ellipse(8, 10, 1, 2)
	true
end

function image_geom_test2()
	#@inferred
	image_geom(nx=16, dx=2, offsets=:dsp, mask=:all_but_edge_xy)
	@test_throws String image_geom(nx=16, dx=1, offsets=:bad)
	@test_throws String image_geom(nx=16, dx=1, mask=:bad) # mask type
	@test_throws String image_geom(nx=16, dx=1, mask=trues(2,2)) # mask size
	ig = image_geom(nx=16, dx=2)
	@test ig.mask == Trues(ig.dims)
	show(isinteractive() ? stdout : devnull, ig)
	show(isinteractive() ? stdout : devnull, MIME("text/plain"), ig)
	image_geom_test2(ig)
	ig = image_geom(nx=16, dx=2, mask=:circ)
	ig.over(2)
	ig.down(3) # test both even and non-even factors
	ig.help
	true
end

function image_geom_test3(ig::ImageGeom)
	@test image_geom_test2(ig)
	ig.wz
	ig.zg
	ig.mask_or
	@inferred ig.expand_nz(2)
	@inferred cbct(ig)
	@test_throws String image_geom(nx=1, dx=1, nz=2, offset_z=-1, offsets=:dsp)
	true
end

function image_geom_test3()
	#@inferred
	image_geom(nx=16, nz=4, dx=2, zfov=1)
	ig = image_geom(nx=16, nz=4, dx=2, dz=3)
	image_geom_test3(ig)
	true
end

@testset "2d" begin
	@test all(ImageGeom{1}((2,), (3,), (0,)).mask)
	@test image_geom_test2()
end
@testset "3d" begin
	@test image_geom_test3()
end


text:
import Pkg; Pkg.activate(".")
using Distributed
@time @everywhere include("src/required.jl")

## Get & prepare data
@time begin
    # Define coordinates range
    lon_range = (-145.0, -50.0)
    lat_range = (20.0, 75.0)
end

## Get environmental data data
# WorldClim data
@time wc_vars = pmap(x -> worldclim(x, resolution = "10")[lon_range, lat_range], [1,12]);
# Landcover data
@time lc_vars = landcover(1:10, resolution = "10")[lon_range, lat_range]
# Combine environmental data
env_vars = vcat(wc_vars, lc_vars)

# Plot wcvars1
wc_plot = plotSDM(wc_vars[1], c=:auto)
heatmap!(wc_plot, clim=(-10.0, 30.0), colorbar_title="Annual Mean Temperature (°C)", dpi=300)
# Plot lcvars8 (urban)
lc_plot = plotSDM(lc_vars[2], c=:auto)
heatmap!(lc_plot, colorbar_title="Crops land cover (%)", dpi=300)

## Load data
@load "data/jld2/raw-pres-abs.jld2" pres_abs
@load "data/jld2/sdm-predictions-landcover.jld2" predictions

## Plot Single species
replace!(pres_abs[13].grid, 0.0 => NaN)
singlesp_raw = plotSDM(pres_abs[13], c=:BuPu)
title!(singlesp_raw, "")
heatmap!(singlesp_raw, clim=(0.0, 1.0), colorbar_title="Probability of occurrence", dpi=300)
heatmap!(colorbar=:none)
scatter!(singlesp_raw, [NaN], label="Occurrence", color=:purple, markershape=:rect, markersize=2,
                        legend=:bottomright, legendfontsize=4)
singlesp_sdm = plotSDM(predictions[13], c=:BuPu)
title!(singlesp_sdm, "")
heatmap!(singlesp_sdm, clim=(0.0, 1.0), colorbar_title="Probability of occurrence", dpi=300)

# Combine figures
singlesp_plots = plot(singlesp_raw, singlesp_sdm)
                      #layout=(2,1), size=(300,480))
## Raw
# Richness
@time include("raw/03_raw_richness.jl")
richness_raw = richness_plot
title!(richness_raw, "")
heatmap!(richness_raw, clim=(0.0, 60.0), colorbar_title = "Number of species per site", dpi=300)

# LCBD
@time include("raw/05_raw_lcbd.jl")
lcbd_raw = lcbd_plot2
title!(lcbd_raw, "")
heatmap!(lcbd_raw, colorbar_title = "LCBD value (relative to maximum)", dpi=300)

# Relationship
@time include("raw/06_raw_relation-lcbd-richness.jl")
relation_raw = relationtr_plot
title!(relation_raw, "")
xaxis!(relation_raw, (0.0, 1.0))
yaxis!(relation_raw, (0.0, 1.0))

# Keep richness & lcbd values for combined relationship plot
rel_richness_raw = vec(rel_richness)
rel_lcbd_raw = vec(LCBD[2].grid)

## SDM

# Richness
@time include("sdm/03_sdm_richness.jl")
richness_sdm = richness_plot
title!(richness_sdm, "")
heatmap!(richness_sdm, clim=(0.0, 60.0), colorbar_title = "Number of species per site", dpi=300)

# LCBD
@time include("sdm/05_sdm_lcbd.jl")
lcbd_sdm = lcbd_plot1
title!(lcbd_sdm, "")
heatmap!(lcbd_sdm, colorbar_title = "LCBD value (relative to maximum)", dpi=300)

# Relationship
@time include("sdm/06_sdm_relation-lcbd-richness.jl")
relation_sdm = relation_plot
title!(relation_sdm, "")
xaxis!(relation_sdm, (0.0, 1.0))
yaxis!(relation_sdm, (0.0, 1.0))

# Keep richness & lcbd values for combined relationship plot
rel_richness_sdm = vec(rel_richness)
rel_lcbd_sdm = vec(LCBD[1].grid)

## Combine plots
# Side by side
richness_plots = plot(richness_raw, richness_sdm)
lcbd_plots = plot(lcbd_raw, lcbd_sdm)
relation_plots = plot(relation_raw, relation_sdm)

# Relationships on same plot
relation_oneplot = scatter(rel_richness_raw, rel_lcbd_raw,
        markersize = 1,
        c = :skyblue,
        label = "Raw occurrence data",
        legend = :topright,
        xlims = (0.0, 1.0), ylims = (0.0, 1.0),
        yticks = 0.0:0.20:1.0,
        xlabel = "Species richness (\\alpha\\/\\gamma)", ylabel = "LCBD (relative to maximum)",
        grid=:none)
scatter!(relation_oneplot, rel_richness_sdm, rel_lcbd_sdm,
         markersize = 1, color=:orange, label = "SDM predictions")
# plot!(relation_oneplot, aspectratio=1)

relation_oneplot_empty = scatter(
        markersize = 1,
        c = :skyblue,
        label = "Raw occurrence data",
        legend = :bottomright,
        xlims = (0.0, 1.0), ylims = (0.0, 1.0),
        yticks = 0.0:0.20:1.0,
        xlabel = "Species richness (\\alpha\\/\\gamma)", ylabel = "LCBD (relative to maximum)",
        grid=:none)

# Mean LCBD / richness
rel_df_raw = DataFrame(richness = rel_richness_raw, lcbd = rel_lcbd_raw)
rel_mean_raw = sort(by(rel_df_raw, :richness, lcbd_mean = :lcbd => mean), :richness)
rel_mean_plot = scatter(rel_mean_raw.richness, rel_mean_raw.lcbd_mean)

rel_df_sdm = DataFrame(richness = rel_richness_sdm, lcbd = rel_lcbd_sdm)
rel_mean_sdm = sort(by(rel_df_sdm, :richness, lcbd_mean = :lcbd => mean), :richness)
scatter!(rel_mean_plot, rel_mean_sdm.richness, rel_mean_sdm.lcbd_mean)

# Median LCBD / richness
rel_df_raw = DataFrame(richness = rel_richness_raw, lcbd = rel_lcbd_raw)
rel_median_raw = sort(by(rel_df_raw, :richness, lcbd_median = :lcbd => median), :richness)
rel_median_plot = plot(rel_median_raw.richness, rel_median_raw.lcbd_median, smooth=true)

rel_df_sdm = DataFrame(richness = rel_richness_sdm, lcbd = rel_lcbd_sdm)
rel_median_sdm = sort(by(rel_df_sdm, :richness, lcbd_median = :lcbd => median), :richness)
plot!(rel_median_plot, rel_median_sdm.richness, rel_median_sdm.lcbd_median, smooth=true)

# Density
#=
using StatsPlots
density_richness_raw = density(filter(!isnan, rel_richness_raw))
density_richness_sdm = density(filter(!isnan, rel_richness_sdm))

density_lcbd_raw = density(filter(!isnan, rel_lcbd_raw))
density_lcbd_sdm = density(filter(!isnan, rel_lcbd_sdm))
=#

## Export figures

# env vars
savefig(wc_plot, "doc/fig/wc_temp.png")
savefig(lc_plot, "doc/fig/lc_temp.png")

# Raw
savefig(singlesp_raw, "doc/fig/01_raw_singlesp.png")
savefig(richness_raw, "doc/fig/03_raw_richness.png")
savefig(lcbd_raw, "doc/fig/05_raw_lcbd-transf.png")
savefig(relation_raw, "doc/fig/06_raw_relation.png")

# SDM
savefig(singlesp_sdm, "doc/fig/01_sdm_singlesp.png")
savefig(richness_sdm, "doc/fig/03_sdm_richness.png")
savefig(lcbd_sdm, "doc/fig/05_sdm_lcbd.png")
savefig(relation_sdm, "doc/fig/06_sdm_relation.png")

# Combined
savefig(richness_plots, "doc/fig/03_cmb_richness.png")
savefig(lcbd_plots, "doc/fig/05_cmb_lcbd.png")
savefig(relation_plots, "doc/fig/06_cmb_relation.png")
savefig(relation_oneplot, "doc/fig/06_cmb_relation-oneplot.png")
savefig(relation_oneplot_empty, "doc/fig/06_cmb_relation-oneplot-empty.png")

## Single species SDM with threshold (run if needed, needs to load whole dataset, takes a while to run)
#=
# Option 1: sdm_single-species script, fast but not exactly the same
@time include("sdm/sdm_single-species.jl")
singlesp_sdm_threshold = maps[2]
title!(singlesp_sdm_threshold, "")
heatmap!(singlesp_sdm_threshold, clim=(0.0, 1.0), colorbar_title="Probability of occurrence", dpi=300)
savefig(singlesp_sdm_threshold, "doc/fig/01_sdm_singlesp-threshold.png")

# Option 2: sdm_predictions script, exact same but soooo long since not parallelized
## Exact same prediction, but soooooo long since not parallelized
## Get & prepare data
@time begin
    # Load data from CSV files
    df = CSV.read("data/proc/ebd_warblers_prep.csv", header=true, delim="\t")
    # Separate species
    warblers_occ = [df[df.species .== u,:] for u in unique(df.species)]

    # Define coordinates range
    lon_range = (-145.0, -50.0)
    lat_range = (20.0, 75.0)
    # Observed coordinates range
    lon_range_obs = extrema(df.longitude)
    lat_range_obs = extrema(df.latitude)
end

## Get environmental data
# WorldClim data
@time wc_vars = pmap(x -> worldclim(x, resolution = "10")[lon_range, lat_range], [1,12]);
# WorldClim data with different training resolutions
#=
@time wc_vars_pred = pmap(x -> worldclim(x, resolution = "10")[lon_range, lat_range], 1:19);
@time wc_vars_train = pmap(x -> worldclim(x, resolution = "5")[lon_range_obs, lat_range_obs], 1:19);
=#
# Landcover data
@time lc_vars = load_landcover(lon_range, lat_range)
# Combine environmental data
env_vars = vcat(wc_vars, lc_vars)

## Make prediction for Yellow Warbler
@time singlesp_pred = species_bclim(warblers_occ[13], env_vars, with_threshold=true)
singlesp_sdm_threshold = plotSDM(singlesp_pred, c=:BuPu)
title!(singlesp_sdm_threshold, "")
heatmap!(singlesp_sdm_threshold, clim=(0.0, 1.0), colorbar_title="Probability of occurrence", dpi=300)
savefig(singlesp_sdm_threshold, "doc/fig/01_sdm_singlesp-threshold.png")
=#

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
<!DOCTYPE html>

    <title>Title</title>

Hello!

max_stars_repo_name: piecubed/angel

max_issues_repo_name: piecubed/angel

max_forks_repo_name: piecubed/angel


text:
using Quon
using ZXCalculus

q = Tait(2)
push_gate!(q, Val(:CZ), 1, 2)
push_gate!(q, Val(:CZ), 1, 2)
Quon.simplify!(q)

plot(q)
plot(simplify(q))
plot(q; backend = :compose)
zxg = ZXDiagram(q) |> ZXGraph |> full_reduction
zxg |> ancilla_extraction |> plot

using ZXCalculus
# s = tait_swap()
# zxd = ZXDiagram(s)
# zxg = ZXGraph(zxd)
# full_reduction(zxg)
# circ = ancilla_extraction(zxg)
# circ |> plot

# zxg = ZXGraph(zxd)
# YaoPlots.plot(zxg)
# clifford_simplification(zxg)
# YaoPlots.plot(zxg)
# ZXCalculus.generate_layout!(zxg)

# rz = tait_rz(pi*im)
# zxd_rz = ZXDiagram(rz)
# zxd_h = ZXDiagram(tait_hadamard())
# zxd_cz = ZXDiagram(tait_cz())
zxg_cnot(i, j) = tait_cnot(i, j) |> ZXDiagram |> ZXGraph |> full_reduction |> ancilla_extraction
zxg_cnot(1, 30) |> plot


text:
using PrettyTables

header = [ "Header 1" "Header 2" "Header 3" "Header 4";
              "Sub 1"    "Sub 2"    "Sub 3"    "Sub 4"]
data = [ true  100.0 0x8080 "String";
         false 200.0 0x0808 "String";
         true  300.0 0x1986 "String";
         false 400.0 0x1987 "String"; ]

for (filename, format) in
    (("./html_format_default.html",    html_default),
     ("./html_format_dark.html",       html_dark),
     ("./html_format_minimalist.html", html_minimalist),
     ("./html_format_simple.html",     html_simple))

    open(filename, "w") do f
        pretty_table(f, data, header, backend = Val(:html), tf = format)
    end
end

open("./html_format_matrix.html", "w") do f
    pretty_table(f, data, tf = html_matrix, noheader = true)
end


text:
module Mifid

import ..Clarus

function dsb(;params...)
  return Clarus.api_request("MIFID","DSB";params...)
end

function firds(;params...)
  return Clarus.api_request("MIFID","FIRDS";params...)
end

end #MODULE-END


text:
module Persistence

# exports
export ObjectPool, PersistentOID, transaction

# includes
include("lib.jl")
include("transaction.jl")
include("persist.jl")
include("trace.jl")
include("arrays.jl")

# include submodules
using .Lib
using .Transaction
using .Persist
using .Arrays

end # module


text:
mutable struct Instrument
    #ADSR envelope for amplitude
    envelope
    #a function for generating a combination of sine waves or other type of waveforms
    waveform
end
using WAV

function sine(seconds, hz, fs, vib_amp, vib_rate)
    t = 0.0 : inv(fs) : seconds
    freq_vib = hz .+ vib_amp .* sin.(t * 2π * vib_rate)
    phase_vib = zeros(length(t))
    for i in 2:length(t)
        phase_vib[i] = phase_vib[i-1] + 2π* freq_vib[i-1] / fs
    end
    sin.(phase_vib)
end

function sin_weird(amplitude,hz)
    amplitude*sine(1,hz,41000,log2(hz*(rand()+rand()))/1+2,log2(hz*(rand()+rand()))/2+4)
end
function sin_weird2(amplitude,hz)
    amplitude*sine(1,hz,41000,time()%5+3,time()%3+3)

end
function sine_source(wave_count,n, s)

    freqs = 0:1:(wave_count-1)

    hz ->  sum(s.(1 ./((freqs.+1) .^ n ), hz * (2 .^ freqs)))

end

function norm(a)
    a/maximum(a)
end
function compound(x)
    f1 = sine_source(24, 1.2, sin_weird)
    norm(f1(x))
end

function compound2(x)
    f1 = sine_source(25, 1.2, sin_weird2)
    norm(f1(x))
end
function test()
   wavplay(compound(440), 41000)
end

max_stars_repo_name: justinbroce/julia-audio-stuffs

max_issues_repo_name: justinbroce/julia-audio-stuffs

max_forks_repo_name: justinbroce/julia-audio-stuffs


text:
module FourierAcceleration

using FFTW
using LinearAlgebra

using ..Models: AbstractModel
using ..Utilities: get_index, reshaped

export FourierAccelerator, update_Q!, update_M!, fourier_accelerate!

struct FourierAccelerator{T<:AbstractFloat,Tfft<:AbstractFFTs.Plan,Tifft<:AbstractFFTs.Plan}

    "Copy of input vector if input vector is real."
    vin::Vector{Complex{T}}

    "Copy of output vector if output vector is real."
    vout::Vector{Complex{T}}

    "Vector to store data associated with FFT of single time slice"
    u::Vector{Complex{T}}

    "Vector representing diagonal acceleration matrix using old convention."
    Q::Vector{T}

    "Vector representing c1 matrix using new convetion."
    M::Vector{T}

    "Performs forward fourier transformation"
    pfft::Tfft

    "Performs forward fourier transformation"
    pifft::Tifft

    "Number of Phonons in Lattices."
    N::Int

    "Length of imagniary time axis"
    L::Int

    #######################
    ## INNER CONSTRUCTOR ##
    #######################

    """
    Constructor for FourierAccelerator type.
    """
    function FourierAccelerator(model::AbstractModel{T1,T2,T3}) where {T1,T2,T3}

        # getting number of sites in lattice
        N = model.Nph

        # length of imaginary time axis
        L = model.Lτ

        # constructing Q and M
        Q = zeros(T1,N*L)

        # constructing M
        M = zeros(T1,N*L)

        # declaring a temporary vector for storage purpose if input vectors are real
        vin = zeros(Complex{T1},N*L)

        # declaring a temporary vector for storage purpose if output vectors are real
        vout = zeros(Complex{T1},N*L)

        # declaring temporary storage vector to represent vectors in frequency space.
        u = zeros(Complex{T1},N*L)

        # vector for planning FFT
        v = zeros(Complex{T1},L,N)

        # planning forward FFT
        pfft  = plan_fft(v, (1,), flags=FFTW.PATIENT)

        # planning inverse FFT
        pifft = plan_ifft(v, (1,), flags=FFTW.PATIENT)

        new{T1,typeof(pfft),typeof(pifft)}(vin,vout,u,Q,M,pfft,pifft,N,L)
    end

end

#######################################################
## FUNCITONS ASSOCIATED WITH FourierAccelerator TYPE ##
#######################################################

"""
Accelerate vector by multiplying with Q matrix.
"""
function fourier_accelerate!(v′::AbstractVector{Complex{T}}, fa::FourierAccelerator{T}, v::AbstractVector{Complex{T}}, power::T; use_mass::Bool=false) where {T<:AbstractFloat}

    N = fa.N # number of sites in lattice
    L = fa.L # length of imaginary time axis
    u = fa.u::Vector{Complex{T}}

    a  = reshaped(v,L,N)
    a′ = reshaped(v′,L,N)
    u′ = reshaped(u,L,N)

    # τ → ω
    mul!(u′,fa.pfft,a)

    # apply diagonal acceleration matrix Q^power
    if use_mass
        @. u *= fa.M^power
    else
        @. u *= fa.Q^power
    end

    # ω → τ
    mul!(a′,fa.pifft,u′)

    return nothing
end

function fourier_accelerate!(v′::AbstractVector{Complex{T}}, fa::FourierAccelerator{T}, v::AbstractVector{T}, power::T; use_mass::Bool=false) where {T<:AbstractFloat}

    copyto!(fa.vin,v)
    fourier_accelerate!(v′,fa,fa.vin,power,use_mass=use_mass)
    return nothing
end

function fourier_accelerate!(v′::AbstractVector{T}, fa::FourierAccelerator{T}, v::AbstractVector{Complex{T}}, power::T; use_mass::Bool=false) where {T<:AbstractFloat}

    fourier_accelerate!(fa.vout,fa,v,power,use_mass=use_mass)
    @. v′ = real(fa.vout)
    return nothing
end

function fourier_accelerate!(v′::AbstractVector{T}, fa::FourierAccelerator{T}, v::AbstractVector{T}, power::T; use_mass::Bool=false) where {T<:AbstractFloat}

    copyto!(fa.vin,v)
    fourier_accelerate!(fa.vout,fa,fa.vin,power,use_mass=use_mass)
    @. v′ = real(fa.vout)
    return nothing
end

function fourier_accelerate!(v::AbstractVector, fa::FourierAccelerator{T}, power::T; use_mass::Bool=false) where {T<:AbstractFloat}

    fourier_accelerate!(v,fa,v,power,use_mass=use_mass)
    return nothing
end

"""
Updates the fourier acceleration matrix for sites with phonon frequencies withing the specified range.
"""
function update_Q!(fa::FourierAccelerator{T1},model::AbstractModel{T1,T2},ω_min::T1,ω_max::T1,m::T1) where {T1<:AbstractFloat,T2<:Number}

    # updating the acceleration matrix for sites with a phonon frequency withing the specified range
    update_Q!(fa.Q,model,ω_min,ω_max,m)

    return nothing
end

"""
Updates the c1 matrix for sites with phonon frequencies withing the specified range.
"""
function update_M!(fa::FourierAccelerator{T1},model::AbstractModel{T1,T2},ω_min::T1,ω_max::T1,m0::T1,c::T1=0.0) where {T1<:AbstractFloat,T2<:Number}

    # updating the acceleration matrix for sites with a phonon frequency withing the specified range
    update_M!(fa.M,model,ω_min,ω_max,m0,c)

    return nothing
end

#######################
## PRIVATE FUNCTIONS ##
#######################

"""
Updates the fourier acceleration matrix for sites with phonon frequencies withing the specified range.
"""
function update_Q!(Q::Vector{T1},model::AbstractModel{T1,T2},ω_min::T1,ω_max::T1,m::T1) where {T1<:AbstractFloat,T2<:Number}

    N  = model.Nph::Int
    L  = model.Lτ
    Δτ = model.Δτ::T1
    ω  = model.ω
    # iterating over phonons in lattice
    for phonon in 1:N
        # if phonon frequncy falls withing specified range
        if ω_min < ω[phonon] < ω_max
            # get a view into Q matrix for current lattice phonon
            Qi = @view Q[get_index(1,phonon,L):get_index(L,phonon,L)]
            # define Q matrix just for current phonon
            construct_Qi!( Qi , ω[phonon] , Δτ , m )
        end
    end
    return nothing
end

"""
Calculates acceleration matrix for specified phonon frequency `ω`, discretization `Δτ`.
Obeys the FFTW convention for the ordering of the momentum values.
"""
function construct_Qi!(Qi::AbstractVector{T},ω::T,Δτ::T,m::T) where {T<:AbstractFloat}

    L = length(Qi)
    for k in 0:L-1
        Qi[k+1] = element_Qi(k,ω,Δτ,m,L)
    end
    return nothing
end

"""
Calculates a specified matrix element of the acceleration matrix for a given mode k.
"""
function element_Qi(k::Int,ω::T,Δτ::T,m::T,L::Int)::T where {T<:Number}

    val = (m^2 + Δτ*ω*ω + 4.0/Δτ) / (m^2 + Δτ*ω*ω + (2-2*cos(2*π*k/L))/Δτ)
    return val
end

"""
Updates the c1 matrix for sites with phonon frequencies withing the specified range.
"""
function update_M!(M::Vector{T1},model::AbstractModel{T1,T2},ω_min::T1,ω_max::T1,m0::T1,c::T1) where {T1<:AbstractFloat,T2<:Number}

    N  = model.Nph::Int
    L  = model.Lτ
    Δτ = model.Δτ::T1
    ω  = model.ω
    # iterating over phonons in lattice
    for phonon in 1:N
        # if phonon frequncy on site falls withing specified range
        if ω_min < ω[phonon] < ω_max
            # get a view into M matrix for current phonon
            Mi = @view M[get_index(1,phonon,L):get_index(L,phonon,L)]
            # define M matrix just for current phonon
            construct_Mi!( Mi , ω[phonon] , Δτ , m0 , c )
        end
    end
    return nothing
end

"""
Calculates matrix matrix for specified phonon frequency `ω`, discretization `Δτ`.
Obeys the FFTW convention for the ordering of the momentum values.
"""
function construct_Mi!(Mi::AbstractVector{T},ω::T,Δτ::T,m0::T,c::T) where {T<:AbstractFloat}

    L = length(Mi)
    for k in 0:L-1
        Mi[k+1] = element_Mi(k,ω,Δτ,m0,c,L)
    end
    return nothing
end

"""
Calculates a specified matrix element of the c1 matrix for a given mode k.
"""
function element_Mi(k::Int,ω::T,Δτ::T,m0::T,c::T,L::Int)::T where {T<:Number}

    k′  = min(k,L-k)
    m   = m0 * exp(-(c*k′/L)^2)
    val = Δτ * (m^2 + ω^2 + (2-2*cos(2*π*k′/L))/Δτ^2 ) / (m^2 + ω^2)
    return val
end

end

max_stars_repo_name: cohensbw/Langevin

max_issues_repo_name: cohensbw/ElPhDynamics

max_forks_repo_name: cohensbw/ElPhDynamics


text:
"""
    heideldiag(
        x::AbstractVector{<:Real}; alpha::Real=0.05, eps::Real=0.1, start::Int=1, kwargs...
    )

Compute the Heidelberger and Welch diagnostic [^Heidelberger1983]. This diagnostic tests for
non-convergence (non-stationarity) and whether ratios of estimation interval halfwidths to
means are within a target ratio. Stationarity is rejected (0) for significant test p-values.
Halfwidth tests are rejected (0) if observed ratios are greater than the target, as is the
case for `s2` and `beta[1]`.

[^Heidelberger1983]: Heidelberger, P., & Welch, P. D. (1983). Simulation run length control in the presence of an initial transient. Operations Research, 31(6), 1109-1144.
"""
function heideldiag(
    x::AbstractVector{<:Real}; alpha::Real=0.05, eps::Real=0.1, start::Int=1, kwargs...
)
    n = length(x)
    delta = trunc(Int, 0.10 * n)
    y = x[trunc(Int, n / 2):end]
    S0 = length(y) * mcse(y; kwargs...)^2
    i, pvalue, converged, ybar = 1, 1.0, false, NaN
    while i < n / 2
        y = x[i:end]
        m = length(y)
        ybar = Statistics.mean(y)
        B = cumsum(y) - ybar * collect(1:m)
        Bsq = (B .* B) ./ (m * S0)
        I = sum(Bsq) / m
        pvalue = 1.0 - pcramer(I)
        converged = pvalue > alpha
        if converged
            break
        end
        i += delta
    end
    halfwidth = sqrt(2) * SpecialFunctions.erfcinv(alpha) * mcse(y; kwargs...)
    passed = halfwidth / abs(ybar) <= eps
    return (
        burnin=i + start - 2,
        stationarity=converged,
        pvalue=pvalue,
        mean=ybar,
        halfwidth=halfwidth,
        test=passed,
    )
end

## Csorgo S and Faraway JJ. The exact and asymptotic distributions of the
## Cramer-von Mises statistic. Journal of the Royal Statistical Society,
## Series B, 58: 221-234, 1996.
function pcramer(q::Real)
    p = 0.0
    for k in 0:3
        c1 = 4.0 * k + 1.0
        c2 = c1^2 / (16.0 * q)
        p +=
            SpecialFunctions.gamma(k + 0.5) / factorial(k) *
            sqrt(c1) *
            exp(-c2) *
            SpecialFunctions.besselk(0.25, c2)
    end
    return p / (pi^1.5 * sqrt(q))
end


text:
# syntax: proto2
using Compat
using ProtoBuf
import ProtoBuf.meta
import ..hadoop

struct __enum_SchedulerResourceTypes <: ProtoEnum
    MEMORY::Int32
    CPU::Int32
    __enum_SchedulerResourceTypes() = new(0,1)
end #struct __enum_SchedulerResourceTypes
const SchedulerResourceTypes = __enum_SchedulerResourceTypes()

struct __enum_ApplicationsRequestScopeProto <: ProtoEnum
    ALL::Int32
    VIEWABLE::Int32
    OWN::Int32
    __enum_ApplicationsRequestScopeProto() = new(0,1,2)
end #struct __enum_ApplicationsRequestScopeProto
const ApplicationsRequestScopeProto = __enum_ApplicationsRequestScopeProto()

mutable struct RegisterApplicationMasterRequestProto <: ProtoType
    host::AbstractString
    rpc_port::Int32
    tracking_url::AbstractString
    RegisterApplicationMasterRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct RegisterApplicationMasterRequestProto

mutable struct FinishApplicationMasterRequestProto <: ProtoType
    diagnostics::AbstractString
    tracking_url::AbstractString
    final_application_status::Int32
    FinishApplicationMasterRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct FinishApplicationMasterRequestProto

mutable struct FinishApplicationMasterResponseProto <: ProtoType
    isUnregistered::Bool
    FinishApplicationMasterResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct FinishApplicationMasterResponseProto
const __val_FinishApplicationMasterResponseProto = Dict(:isUnregistered => false)
meta(t::Type{FinishApplicationMasterResponseProto}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, __val_FinishApplicationMasterResponseProto, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct AllocateRequestProto <: ProtoType
    ask::Base.Vector{ResourceRequestProto}
    release::Base.Vector{ContainerIdProto}
    blacklist_request::ResourceBlacklistRequestProto
    response_id::Int32
    progress::Float32
    increase_request::Base.Vector{ContainerResourceIncreaseRequestProto}
    AllocateRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct AllocateRequestProto

mutable struct NMTokenProto <: ProtoType
    nodeId::NodeIdProto
    token::hadoop.common.TokenProto
    NMTokenProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct NMTokenProto

mutable struct RegisterApplicationMasterResponseProto <: ProtoType
    maximumCapability::ResourceProto
    client_to_am_token_master_key::Array{UInt8,1}
    application_ACLs::Base.Vector{ApplicationACLMapProto}
    containers_from_previous_attempts::Base.Vector{ContainerProto}
    queue::AbstractString
    nm_tokens_from_previous_attempts::Base.Vector{NMTokenProto}
    scheduler_resource_types::Base.Vector{Int32}
    RegisterApplicationMasterResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct RegisterApplicationMasterResponseProto

mutable struct AllocateResponseProto <: ProtoType
    a_m_command::Int32
    response_id::Int32
    allocated_containers::Base.Vector{ContainerProto}
    completed_container_statuses::Base.Vector{ContainerStatusProto}
    limit::ResourceProto
    updated_nodes::Base.Vector{NodeReportProto}
    num_cluster_nodes::Int32
    preempt::PreemptionMessageProto
    nm_tokens::Base.Vector{NMTokenProto}
    increased_containers::Base.Vector{ContainerResourceIncreaseProto}
    decreased_containers::Base.Vector{ContainerResourceDecreaseProto}
    am_rm_token::hadoop.common.TokenProto
    AllocateResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct AllocateResponseProto

mutable struct GetNewApplicationRequestProto <: ProtoType
    GetNewApplicationRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetNewApplicationRequestProto

mutable struct GetNewApplicationResponseProto <: ProtoType
    application_id::ApplicationIdProto
    maximumCapability::ResourceProto
    GetNewApplicationResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetNewApplicationResponseProto

mutable struct GetApplicationReportRequestProto <: ProtoType
    application_id::ApplicationIdProto
    GetApplicationReportRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetApplicationReportRequestProto

mutable struct GetApplicationReportResponseProto <: ProtoType
    application_report::ApplicationReportProto
    GetApplicationReportResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetApplicationReportResponseProto

mutable struct SubmitApplicationRequestProto <: ProtoType
    application_submission_context::ApplicationSubmissionContextProto
    SubmitApplicationRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct SubmitApplicationRequestProto

mutable struct SubmitApplicationResponseProto <: ProtoType
    SubmitApplicationResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct SubmitApplicationResponseProto

mutable struct KillApplicationRequestProto <: ProtoType
    application_id::ApplicationIdProto
    KillApplicationRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct KillApplicationRequestProto

mutable struct KillApplicationResponseProto <: ProtoType
    is_kill_completed::Bool
    KillApplicationResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct KillApplicationResponseProto
const __val_KillApplicationResponseProto = Dict(:is_kill_completed => false)
meta(t::Type{KillApplicationResponseProto}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, __val_KillApplicationResponseProto, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct GetClusterMetricsRequestProto <: ProtoType
    GetClusterMetricsRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetClusterMetricsRequestProto

mutable struct GetClusterMetricsResponseProto <: ProtoType
    cluster_metrics::YarnClusterMetricsProto
    GetClusterMetricsResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetClusterMetricsResponseProto

mutable struct MoveApplicationAcrossQueuesRequestProto <: ProtoType
    application_id::ApplicationIdProto
    target_queue::AbstractString
    MoveApplicationAcrossQueuesRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct MoveApplicationAcrossQueuesRequestProto
const __req_MoveApplicationAcrossQueuesRequestProto = Symbol[:application_id,:target_queue]
meta(t::Type{MoveApplicationAcrossQueuesRequestProto}) = meta(t, __req_MoveApplicationAcrossQueuesRequestProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct MoveApplicationAcrossQueuesResponseProto <: ProtoType
    MoveApplicationAcrossQueuesResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct MoveApplicationAcrossQueuesResponseProto

mutable struct GetApplicationsRequestProto <: ProtoType
    application_types::Base.Vector{AbstractString}
    application_states::Base.Vector{Int32}
    users::Base.Vector{AbstractString}
    queues::Base.Vector{AbstractString}
    limit::Int64
    start_begin::Int64
    start_end::Int64
    finish_begin::Int64
    finish_end::Int64
    applicationTags::Base.Vector{AbstractString}
    scope::Int32
    GetApplicationsRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetApplicationsRequestProto
const __val_GetApplicationsRequestProto = Dict(:scope => ApplicationsRequestScopeProto.ALL)
meta(t::Type{GetApplicationsRequestProto}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, __val_GetApplicationsRequestProto, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct GetApplicationsResponseProto <: ProtoType
    applications::Base.Vector{ApplicationReportProto}
    GetApplicationsResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetApplicationsResponseProto

mutable struct GetClusterNodesRequestProto <: ProtoType
    nodeStates::Base.Vector{Int32}
    GetClusterNodesRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetClusterNodesRequestProto

mutable struct GetClusterNodesResponseProto <: ProtoType
    nodeReports::Base.Vector{NodeReportProto}
    GetClusterNodesResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetClusterNodesResponseProto

mutable struct GetQueueInfoRequestProto <: ProtoType
    queueName::AbstractString
    includeApplications::Bool
    includeChildQueues::Bool
    recursive::Bool
    GetQueueInfoRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetQueueInfoRequestProto

mutable struct GetQueueInfoResponseProto <: ProtoType
    queueInfo::QueueInfoProto
    GetQueueInfoResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetQueueInfoResponseProto

mutable struct GetQueueUserAclsInfoRequestProto <: ProtoType
    GetQueueUserAclsInfoRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetQueueUserAclsInfoRequestProto

mutable struct GetQueueUserAclsInfoResponseProto <: ProtoType
    queueUserAcls::Base.Vector{QueueUserACLInfoProto}
    GetQueueUserAclsInfoResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetQueueUserAclsInfoResponseProto

mutable struct GetNodesToLabelsRequestProto <: ProtoType
    GetNodesToLabelsRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetNodesToLabelsRequestProto

mutable struct GetNodesToLabelsResponseProto <: ProtoType
    nodeToLabels::Base.Vector{NodeIdToLabelsProto}
    GetNodesToLabelsResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetNodesToLabelsResponseProto

mutable struct GetLabelsToNodesRequestProto <: ProtoType
    nodeLabels::Base.Vector{AbstractString}
    GetLabelsToNodesRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetLabelsToNodesRequestProto

mutable struct GetLabelsToNodesResponseProto <: ProtoType
    labelsToNodes::Base.Vector{LabelsToNodeIdsProto}
    GetLabelsToNodesResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetLabelsToNodesResponseProto

mutable struct GetClusterNodeLabelsRequestProto <: ProtoType
    GetClusterNodeLabelsRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetClusterNodeLabelsRequestProto

mutable struct GetClusterNodeLabelsResponseProto <: ProtoType
    nodeLabels::Base.Vector{AbstractString}
    GetClusterNodeLabelsResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetClusterNodeLabelsResponseProto

mutable struct StartContainerRequestProto <: ProtoType
    container_launch_context::ContainerLaunchContextProto
    container_token::hadoop.common.TokenProto
    StartContainerRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct StartContainerRequestProto

mutable struct StartContainerResponseProto <: ProtoType
    services_meta_data::Base.Vector{StringBytesMapProto}
    StartContainerResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct StartContainerResponseProto

mutable struct StopContainerRequestProto <: ProtoType
    container_id::ContainerIdProto
    StopContainerRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct StopContainerRequestProto

mutable struct StopContainerResponseProto <: ProtoType
    StopContainerResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct StopContainerResponseProto

mutable struct GetContainerStatusRequestProto <: ProtoType
    container_id::ContainerIdProto
    GetContainerStatusRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetContainerStatusRequestProto

mutable struct GetContainerStatusResponseProto <: ProtoType
    status::ContainerStatusProto
    GetContainerStatusResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetContainerStatusResponseProto

mutable struct StartContainersRequestProto <: ProtoType
    start_container_request::Base.Vector{StartContainerRequestProto}
    StartContainersRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct StartContainersRequestProto

mutable struct ContainerExceptionMapProto <: ProtoType
    container_id::ContainerIdProto
    exception::SerializedExceptionProto
    ContainerExceptionMapProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ContainerExceptionMapProto

mutable struct StartContainersResponseProto <: ProtoType
    services_meta_data::Base.Vector{StringBytesMapProto}
    succeeded_requests::Base.Vector{ContainerIdProto}
    failed_requests::Base.Vector{ContainerExceptionMapProto}
    StartContainersResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct StartContainersResponseProto

mutable struct StopContainersRequestProto <: ProtoType
    container_id::Base.Vector{ContainerIdProto}
    StopContainersRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct StopContainersRequestProto

mutable struct StopContainersResponseProto <: ProtoType
    succeeded_requests::Base.Vector{ContainerIdProto}
    failed_requests::Base.Vector{ContainerExceptionMapProto}
    StopContainersResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct StopContainersResponseProto

mutable struct GetContainerStatusesRequestProto <: ProtoType
    container_id::Base.Vector{ContainerIdProto}
    GetContainerStatusesRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetContainerStatusesRequestProto

mutable struct GetContainerStatusesResponseProto <: ProtoType
    status::Base.Vector{ContainerStatusProto}
    failed_requests::Base.Vector{ContainerExceptionMapProto}
    GetContainerStatusesResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetContainerStatusesResponseProto

mutable struct GetApplicationAttemptReportRequestProto <: ProtoType
    application_attempt_id::ApplicationAttemptIdProto
    GetApplicationAttemptReportRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetApplicationAttemptReportRequestProto

mutable struct GetApplicationAttemptReportResponseProto <: ProtoType
    application_attempt_report::ApplicationAttemptReportProto
    GetApplicationAttemptReportResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetApplicationAttemptReportResponseProto

mutable struct GetApplicationAttemptsRequestProto <: ProtoType
    application_id::ApplicationIdProto
    GetApplicationAttemptsRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetApplicationAttemptsRequestProto

mutable struct GetApplicationAttemptsResponseProto <: ProtoType
    application_attempts::Base.Vector{ApplicationAttemptReportProto}
    GetApplicationAttemptsResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetApplicationAttemptsResponseProto

mutable struct GetContainerReportRequestProto <: ProtoType
    container_id::ContainerIdProto
    GetContainerReportRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetContainerReportRequestProto

mutable struct GetContainerReportResponseProto <: ProtoType
    container_report::ContainerReportProto
    GetContainerReportResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetContainerReportResponseProto

mutable struct GetContainersRequestProto <: ProtoType
    application_attempt_id::ApplicationAttemptIdProto
    GetContainersRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetContainersRequestProto

mutable struct GetContainersResponseProto <: ProtoType
    containers::Base.Vector{ContainerReportProto}
    GetContainersResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GetContainersResponseProto

mutable struct UseSharedCacheResourceRequestProto <: ProtoType
    applicationId::ApplicationIdProto
    resourceKey::AbstractString
    UseSharedCacheResourceRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct UseSharedCacheResourceRequestProto

mutable struct UseSharedCacheResourceResponseProto <: ProtoType
    path::AbstractString
    UseSharedCacheResourceResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct UseSharedCacheResourceResponseProto

mutable struct ReleaseSharedCacheResourceRequestProto <: ProtoType
    applicationId::ApplicationIdProto
    resourceKey::AbstractString
    ReleaseSharedCacheResourceRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ReleaseSharedCacheResourceRequestProto

mutable struct ReleaseSharedCacheResourceResponseProto <: ProtoType
    ReleaseSharedCacheResourceResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ReleaseSharedCacheResourceResponseProto

mutable struct ReservationSubmissionRequestProto <: ProtoType
    queue::AbstractString
    reservation_definition::ReservationDefinitionProto
    ReservationSubmissionRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ReservationSubmissionRequestProto

mutable struct ReservationSubmissionResponseProto <: ProtoType
    reservation_id::ReservationIdProto
    ReservationSubmissionResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ReservationSubmissionResponseProto

mutable struct ReservationUpdateRequestProto <: ProtoType
    reservation_definition::ReservationDefinitionProto
    reservation_id::ReservationIdProto
    ReservationUpdateRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ReservationUpdateRequestProto

mutable struct ReservationUpdateResponseProto <: ProtoType
    ReservationUpdateResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ReservationUpdateResponseProto

mutable struct ReservationDeleteRequestProto <: ProtoType
    reservation_id::ReservationIdProto
    ReservationDeleteRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ReservationDeleteRequestProto

mutable struct ReservationDeleteResponseProto <: ProtoType
    ReservationDeleteResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ReservationDeleteResponseProto

mutable struct RunSharedCacheCleanerTaskRequestProto <: ProtoType
    RunSharedCacheCleanerTaskRequestProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct RunSharedCacheCleanerTaskRequestProto

mutable struct RunSharedCacheCleanerTaskResponseProto <: ProtoType
    accepted::Bool
    RunSharedCacheCleanerTaskResponseProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct RunSharedCacheCleanerTaskResponseProto

export SchedulerResourceTypes, ApplicationsRequestScopeProto, RegisterApplicationMasterRequestProto, RegisterApplicationMasterResponseProto, FinishApplicationMasterRequestProto, FinishApplicationMasterResponseProto, AllocateRequestProto, NMTokenProto, AllocateResponseProto, GetNewApplicationRequestProto, GetNewApplicationResponseProto, GetApplicationReportRequestProto, GetApplicationReportResponseProto, SubmitApplicationRequestProto, SubmitApplicationResponseProto, KillApplicationRequestProto, KillApplicationResponseProto, GetClusterMetricsRequestProto, GetClusterMetricsResponseProto, MoveApplicationAcrossQueuesRequestProto, MoveApplicationAcrossQueuesResponseProto, GetApplicationsRequestProto, GetApplicationsResponseProto, GetClusterNodesRequestProto, GetClusterNodesResponseProto, GetQueueInfoRequestProto, GetQueueInfoResponseProto, GetQueueUserAclsInfoRequestProto, GetQueueUserAclsInfoResponseProto, GetNodesToLabelsRequestProto, GetNodesToLabelsResponseProto, GetLabelsToNodesRequestProto, GetLabelsToNodesResponseProto, GetClusterNodeLabelsRequestProto, GetClusterNodeLabelsResponseProto, StartContainerRequestProto, StartContainerResponseProto, StopContainerRequestProto, StopContainerResponseProto, GetContainerStatusRequestProto, GetContainerStatusResponseProto, StartContainersRequestProto, ContainerExceptionMapProto, StartContainersResponseProto, StopContainersRequestProto, StopContainersResponseProto, GetContainerStatusesRequestProto, GetContainerStatusesResponseProto, GetApplicationAttemptReportRequestProto, GetApplicationAttemptReportResponseProto, GetApplicationAttemptsRequestProto, GetApplicationAttemptsResponseProto, GetContainerReportRequestProto, GetContainerReportResponseProto, GetContainersRequestProto, GetContainersResponseProto, UseSharedCacheResourceRequestProto, UseSharedCacheResourceResponseProto, ReleaseSharedCacheResourceRequestProto, ReleaseSharedCacheResourceResponseProto, ReservationSubmissionRequestProto, ReservationSubmissionResponseProto, ReservationUpdateRequestProto, ReservationUpdateResponseProto, ReservationDeleteRequestProto, ReservationDeleteResponseProto, RunSharedCacheCleanerTaskRequestProto, RunSharedCacheCleanerTaskResponseProto


text:
using SuperTB, PyPlot

plt.close("all")

tb = TBModel(hopHoneycomb)
tb.pos = buildHoneycombRestriction(x->circularRestriction(x,10), 10)
tb.periodicity = falses(2)

plot(tb, plothopping=true)

# crude dos
ε = solve(tb)
dos(eeig,espan) = map(e -> sum(imag(inv.((e - .01im) .- eeig))), espan)
εˢ = range(minimum(ε), maximum(ε), length=250)

figure()
plot(εˢ, dos(ε, εˢ), "-k")

max_issues_repo_name: thchr/super-tb

max_forks_repo_name: thchr/super-tb


text:
module TestDatasets

# using Revise
using Test
using MLJBase

X, y = @load_boston
@test schema(X).names  == (:Crim, :Zn, :Indus, :NOx, :Rm, :Age,
                           :Dis, :Rad, :Tax, :PTRatio, :Black, :LStat)
@test scitype(y) <: AbstractVector{Continuous}

X, y = @load_ames
schema(X).names == (:MSSubClass, :MSZoning, :LotFrontage, :LotArea,
                    :Street, :LotShape, :LandContour, :LotConfig,
                    :LandSlope, :Neighborhood, :Condition1,
                    :Condition2, :BldgType, :HouseStyle,
                    :OverallQual, :OverallCond, :YearBuilt,
                    :YearRemodAdd, :RoofStyle, :RoofMatl,
                    :Exterior1st, :Exterior2nd, :MasVnrType,
                    :MasVnrArea, :ExterQual, :ExterCond,
                    :Foundation, :BsmtQual, :BsmtCond, :BsmtExposure,
                    :BsmtFinType1, :BsmtFinSF1, :BsmtFinType2,
                    :BsmtFinSF2, :BsmtUnfSF, :TotalBsmtSF, :Heating,
                    :HeatingQC, :CentralAir, :Electrical, :x1stFlrSF,
                    :x2ndFlrSF, :LowQualFinSF, :GrLivArea, :BsmtFullBath,
                    :BsmtHalfBath, :FullBath, :HalfBath,
                    :BedroomAbvGr, :KitchenAbvGr, :KitchenQual,
                    :TotRmsAbvGrd, :Functional, :Fireplaces,
                    :FireplaceQu, :GarageType, :GarageYrBlt,
                    :GarageFinish, :GarageCars, :GarageArea,
                    :GarageQual, :GarageCond, :PavedDrive,
                    :WoodDeckSF, :OpenPorchSF, :EnclosedPorch,
                    :x3SsnPorch, :ScreenPorch, :PoolArea, :MiscVal,
                    :MoSold, :YrSold, :SaleType, :SaleCondition)
@test scitype(y) <: AbstractVector{Continuous}

X, y = @load_reduced_ames
schema(X).names == (:OverallQual, :GrLivArea, :Neighborhood,
                    :x1stFlrSF, :TotalBsmtSF, :BsmtFinSF1, :LotArea,
                    :GarageCars, :MSSubClass, :GarageArea, :YearRemodAdd,
                    :YearBuilt)
@test scitype(y) <: AbstractVector{Continuous}

X, y = @load_iris
@test schema(X).names == (:sepal_length, :sepal_width, :petal_length,
                          :petal_width)
@test scitype(y) <: AbstractVector{<:Multiclass}

X, y = @load_crabs
@test schema(X).names ==  (:FL, :RW, :CL, :CW, :BD)
@test scitype(y) <: AbstractVector{<:Multiclass}

X, y = @load_smarket
@test schema(X).names == (:Year, :Lag1, :Lag2, :Lag3, :Lag4, :Lag5, :Volume, :Today)
@test scitype(y) == AbstractVector{Multiclass{2}}

X = @load_sunspots
@test schema(X).names == (:sunspot_number, )

end # module
true


text:
export zero_problem

struct ZeroFunc{U,P,F}
    f!::F
end

_lift(T, val) = val
_lift(::Type{<:Number}, val) = val[1]

(af::ZeroFunc{U,P})(res, (u, p)) where {U,P} = af.f!(res, _lift(U, u), _lift(P, p))

"""
Construct a zero problem of the form

```math
    0 = f(u, p),
```

where `u` is the state and `p` is the parameter(s). The function can operate on scalars or
vectors, and be in-place or not. It assumes that the function output is of the same
dimension as `u`.

# Parameters

* `name::String` : the name of the zero problem.
* `f` : the function to use for the zero problem. It takes either two arguments
  (`u` and `p`) or three arguments for an in-place version (`res`, `u`, and `p`).
* `u0` : the initial state value (either scalar- or vector-like).
* `p0` : the initial parameter value (either scalar- or vector-like).
* `pnames` : (keyword, optional) the names of the parameters. If not specified,
  auto-generated names will be used.

# Example

```
prob = zero_problem("cubic", (u, p) -> u^3 - p, 1.5, 1)  # u0 = 1.5, p0 = 1
```
"""
function zero_problem(name::String, f, u0, p0; pnames=nothing)
    # Determine whether f is in-place or not
    if any(method.nargs == 4 for method in methods(f))
        f! = f
    else
        f! = (res, u, p) -> res .= f(u, p)
    end
    # Generate parameter names
    _pnames = parameter_names(p0, pnames)
    # Give the user-provided function the input expected
    U = u0 isa Number ? Number : Vector
    P = p0 isa Number ? Number : Vector
    alg = ZeroFunc{U,P,typeof(f!)}(f!)
    # Create the necessary continuation variables and add the function
    u = Var("u"; initial_u=(U === Number ? [u0] : u0))
    p = Var("p"; initial_u=(P === Number ? [p0] : p0))
    func = Func("f", alg; initial_dim=length(u0))
    push!(func, u)
    push!(func, p)
    # TODO: Should also add the parameters as monitor functions
    return push!(Problem(name), func)
end


text:
using Oceananigans.Fields: AbstractField
using Oceananigans.Grids: znode
using Oceananigans.Operators: Δzᵃᵃᶠ, Δzᵃᵃᶜ

const c = Center()
const f = Face()

""" Return the geopotential height at `i, j, k` at cell centers. """
@inline Zᵃᵃᶜ(i, j, k, grid) =
    ifelse(k < 1,       znode(c, c, c, i, j,       1, grid) + (1 - k) * Δzᵃᵃᶠ(i, j, 1, grid),
    ifelse(k > grid.Nz, znode(c, c, c, i, j, grid.Nz, grid) - (k - grid.Nz) * Δzᵃᵃᶠ(i, j, grid.Nz, grid),
                        znode(c, c, c, i, j,       k, grid)))

""" Return the geopotential height at `i, j, k` at cell z-interfaces. """
@inline Zᵃᵃᶠ(i, j, k, grid) =
    ifelse(k < 1,           znode(c, c, f, i, j,           1, grid) + (1 - k) * Δzᵃᵃᶜ(i, j, 1, grid),
    ifelse(k > grid.Nz + 1, znode(c, c, f, i, j, grid.Nz + 1, grid) - (k - grid.Nz + 1) * Δzᵃᵃᶜ(i, j, k, grid),
                            znode(c, c, f, i, j,           k, grid)))

# Dispatch shenanigans
@inline θ_and_sᴬ(i, j, k, θ::AbstractArray, sᴬ::AbstractArray) = @inbounds θ[i, j, k], sᴬ[i, j, k]
@inline θ_and_sᴬ(i, j, k, θ::Number,        sᴬ::AbstractArray) = @inbounds θ, sᴬ[i, j, k]
@inline θ_and_sᴬ(i, j, k, θ::AbstractArray, sᴬ::Number)        = @inbounds θ[i, j, k], sᴬ
@inline θ_and_sᴬ(i, j, k, θ::Number,        sᴬ::Number)        = @inbounds θ, sᴬ

# Basic functionality
@inline ρ′(i, j, k, grid, eos, θ, sᴬ) = @inbounds ρ′(θ_and_sᴬ(i, j, k, θ, sᴬ)..., Zᵃᵃᶜ(i, j, k, grid), eos)

@inline thermal_expansionᶜᶜᶜ(i, j, k, grid, eos, θ, sᴬ) = thermal_expansion(θ_and_sᴬ(i, j, k, θ, sᴬ)..., Zᵃᵃᶜ(i, j, k, grid), eos)
@inline thermal_expansionᶠᶜᶜ(i, j, k, grid, eos, θ, sᴬ) = thermal_expansion(ℑxᶠᵃᵃ(i, j, k, grid, θ), ℑxᶠᵃᵃ(i, j, k, grid, sᴬ), Zᵃᵃᶜ(i, j, k, grid), eos)
@inline thermal_expansionᶜᶠᶜ(i, j, k, grid, eos, θ, sᴬ) = thermal_expansion(ℑyᵃᶠᵃ(i, j, k, grid, θ), ℑyᵃᶠᵃ(i, j, k, grid, sᴬ), Zᵃᵃᶜ(i, j, k, grid), eos)
@inline thermal_expansionᶜᶜᶠ(i, j, k, grid, eos, θ, sᴬ) = thermal_expansion(ℑzᵃᵃᶠ(i, j, k, grid, θ), ℑzᵃᵃᶠ(i, j, k, grid, sᴬ), Zᵃᵃᶠ(i, j, k, grid), eos)

@inline haline_contractionᶜᶜᶜ(i, j, k, grid, eos, θ, sᴬ) = haline_contraction(θ_and_sᴬ(i, j, k, θ, sᴬ)..., Zᵃᵃᶜ(i, j, k, grid), eos)
@inline haline_contractionᶠᶜᶜ(i, j, k, grid, eos, θ, sᴬ) = haline_contraction(ℑxᶠᵃᵃ(i, j, k, grid, θ), ℑxᶠᵃᵃ(i, j, k, grid, sᴬ), Zᵃᵃᶜ(i, j, k, grid), eos)
@inline haline_contractionᶜᶠᶜ(i, j, k, grid, eos, θ, sᴬ) = haline_contraction(ℑyᵃᶠᵃ(i, j, k, grid, θ), ℑyᵃᶠᵃ(i, j, k, grid, sᴬ), Zᵃᵃᶜ(i, j, k, grid), eos)
@inline haline_contractionᶜᶜᶠ(i, j, k, grid, eos, θ, sᴬ) = haline_contraction(ℑzᵃᵃᶠ(i, j, k, grid, θ), ℑzᵃᵃᶠ(i, j, k, grid, sᴬ), Zᵃᵃᶠ(i, j, k, grid), eos)


text:
export Midpoint

"""
    Midpoint()

The (explicit) [Midpoint](https://en.wikipedia.org/wiki/List_of_Runge–Kutta_methods#Explicit_midpoint_method)
predictor of order 2.
"""
struct Midpoint <: AbstractPredictor end
struct MidpointCache{T} <: AbstractPredictorCache
    dt::Vector{T}
    mk₂::Vector{T}
end

function cache(::Midpoint, H, x, ẋ, t)
    MidpointCache(promote(dt(H, x, t), copy(ẋ))...)
end
#
function predict!(
    xnext,
    cache::MidpointCache,
    H::HomotopyWithCache,
    x,
    t,
    Δt,
    ẋ,
    Jac::JacobianMonitor,
)
    dt, mk₂ = cache.dt, cache.mk₂
    n = length(xnext)
    @inbounds for i = 1:n
        xnext[i] = x[i] + 0.5Δt * ẋ[i]
    end

    minus_ẋ!(mk₂, H, xnext, t + 0.5Δt, Jac, dt)

    @inbounds for i = 1:n
        xnext[i] = x[i] - Δt * mk₂[i]
    end

    nothing
end

order(::MidpointCache) = 3


text:
module TestSample
    using StreamStats
    using Base.Test
    using StatsBase

    n = 10
    k = 3
    xs = rand(1:100_000, n)
    @assert length(unique(xs)) == n
    counts = Dict(zip(unique(xs), fill(0, n)))
    counts_state = Dict(zip(unique(xs), fill(0, n)))

    for itr in 1:100_000
        stat = StreamStats.Sample(Int, k)
        for x in xs
            update!(stat, x)
        end
        for sample in sample(stat)
            @test sample in xs
            counts[sample] += 1
        end
        for sample in state(stat)
            @test sample in xs
            counts_state[sample] += 1
        end
    end
    for (elt, count) in counts
        # TODO: Write a sharper test based on SEM's
        @test 29_500 <= count <= 30_500
    end
    for (elt, count) in counts_state
        # TODO: Write a sharper test based on SEM's
        @test 29_500 <= count <= 30_500
    end
end


text:
"""
    In order to send images to Python, we have to convert them to a primtive
    array. We overload view_image here for specific types which would be difficult
    to display in Napari otherwise
"""
module Image

# The main reason to have this in a discrete module is really to contain the imports
# until we can be more explicit about what we are importing
#
# In lieu of PR 876 in PyCall, we use NumPyArrays, to perform no-copy conversions into Python friendly arrays
# https://github.com/JuliaPy/PyCall.jl/pull/876
using Images, FixedPointNumbers, PyCall, AxisArrays, ImageMetadata, NumPyArrays

import Napari: napari, view_image, add_image

# This module

view_image(A::AbstractArray{C}, args...; kwargs...) where C <: Colorant{T,3} where T <: FixedPoint{F} where F =
    view_image(primitive_array(A), args...; kwargs...)
add_image(viewer, A::AbstractArray{C}, args...; kwargs...) where C <: Colorant{T,3} where T <: FixedPoint{F} where F =
    add_image(viewer, primitive_array(A), args...; kwargs...)

function primitive_array(A::AbstractArray{C}) where C <: Colorant{T,3} where T <: FixedPoint{F} where F
    # 1. Grab channel view which gets us Array{T,3}. This is a Base.ReinterpretArray
    # 2. Convert to Array{F,3}, still a Base.ReinterpretArray
    # 3. Permute dims such that the channels are the last dimension, without copying
    PermutedDimsArray( reinterpret(F,channelview(A)), (2,3,1))
end

view_image(A::AbstractArray{C}, args...; kwargs...) where C <: Colorant{T,4} where T <: FixedPoint{F} where F =
    view_image(primitive_array(A), args...; kwargs...)
add_image(viewer, A::AbstractArray{C}, args...; kwargs...) where C <: Colorant{T,4} where T <: FixedPoint{F} where F =
    add_image(viewer, primitive_array(A), args...; kwargs...)

function primitive_array(A::AbstractArray{C}) where C <: Colorant{T,4} where T <: FixedPoint{F} where F
    # 1. Grab channel view which gets us Array{T,3}. This is a Base.ReinterpretArray
    # 2. Convert to Array{F,3}, still a Base.ReinterpretArray
    # 3. Permute dims such that the channels are the last dimension, without copying
    PermutedDimsArray( reinterpret(F,channelview(A)), (2,3,1))
end

view_image(A::AbstractArray{C}, args...; kwargs...) where C <: Colorant{T,1} where T =
    view_image(primitive_array(A), args...; kwargs...)
add_image(viewer, A::AbstractArray{C}, args...; kwargs...) where C <: Colorant{T,1} where T =
    add_image(viewer, primitive_array(A), args...; kwargs...)

function primitive_array(A::AbstractArray{C}) where C <: Colorant{T,1} where T <: FixedPoint{F} where F
    # 1. Grab channel view which gets us Array{T,3}. This is a Base.ReinterpretArray
    # 2. Convert to Array{F,3}, still a Base.ReinterpretArray
    reinterpret(F,channelview(A))
end
function primitive_array(A::AbstractArray{C}) where C <: Colorant{F, 1} where F <: AbstractFloat
    reinterpret(F,channelview(A))
end

view_image(A::AbstractArray{C}, args...; kwargs...) where C <: TransparentColor{Gray{F}, F} where F =
    view_image(primitive_array(A), args...; kwargs...)
add_image(A::AbstractArray{C}, args...; kwargs...) where C <: TransparentColor{Gray{F}, F} where F =
    add_image(primitive_array(A), args...; kwargs...)

function primitive_array(A::AbstractArray{C}) where C <: TransparentColor{Gray{F}, F} where F
    # Napari does not know how to combine alpha channel and gray
    primitive_array( RGBA.(A) )
end

"""
    view_image(img::ImageMeta, ...)

    This function adds the following keywords if not provided:
        channel_axis = axis with name :channel
        order = gives axes in order of :time, :channel, :y, :x
        axis_labels = axisnames(img)
        metadata = properties(img)
"""
function view_image(img::ImageMeta, args...; kwargs...)
    # ImageMeta might be from OMETIFF
    # For example, testimage("multi-channel-time-series.ome.tif")
    kwdict = Dict{Symbol,Any}(kwargs)

    if !haskey(kwdict, :metadata)
        kwdict[:metadata] = properties(img)
    end

    view_image( arraydata(img) , args... ; kwdict...)
end

"""
    add_image(viewer, img::ImageMeta, ...)

    This function adds the following keywords if not provided:
        channel_axis = axis with name :channel
        metadata = properties(img)

    Unlike view_image(img::ImageMeta, ...) it does change the
    order or axis_labels
"""
function add_image(viewer, img::ImageMeta, args...; kwargs...)
    # ImageMeta might be from OMETIFF
    # For example, testimage("multi-channel-time-series.ome.tif")
    kwdict = Dict{Symbol,Any}(kwargs)

    if !haskey(kwdict, :metadata)
        kwdict[:metadata] = properties(img)
    end

    add_image(viewer, arraydata(img) , args... ; kwdict...)
end

"""
    view_image(img::AxisArray{Gray{T}}, ...)

    This function adds the following keywords if not provided:
        channel_axis = axis with name :channel
        order = gives axes in order of :time, :channel, :y, :x
        axis_labels = axisnames(img)
"""
function view_image(img::AxisArray, args...; kwargs...)
    A = primitive_array( img )
    dims = collect( 1:ndims( img ) )
    dim_names = axisnames(img)
    time_dim = findfirst(in((:time, :T)), dim_names)
    channel_dim = findfirst(in((:channel, :C)), dim_names)
    kwdict = Dict{Symbol,Any}(kwargs)
    # If we have an image with (:y, :x, :channel, :time)
    # transform it to (:time, :channel, :y, :x)
    if !isnothing(time_dim)
        popat!(dims, time_dim)
        if !isnothing(channel_dim)
            popat!(dims, channel_dim)
            pushfirst!(dims, channel_dim)
            if !haskey(kwdict, :channel_axis)
                kwdict[:channel_axis] = channel_dim - 1
            end
        end
        pushfirst!(dims, time_dim)
        if !haskey(kwdict, :order)
           kwdict[:order] = dims .- 1
        end
    else
        # No time axis found
        if !isnothing(channel_dim)
            popat!(dims, channel_dim)
            pushfirst!(dims, channel_dim)
            if !haskey(kwdict, :channel_axis)
                kwdict[:channel_axis] = channel_dim - 1
            end
        end
        if !haskey(kwdict, :order)
            order = dims .- 1
        end
    end
    if !haskey(kwdict, :axis_labels)
        kwdict[:axis_labels] = String.(dim_names)
    end
    @info kwdict

    view_image(A, args... ; kwdict...)
end

"""
    add_image(viewer, img::AxisArrays{Gray{T},...})

    This function adds the channel_axis keyword if not provided.

    Unlike view_image(viewer, img::AxisArrays{Gray{T},...}) it does
    not change the order or add axis_labels.
"""
function add_image(viewer, img::AxisArray, args...; kwargs...)
    A = primitive_array( img )
    dims = collect( 1:ndims( img ) )
    dim_names = axisnames(img)
    time_dim = findfirst(in((:time, :T)), dim_names)
    channel_dim = findfirst(in((:channel, :C)), dim_names)
    kwdict = Dict{Symbol,Any}(kwargs)
    # If we have an image with (:y, :x, :channel, :time)
    # transform it to (:time, :channel, :y, :x)
    if !isnothing(time_dim)
        popat!(dims, time_dim)
        if !isnothing(channel_dim)
            popat!(dims, channel_dim)
            pushfirst!(dims, channel_dim)
            if !haskey(kwdict, :channel_axis)
                # We are going to permute the array
                kwdict[:channel_axis] = 1
            end
        end
        pushfirst!(dims, time_dim)
        # We permute the array since we are adding to an existing viewer
        @info "Permuting Dims"
        A = PermutedDimsArray(A, dims)
    else
        # No time axis found
        if !isnothing(channel_dim)
            popat!(dims, channel_dim)
            pushfirst!(dims, channel_dim)
            if !haskey(kwdict, :channel_axis)
                kwdict[:channel_axis] = channel_dim - 1
            end
        end
    end

    add_image(viewer, A, args... ; kwdict...)
end

function primitive_array(A::AxisArray{Gray{T},N}) where T <: FixedPoint{F} where {F,N}
    reinterpret(F, arraydata(A))
end
function primitive_array(A::AxisArray{T}) where T
    reinterpret(T, arraydata(A))
end

# Use NumPyArray for types with strides and a parent

view_image(img::PermutedDimsArray, args...; kwargs...) = view_image(NumPyArray(img), args...; kwargs...)
add_image(viewer, img::PermutedDimsArray, args...; kwargs...) = add_image(viewer, NumPyArray(img), args...; kwargs...)
view_image(img::SubArray, args...; kwargs...) = view_image(NumPyArray(img), args...; kwargs...)
add_image(viewer, img::SubArray, args...; kwargs...) = add_image(viewer, NumPyArray(img), args...; kwargs...)
view_image(img::Base.ReinterpretArray, args...; kwargs...) = view_image(NumPyArray(img), args...; kwargs...)
add_image(viewer, img::Base.ReinterpretArray, args...; kwargs...) = add_image(viewer, NumPyArray(img), args...; kwargs...)
view_image(img::Base.ReshapedArray, args...; kwargs...) = view_image(NumPyArray(img), args...; kwargs...)
add_image(viewer, img::Base.ReshapedArray, args...; kwargs...) = add_image(viewer, NumPyArray(img), args...; kwargs...)

end # PyConvert module end


text:
module MakieGallery

using AbstractPlotting
using Documenter
using Documenter: Selectors, Expanders
using Markdown
using Markdown: Link, Paragraph
using ImageTransformations, FileIO
using ImageFiltering  # needed for Gaussian-filtering images during resize
using Random
using Documenter.Writers
using Documenter.Writers.HTMLWriter
using Test, Statistics
using FFMPEG
using BinaryProvider

include("documenter_extension.jl")
include("database.jl")
include("io.jl")
include("gallery.jl")
include("visualregression.jl")
include("testruns.jl")

function event_path(entry, ending)
    joinpath(@__DIR__, "..", "examples", "recorded_events", string(entry.unique_name, ".jls"))
end

"""
Loads the full example database and returns it!
"""
function load_database()
    return load_database([
        "tutorials.jl",
        "examples2d.jl",
        "examples3d.jl",
        "interactive.jl",
        "documentation.jl",
        "diffeq.jl",
        "implicits.jl",
        "short_tests.jl",
        "recipes.jl",
        "bigdata.jl",
        "layouting.jl",
        "statsmakie.jl",
    ])
end

"""
    load_tests()
Loads a database of minimal tests and returns it!
"""
function load_tests()
    return load_database([
        "tutorials.jl",
        "short_tests.jl",
        "layouting.jl",
    ])
end

"""
    `load_database(files::Array{String, 1})`

Loads a database with the files given as a String array.
They must be filenames `filename.jl`, and must be located
within `joinpath(dirname(pathof(MakieGallery)), "..")`.
(On Linux, this would translate to `MakieGallery/examples/`).
"""
function load_database(files::AbstractVector{<: AbstractString})
    empty!(unique_names)
    empty!(database)
    dir = abspath(joinpath(dirname(pathof(MakieGallery)), ".."))
    nfiles = joinpath.("$dir", "examples", files)
    for file in nfiles
        MakieGallery.eval(:(include($file)))
    end
    database
end

export load_database, eval_example, available_examples, run_example

end