text:
using Images, ImageView, ImageSegmentation

# convert image segments to mean color value
segment_to_image(segments) = map(i->segment_mean(segments, i), labels_map(segments))

# one object
img = load("sample-images/cat-3352842_640.jpg");
imshow(img);

# find segments
segments_1 = fast_scanning(img, 0.05)
segments_2 = fast_scanning(img, 0.15)
segments_3 = fast_scanning(img, 0.2)

# define helper values
color_black = RGB4(0.,0.,0.)
img_width = size(img, 2)

# create a new image which will be a stack
# of 3 different results of felzenszwalb algorithm
new_img = fill(color_black, size(img) .* (1, 3))
new_img[:, 1:img_width] = segment_to_image(segments_1)
new_img[:, img_width+1:img_width*2] = segment_to_image(segments_2)
new_img[:, img_width*2+1:img_width*3] = segment_to_image(segments_3)
new_img[:, img_width] = new_img[:, img_width*2] = RGB4(0.,0.,0.)

# prune segments
deletion_rule = i -> (segment_pixel_count(segments_3,i) < 750)
replacement_rule = (i,j) -> (-segment_pixel_count(segments_3, j))
segments_n = prune_segments(segments_3, deletion_rule, replacement_rule)

# preview the result
segment_to_image(segments_n)
# save("cat-felzenszwalb-123.jpg", restrict(new_img))

# multiple objects
img = load("sample-images/bird-3183441_640.jpg");
#imshow(img);

# find segments
segments_1 = fast_scanning(img, 0.05)
segments_2 = fast_scanning(img, 0.15)
segments_3 = fast_scanning(img, 0.2)

# define helper values
color_black = RGB4(0.,0.,0.)
img_width = size(img, 2)

# create a new image which will be a stack
# of 3 different results of felzenszwalb algorithm
new_img = fill(color_black, size(img) .* (1, 3))
new_img[:, 1:img_width] = segment_to_image(segments_1)
new_img[:, img_width+1:img_width*2] = segment_to_image(segments_2)
new_img[:, img_width*2+1:img_width*3] = segment_to_image(segments_3)
new_img[:, img_width] = new_img[:, img_width*2] = RGB4(0.,0.,0.)

deletion_rule = i -> (segment_pixel_count(segments_2,i) < 750)
replacement_rule = (i,j) -> (-segment_pixel_count(segments_2, j))
segments_n = prune_segments(segments_2, deletion_rule, replacement_rule)

# preview the result

imshow(new_img)


text:
@testset "ArrayDictionary" begin
    @test ArrayDictionary() isa ArrayDictionary{Any, Any}

    d = ArrayDictionary{Int64, Int64}()

    @test issettable(d)
    @test isinsertable(d)
    @test istokenizable(d)
    @test length(d) == 0
    @test isempty(d)
    @test isempty(keys(d))
    @test isequal(copy(d), d)
    @test_throws IndexError d[10]
    @test get(d, 10, 15) == 15
    @test length(unset!(d, 10)) == 0
    io = IOBuffer(); print(io, d); @test String(take!(io)) == "{}"
    if VERSION < v"1.6-"
        io = IOBuffer(); show(io, MIME"text/plain"(), d); @test String(take!(io)) == "0-element ArrayDictionary{Int64,Int64,ArrayIndices{Int64,Array{Int64,1}},Array{Int64,1}}"
    else
        io = IOBuffer(); show(io, MIME"text/plain"(), d); @test String(take!(io)) == "0-element ArrayDictionary{Int64, Int64, ArrayIndices{Int64, Vector{Int64}}, Vector{Int64}}"
    end
    @test_throws IndexError d[10] = 11
    @test_throws IndexError delete!(d, 10)

    insert!(d, 10, 11)

    @test d[10] == 11
    @test haskey(d, 10)
    @test !haskey(d, 100)
    @test get(d, 10, 15) == 11
    @test get!(d, 10, 15) == 11
    @test length(d) == 1
    @test keys(d) === keys(keys(d))
    @test !isempty(d)
    @test isequal(copy(d), d)
    @test_throws IndexError insert!(d, 10, 12)
    @test d[10] == 11
    set!(d, 10, 12)
    @test length(d) == 1
    @test d[10] == 12
    d[10] = 13
    @test length(d) == 1
    @test d[10] == 13
    io = IOBuffer(); print(io, d); @test String(take!(io)) == "{10 = 13}"
    if VERSION < v"1.6-"
        io = IOBuffer(); show(io, MIME"text/plain"(), d); @test String(take!(io)) == "1-element ArrayDictionary{Int64,Int64,ArrayIndices{Int64,Array{Int64,1}},Array{Int64,1}}\n 10 │ 13"
    else
        io = IOBuffer(); show(io, MIME"text/plain"(), d); @test String(take!(io)) == "1-element ArrayDictionary{Int64, Int64, ArrayIndices{Int64, Vector{Int64}}, Vector{Int64}}\n 10 │ 13"
    end
    @test !isequal(d, empty(d))
    @test isequal(d, copy(d))
    @test isempty(empty(d))

    delete!(d, 10)
    @test isequal(d, ArrayDictionary{Int64, Int64}())

    @test get!(d, 10, 14) == 14
    @test d[10] == 14
    delete!(d, 10)

    for i = 2:2:100
        insert!(d, i, i+1)
        @test d[i] == i + 1 ? true : (@show i; false)
    end
    @test all(in(i, d) == !iseven(i) for i in 2:2:100)
    @test all(in(i, keys(d)) == iseven(i) for i in 2:2:100)
    @test isempty(empty!(d))

    @testset "rand" begin
        dict = ArrayDictionary(["a", "b", "c", "d", "e"], 1:5)
        for i = 1:100
            @test rand(dict) in 1:5
        end
    end
end

ext: jl


text:
#!/usr/bin/env julia
using PlasmaTools
using Base.Test

# Testing gs2ls
g = GS2geometry(equilibrium_option="eik",rhoc=0.6,irho=2,qinp=2.0,shat=1.0,shift=-0.1,akappa=0.2,rmaj=3.0,r_geo=2.9)

r = GS2resolution(nperiod=3,ky=0.2)

p = GS2params(runname="gs2ls_test",beta=0.003)

specs = GS2species[]
push!(specs, GS2species(dens=0.5,temp=1.0,mass=1.0,z=1.0,tprim=3.0,fprim=1.0,vnewk=0.01))
push!(specs, GS2species(dens=0.5,temp=1.0,mass=1.5,z=1.0,tprim=3.0,fprim=1.0,vnewk=0.01))
push!(specs, GS2species(dens=1.0,temp=1.1,mass=0.0002,z=-1.0,tprim=2.5,fprim=1.0,vnewk=0.6))

createInputFile("gs2ls_test.in",p,r,g,specs)
info("Wrote a gs2 input file: gs2ls_test.in. Check for correctness.")

# Resolution for tests
Nv = 800

# Test conditions
n_cgs = 1.0e14
T_eV = 10.0e3

T = T_eV*el
n = n_cgs*1.0e6

deuterium = SpeciesData(density=0.5*n,temperature=T,charge=el,mass=2.0*mp)
tritium = SpeciesData(density=0.5*n,temperature=T,charge=el,mass=3.0*mp)
electron = SpeciesData(density=n,temperature=T,charge=-el,mass=me)
helium3 = SpeciesData(density=0.01*n,temperature=T,charge=2.0*el,mass=3.0*mp)

bulkspecs = [deuterium,tritium,electron]

vti = sqrt(2.0*T/mp)
vte = sqrt(2.0*T/me)
vmax = 5.0*vti
vgrid = createFVgrid(Nv,0.0,vmax,[1,0])

f0trace = helium3.density*(2.0*pi*helium3.temperature/helium3.mass)^(-1.5)*exp(-0.5*helium3.mass*vgrid.ctr.^2/helium3.temperature)

# Test value of lnLambda, from NRL plasma formulary
lnLambda_ee = 23.5 - log(sqrt(n_cgs)*T_eV^-1.25) - sqrt(1e-5 + (log(T_eV)-2)^2/16)
@test_approx_eq(lnLambda(electron,electron),lnLambda_ee)
lnLambda_ei = 24.0-log(sqrt(n_cgs)/T_eV)
@test_approx_eq(lnLambda(electron,helium3),lnLambda_ei)
@test_approx_eq(lnLambda(helium3,electron),lnLambda_ei)
lnLambda_ii = 23.0-log( (2/T_eV)*sqrt( (0.01*n_cgs*4/T_eV) + (0.5*n_cgs/T_eV)) )
@test_approx_eq(lnLambda(helium3,deuterium),lnLambda_ii)
@test_approx_eq(lnLambda(deuterium,helium3),lnLambda_ii)

# Test value of nuhat against deuterium
@test_approx_eq(nuhat(helium3,deuterium),10.411163139347321*lnLambda_ii)

# Test low-v limit of Chandrasekhar function
x = 1.0e-4
@test_approx_eq_eps(gfunc(x),2*x/(3.0*sqrt(pi)),1.0e-12)

# Test high-v limit of Chandrasekhar function
x = 8.0
@test_approx_eq(gfunc(x),0.5*x^-2)

# Test value of nu_s

# Test value of nu_par

# Test slowing-down time

# Test solver for simple problem, advection only
f0anal = exp(-vgrid.ctr.^2/vti^2)
A0 = 5.0
D0 = 1.0
onefunc(x) = 1.0
Afunc(x) = A0
Dfunc(x) = D0
source = f0anal.*( (D0/vti^2)*(4.0*(vgrid.ctr.^2/vti^2) - 2.0) - 2.0*vgrid.ctr*A0/vti^2)

operator = createAdvecDiffOperator(vgrid,onefunc,Afunc,Dfunc,[0.0,0.0],source)
f0test = operator\source
#println(f0test)
#println(f0anal)
@test_approx_eq_eps(f0test,f0anal,1.0e-3)

# Test solver for simple problem, diffusion only

# Test that collision operator vanishes for Maxwellian
collop = collisionOperator(vgrid,helium3,bulkspecs)
testfunc = collop*f0trace
@test_approx_eq_eps(testfunc[1:Nv],zeros(Nv),1.0e-3)

# Test that solving collision operator gives Maxwellian

println("All tests passed.")

max_stars_repo_name: gjwilkie/PlasmaTools

max_issues_repo_name: gjwilkie/PlasmaTools

max_forks_repo_name: gjwilkie/PlasmaTools


text:
module SurfaceFlow

using DocStringExtensions
using UnPack
using ..Land
using ..VariableTemplates
using ..BalanceLaws
import ..BalanceLaws:
    BalanceLaw,
    prognostic_vars,
    flux,
    source,
    precompute,
    eq_tends,
    vars_state,
    Prognostic,
    Auxiliary,
    Gradient,
    GradientFlux

using ...DGMethods: LocalGeometry, DGModel
using StaticArrays: SVector

export OverlandFlowModel,
    NoSurfaceFlowModel,
    surface_boundary_flux!,
    surface_boundary_state!,
    calculate_velocity,
    Precip,
    VolumeAdvection,
    SurfaceWaterHeight

"""
    SurfaceWaterHeight <: AbstractPrognosticVariable

The prognostic variable type for the 2d overland flow model. Used only for
dispatching on.
"""
struct SurfaceWaterHeight <: AbstractPrognosticVariable end

"""
    NoSurfaceFlowModel <: BalanceLaw

The default surface flow model, which does not add any prognostic variables
to the land model and therefore does not model surface flow.
"""
struct NoSurfaceFlowModel <: BalanceLaw end

"""
    OverlandFlowModel{Sx,Sy,M} <: BalanceLaw

The 2D overland flow model, with a prognostic variable
equal to the height of the surface water.

This model simulates the depth-averaged shallow water equation under the
kinematic approximation, and employs Manning's relationship to relate
velocity to the height of the water.
# Fields
$(DocStringExtensions.FIELDS)
"""
struct OverlandFlowModel{Sx, Sy, M} <: BalanceLaw
    "Slope in x direction; field(x,y), unitless"
    slope_x::Sx
    "Slope in y direction; field(x,y), unitless"
    slope_y::Sy
    "Mannings coefficient; field(x,y), units of s/m^(1/3)"
    mannings::M
end

function OverlandFlowModel(
    slope_x::Function,
    slope_y::Function;
    mannings::Function = (x, y) -> convert(eltype(x), 0.03),
)
    args = (slope_x, slope_y, mannings)
    return OverlandFlowModel{typeof.(args)...}(args...)
end

"""
    calculate_velocity(surface, x::Real, y::Real, h::Real)

Given the surface flow model, calculate the velocity of the flow
based on height, slope, and Manning's coefficient.
"""
function calculate_velocity(surface, x::Real, y::Real, h::Real)
    FT = eltype(h)
    sx = FT(surface.slope_x(x, y))
    sy = FT(surface.slope_y(x, y))
    mannings_coeff = FT(surface.mannings(x, y))
    coeff = h^FT(2 / 3) / mannings_coeff
    #The velocity direction is opposite the slope vector (∂_x z, ∂_y z)
    return SVector(
        -sign(sx) * coeff * sqrt(abs(sx)),
        -sign(sy) * coeff * sqrt(abs(sy)),
        zero(FT),
    )
end

vars_state(surface::OverlandFlowModel, st::Prognostic, FT) = @vars(height::FT)

function Land.land_init_aux!(
    land::LandModel,
    surface::Union{NoSurfaceFlowModel, OverlandFlowModel},
    aux,
    geom::LocalGeometry,
) end

function Land.land_nodal_update_auxiliary_state!(
    land::LandModel,
    surface::Union{NoSurfaceFlowModel, OverlandFlowModel},
    state,
    aux,
    t,
) end

"""
    VolumeAdvection <: TendencyDef{Flux{FirstOrder}}

A first order flux type for the overland flow model.
"""
struct VolumeAdvection <: TendencyDef{Flux{FirstOrder}} end

"""
    flux(::SurfaceWaterHeight, ::VolumeAdvection, land::LandModel, args,)

A first order flux method for the OverlandFlow model, adding in advection of water volume.
"""
function flux(::SurfaceWaterHeight, ::VolumeAdvection, land::LandModel, args)
    @unpack state, aux = args
    x = aux.x
    y = aux.y
    height = max(state.surface.height, eltype(aux)(0.0))
    v = calculate_velocity(land.surface, x, y, height)
    return height * v
end

# Boundary Conditions

# General case - to be used with bc::NoBC
function surface_boundary_flux!(
    nf,
    bc::Land.AbstractBoundaryConditions,
    m::Union{NoSurfaceFlowModel, OverlandFlowModel},
    land::LandModel,
    _...,
) end

function surface_boundary_state!(
    nf,
    bc::Land.AbstractBoundaryConditions,
    m::Union{NoSurfaceFlowModel, OverlandFlowModel},
    land::LandModel,
    _...,
) end

"""
    surface_boundary_flux!(
        nf,
        bc::Land.Dirichlet,
        model::OverlandFlowModel,
        land::LandModel,
        _...,
    )

The surface boundary flux function for the OverlandFlow model, which
does nothing if Dirichlet conditions are chosen in order to not
overconstrain the solution.
"""
function surface_boundary_flux!(
    nf,
    bc::Land.Dirichlet,
    model::OverlandFlowModel,
    land::LandModel,
    _...,
) end

"""
    surface_boundary_state!(
        nf,
        bc::Land.Dirichlet,
        model::OverlandFlowModel,
        land::LandModel,
        state⁺::Vars,
        aux⁺::Vars,
        nM,
        state⁻::Vars,
        aux⁻::Vars,
        t,
        _...,
    )

The surface boundary state function for the OverlandFlow model when
Dirichlet conditions are chosen. This should be equivalent to
outflow boundary conditions, also referred to as gradient outlet
conditions.
"""
function surface_boundary_state!(
    nf,
    bc::Land.Dirichlet,
    model::OverlandFlowModel,
    land::LandModel,
    state⁺::Vars,
    aux⁺::Vars,
    nM,
    state⁻::Vars,
    aux⁻::Vars,
    t,
    _...,
)
    bc_function = bc.state_bc
    state⁺.surface.height = bc_function(aux⁻, t)
end

"""
    Precip <: TendencyDef{Source}

A source term for overland flow where a prescribed net precipitation
drives the overland flow.
"""
struct Precip{FT, F} <: TendencyDef{Source}
    precip::F

    function Precip{FT}(precip::F) where {FT, F}
        new{FT, F}(precip)
    end
end

function (p::Precip{FT})(x, y, t) where {FT}
    FT(p.precip(x, y, t))
end

prognostic_vars(::Precip) = (SurfaceWaterHeight(),)

precompute(source_type::Precip, land::LandModel, args, tt::Source) =
    NamedTuple()

function source(::SurfaceWaterHeight, s::Precip, land::LandModel, args)
    @unpack aux, t = args
    return s(aux.x, aux.y, t)
end

end


text:
using BasisFunctions, BasisFunctions.Test, DomainSets, GridArrays

using Test, LinearAlgebra, Plots

delimit("Plots")
@testset begin
    e = random_expansion(Fourier(4))
    plot(e; plot_complex=false)
    plot(e; plot_complex=true)
    plot(e, exp)

    plot(random_expansion(Fourier(4)⊗Fourier(4)))

    plot(Fourier(4))

    plot(Fourier(4)[1])

    plot(Fourier(4) → -3..5)

    plot(MultiDict((Fourier(4),Fourier(4))))
end


text:
#####
##### There's no free surface with a rigid lid
#####

@inline explicit_barotropic_pressure_x_gradient(i, j, k, grid, free_surface::Nothing) = zero(eltype(grid))
@inline explicit_barotropic_pressure_y_gradient(i, j, k, grid, free_surface::Nothing) = zero(eltype(grid))


text:
"""
    function crsp_stocknames(
        conn::Union{LibPQ.Connection, DBInterface.Connection};
        cols::Array{String}=["permno", "cusip", "ncusip", "comnam", "namedt", "nameenddt", "ticker"]
    )

    function crsp_stocknames(
        conn::Union{LibPQ.Connection, DBInterface.Connection},
        cusip::Array{String};
        cols::Array{String}=["permno", "cusip", "ncusip", "comnam", "namedt", "nameenddt", "ticker"],
        cusip_col="cusip", # either "cusip" "ncusip" or "ticker"
    )

    function crsp_stocknames(
        conn::Union{LibPQ.Connection, DBInterface.Connection},
        permno::Array{<:Number};
        cols::Array{String}=["permno", "cusip", "ncusip", "comnam", "namedt", "nameenddt", "ticker"],
    )

Download crsp.stockname data (with non-missing ncusip).
If an array of strings is passed, by default assumes it is a list of cusips, can be a list of
ncusip or tickers (change `cusip_col`). Can also take an array of numbers which is assumed to be
a list of permnos.

"""
function crsp_stocknames(
    conn::Union{LibPQ.Connection, DBInterface.Connection};
    cols::Array{String}=["permno", "cusip", "ncusip", "comnam", "namedt", "nameenddt", "ticker"],
    filters::Dict{String, <:Any}=Dict(
        "ncusip" => missing
    )
)
    col_str = join(cols, ", ")
    fil = create_filter(filters)

    query = """
        SELECT DISTINCT $col_str
        FROM $(default_tables["crsp_stocknames"])
        $fil
    """
    return run_sql_query(conn, query) |> DataFrame
end

function crsp_stocknames(
    conn::Union{LibPQ.Connection, DBInterface.Connection},
    cusip::Array{String};
    cols::Array{String}=["permno", "cusip", "ncusip", "comnam", "namedt", "nameenddt", "ticker"],
    cusip_col="cusip", # either "cusip" "ncusip" or "ticker"
    filters::Dict{String, <:Any}=Dict(
        "ncusip" => missing
    )
)
    if cusip_col ∉ ["cusip", "ncusip", "ticker"]
        @error("`cusip_col` must be one of \"cusip\", \"ncusip\" or \"ticker\"")
    end
    if 0 < length(cusip) <= 1000
        temp_filters = Dict{String, Any}()
        for (key, val) in filters
            temp_filters[key] = val
        end
        filters = temp_filters
        filters[cusip_col] = cusip
    end
    return crsp_stocknames(
        conn;
        cols,
        filters
    )
end

function crsp_stocknames(
    conn::Union{LibPQ.Connection, DBInterface.Connection},
    permno::Array{<:Number};
    cols::Array{String}=["permno", "cusip", "ncusip", "comnam", "namedt", "nameenddt", "ticker"],
    filters::Dict{String, <:Any}=Dict(
        "ncusip" => missing
    )
)
    if 0 < length(permno) <= 1000
        temp_filters = Dict{String, Any}()
        for (key, val) in filters
            temp_filters[key] = val
        end
        filters = temp_filters
        filters["permno"] = permno
    end
    return crsp_stocknames(
        conn;
        cols,
        filters
    )
end

"""
    function crsp_market(
        conn::Union{LibPQ.Connection, DBInterface.Connection},
        dateStart::Union{Date,Int}=1950,
        dateEnd::Union{Date,Int}=Dates.today();
        col::Union{String, Array{String}}="vwretd"
    )

Downloads the data from the daily or monthly stock index file (dsi and msi) for a range of dates with one value
for each day (with various return columns). Available columns are:

- "vwretd": Value weighted return with dividends
- "vwretx": Value weighted return without dividends
- "ewretd": Equal weighted return with dividends
- "ewretx": Equal weighted return without dividends
- "sprtrn": S&P500 return
"""
function crsp_market(
    conn::Union{LibPQ.Connection, DBInterface.Connection},
    dateStart::Union{Date,Int}=1950,
    dateEnd::Union{Date,Int}=Dates.today();
    cols::Union{String, Array{String}}="vwretd"
)

    if typeof(dateStart) == Int
        dateStart = Dates.Date(dateStart, 1, 1)
    end
    if typeof(dateEnd) == Int
        dateEnd = Dates.Date(dateEnd, 12, 31)
    end

    cols = if typeof(cols) <: String
        [cols, "date"]
    elseif "date" ∉ cols
        vcat(cols, ["date"])
    else
        cols
    end

    col_str = join(cols |> unique, ", ")

    query = """
                        select $col_str
                        from $(default_tables["crsp_index"])
                        where date between '$dateStart' and '$dateEnd'
                        """
    return run_sql_query(conn, query) |> DataFrame
end

# This set of functions is for compiling a series of SQL statements
# or for estimating the amount of time it would take to download a table
# of a given size.

date_filter(date_range::StepRange) = "date BETWEEN '$(date_range[1])' AND '$(date_range[end])'"
date_filter(date_range::AbstractArray{<:StepRange}) = "((" * join(date_filter.(date_range), ") OR (") * "))"
date_filter(dates::AbstractArray{Date}) = "date IN ('$(join(dates, "','"))')"

function sql_query_basic(
    permno::Int,
    dates;
    cols::Vector{String}=["permno", "date", "ret"]
)
    "SELECT $(join(cols, ",")) FROM $(default_tables["crsp_stock_data"]) WHERE permno = $permno AND $(date_filter(dates))"
end

function sql_query_full(
    permnos,
    dates;
    cols=["permno", "date", "ret"]
)
    #"(" *
    join(sql_query_basic.(permnos, dates; cols), " UNION ")
    # * ")"
end

function merge_date_ranges(x::AbstractArray{<:StepRange})
    x = sort(x)
    out = eltype(x)[]
    s = x[1][1]
    e = x[1][end]
    for i in 2:length(x)
        if e >= x[i][1]
            e = x[i][end]
        else
            push!(out, s:Day(1):e)
            s = x[i][1]
            e = x[i][end]
        end
    end
    push!(out, s:Day(1):e)
    out
end

"""
    function crsp_data(
        conn::Union{LibPQ.Connection, DBInterface.Connection},
        [permnos::Vector{<:Real},]
        s::Date=Date(1925),
        e::Date=today();
        cols = ["ret", "vol", "shrout"],
        filters::Dict{String, <:Any}=Dict{String, Any}(),
        adjust_crsp_data::Bool=true
    )

    function crsp_data(
        conn::Union{LibPQ.Connection, DBInterface.Connection},
        permnos::Vector{<:Real},
        dates::Vector{Date},
        [dates_end::Vector{Date}];
        cols=["ret", "vol", "shrout"],
        adjust_crsp_data::Bool=true,
        query_size_limit::Int=3000
    )

Downloads data from the crsp stockfiles, which are individual stocks. To download the data from
the monthly stockfile, change the default table to the monthly stockfile:
```julia
WRDSMerger.default_tables["crsp_stock_data"] = "crsp.msf"
WRDSMerger.default_tables["crsp_index"] = "crsp.msi"
WRDSMerger.default_tables["crsp_delist"] = "crsp.msedelist"
```

# Arguments

- `permnos::Vecotr{<:Real}`: A vector of firm IDs, if provided, will only download data for those firms
- `s::Date=Date(1925)` and `e::Date=today()`: Downloads all data between two dates provided
- `dates::Vector{Date}`: Downloads data for a set of permnos on the date provided
    - `dates_end::Vector{Date}`: If provided, then treats the `dates` as the start of a period
      and will download data for the permnos between the two dates
- `adjust_crsp_data::Bool=true`: This will call `crsp_adjust` with all options
    set to true, it will only do the operations that it has the data for.
"""
function crsp_data(
    conn::Union{LibPQ.Connection, DBInterface.Connection},
    permnos::Vector{<:Real},
    date_start::Date=Date(1926),
    date_end::Date=today();
    cols = ["ret", "vol", "shrout"],
    adjust_crsp_data::Bool=true,
    filters::Dict{String, <:Any} = Dict{String, Any}()
)
    @assert all(isinteger.(permnos)) "All of the Permnos must be convertable to an Integer"

    filters["permno"] = Int.(permnos)
    return crsp_data(
        conn,
        date_start,
        date_end;
        cols,
        filters,
        adjust_crsp_data
    )
end

function crsp_data(
    conn::Union{LibPQ.Connection, DBInterface.Connection},
    s::Date=Date(1925),
    e::Date=today();
    cols = ["ret", "vol", "shrout"],
    filters::Dict{String, <:Any}=Dict{String, Any}(),
    adjust_crsp_data::Bool=true
)
    for col in ["permno", "date"]
        if col ∉ cols
            push!(cols, col)
        end
    end
    colString = join(cols, ", ")

    filter_str = create_filter(filters, "WHERE date BETWEEN '$s' AND '$e'")

    query = """
        select $colString
        from $(default_tables["crsp_stock_data"])
        $filter_str
    """
    crsp = run_sql_query(conn, query) |> DataFrame
    if adjust_crsp_data
        crsp = crsp_adjust(conn, crsp)
    end
    return crsp
end

function crsp_data(
    conn::Union{LibPQ.Connection, DBInterface.Connection},
    permnos::Vector{<:Real},
    dates::Vector{Date};
    cols=["ret", "vol", "shrout"],
    adjust_crsp_data::Bool=true,
    query_size_limit::Int=3000
)

    @assert all(isinteger.(permnos))  "All of the Permnos must be convertable to an Integer"
    @assert length(permnos) == length(dates) "Must have same number of dates as Permnos"

    for col in ["permno", "date"]
        if col ∉ cols
            push!(cols, col)
        end
    end

    permnos = Int.(permnos)

    df = DataFrame(
        permno=permnos,
        date=dates
    ) |> unique

    if length(unique(dates)) < nrow(df) / 5
        temp_market = crsp_market(dsn, minimum(dates), maximum(dates); cols=["dates"])
        df = innerjoin(
            df,
            temp_market,
            on=:date,
            validate=(false, true)
        )
    end

    df = combine(
        groupby(df, :permno),
        :date => (x -> [Vector(x)])
    )
    rename!(df, [:permno, :date])

    queries = sql_query_full.(
        collect(Iterators.partition(df.permno, query_size_limit)),
        collect(Iterators.partition(df.date, query_size_limit));
        cols
    )
    crsp = DataFrame()
    for q in queries
        temp = run_sql_query(conn, q)
        if nrow(crsp) == 0
            crsp = temp[:, :]
        else
            crsp = vcat(crsp, temp)
        end
    end

    if adjust_crsp_data
        crsp = crsp_adjust(conn, crsp)
    end
    return crsp
end

function crsp_data(
    conn::Union{LibPQ.Connection, DBInterface.Connection},
    permnos::Vector{<:Real},
    dates_min::Vector{Date},
    dates_max::Vector{Date};
    cols=["ret", "vol", "shrout"],
    adjust_crsp_data::Bool=true,
    query_size_limit::Int=3000
)
    @assert all(isinteger.(permnos))  "All of the Permnos must be convertable to an Integer"
    @assert length(permnos) == length(dates_min) == length(dates_max) "Must have same number of dates as Permnos"

    for col in ["permno", "date"]
        if col ∉ cols
            push!(cols, col)
        end
    end

    permnos = Int.(permnos)

    df = DataFrame(
        permno=permnos
    )
    df[!, :date_range] = [d1:Day(1):d2 for (d1, d2) in zip(dates_min, dates_max)]
    df = combine(
        groupby(df, :permno),
        "date_range" => x -> [merge_date_ranges(x)]
    )
    rename!(df, [:permno, :date_range])

    queries = sql_query_full.(
        collect(Iterators.partition(df.permno, query_size_limit)),
        collect(Iterators.partition(df.date_range, query_size_limit));
        cols
    )
    crsp = DataFrame()
    for q in queries
        temp = run_sql_query(conn, q)
        if nrow(crsp) == 0
            crsp = temp[:, :]
        else
            crsp = vcat(crsp, temp)
        end
    end

    if adjust_crsp_data
        crsp = crsp_adjust(conn, crsp)
    end
    return crsp
end

function crsp_data(
    conn::Union{LibPQ.Connection, DBInterface.Connection},
    df::DataFrame;
    cols = ["ret", "vol", "shrout"],
    pull_method::Symbol=:optimize, # :optimize, :minimize, :stockonly, :alldata,
    date_start::String="dateStart",
    date_end::String="dateEnd",
    adjust_crsp_data::Bool=true
)
    df = df[:, :]
    for col in ["permno", date_start, date_end]
        if col ∉ names(df)
            @error("$col must be in the DataFrame")
        end
    end
    pull_methods = [:optimize, :minimize, :stockonly, :alldata]
    if pull_method ∉ pull_methods
        @error("pull_method must be one of $pull_methods")
    end

    for col in ["permno", "date"]
        if col ∉ cols
            push!(cols, col)
        end
    end

    if pull_method == :optimize || pull_method == :minimize
        return crsp_data(
            conn,
            df.permno,
            df[:, date_start],
            df[:, date_end];
            cols,
            adjust_crsp_data
        )
    elseif pull_method == :stockonly
        return crsp_data(
            conn,
            df.permno |> unique,
            minimum(df[:, date_start]),
            maximum(df[:, date_end]);
            cols,
            adjust_crsp_data
        )
    else
        return crsp_data(
            conn,
            minimum(df[:, date_start]),
            maximum(df[:, date_end]);
            cols,
            adjust_crsp_data
        )
    end
end

"""
    function crsp_delist(
        conn::Union{LibPQ.Connection, DBInterface.Connection};
        cols::Array{String}=[
            "permno",
            "dlstdt",
            "dlret"
        ],
        date_start::Date=Date(1926),
        date_end::Date=today()
    )

Fetches the CRSP delist dataset, typically for the returns
on the day of delisting.
"""
function crsp_delist(
    conn::Union{LibPQ.Connection, DBInterface.Connection},
    date_start::Date=Date(1926),
    date_end::Date=today();
    cols::Array{String}=[
        "permno",
        "dlstdt",
        "dlret"
    ],
    filters::Dict{String, <:Any}=Dict(
        "dlret" => missing
    )
)
    fil = create_filter(filters, "WHERE dlstdt BETWEEN '$date_start' AND '$date_end' AND dlret != 0")
    col_str = join(cols, ", ")
    query = """
    SELECT $col_str FROM $(default_tables["crsp_delist"]) $fil
    """
    return run_sql_query(conn, query) |> DataFrame
end

"""
    function crsp_adjust(
        conn,
        df::DataFrame;
        kwargs...
    )

This makes 4 common adjustments to CRSP data:
1. When the price in CRSP is negative, that means there was not a
   specific close price. Typically, researchers take the absolute value
   to get the actual price.
2. Prices are not adjusted for splits (returns are). CRSP includes
   the number (cfascpr) that will adjust prices to be comparable through time
3. Similar to prices, shares outstanding are not adjusted. The number used to
   adjust is different due to various events.
4. Prices are not adjusted for delisting, so this downloads the
   necessary dataset and adjusts returns accordingly.

# Arguments

## Options on what to adjust

- adjust_prc_negatives::Bool=true: Corresponds to (1) above
- adjust_prc_splits::Bool=true: Corresponds to (2) above
- adjust_shr_splits::Bool=true: Corresponds to (3) above
- adjust_delist::Bool=true: Corresponds to (4) above

## Pre-existing column names

- date::String="date": date column
- idcol::String="permno": primary identifier (permno) column
- prc_col::String="prc": price column
- ret_col::String="ret": return column
- prc_splits_col::String="cfacpr": price split adjustment factor
- shrout_col::String="shrout": shares outstanding
- shrout_splits_col::String="cfacshr": share split adjustment factor

## Options to relabel columns

- adjusted_neg_prc_col::String="prc": relabel prices to this after taking absolute value
- adjusted_prc_col::String="prc_adj": relabel prices to this after adjusting for splits
- adjusted_shrout_col::String="shrout_adj":relabel shares outstanding to this after adjusting for splits

"""
function crsp_adjust(
    conn,
    df::DataFrame;
    adjust_prc_negatives::Bool=true,
    adjust_prc_splits::Bool=true,
    adjust_shr_splits::Bool=true,
    adjust_delist::Bool=true,
    date::String="date",
    idcol::String="permno",
    prc_col::String="prc",
    adjusted_neg_prc_col::String="prc",
    prc_splits_col::String="cfacpr",
    adjusted_prc_col::String="prc_adj",
    shrout_col::String="shrout",
    shrout_splits_col::String="cfacshr",
    adjusted_shrout_col::String="shrout_adj",
    ret_col::String="ret"
)
    df = copy(df)
    if adjust_prc_negatives && prc_col ∈ names(df)
        df[!, adjusted_neg_prc_col] = abs.(df[!, prc_col])
    end

    if adjust_prc_splits && prc_col ∈ names(df) && prc_splits_col ∈ names(df)
        df[!, adjusted_prc_col] = abs.(df[:, prc_col]) ./ df[:, prc_splits_col]
    end

    if adjust_shr_splits && shrout_col ∈ names(df) && shrout_splits_col ∈ names(df)
        df[!, adjusted_shrout_col] = df[:, shrout_col] .* df[:, shrout_splits_col]
    end

    if adjust_delist && ret_col ∈ names(df)
        delist_ret = crsp_delist(
            conn,
            minimum(df[:, date]),
            maximum(df[:, date])
        ) |> dropmissing |> unique
        df = leftjoin(
            df,
            delist_ret,
            on=[idcol => "permno", date => "dlstdt"],
            validate=(true, true)
        )
        df[!, :dlret] = coalesce.(df.dlret, 0.0)
        df[!, ret_col] = (1 .+ df[:, ret_col]) .* (1 .+ df.dlret) .- 1
        select!(df, Not(:dlret))
    end
    return df
end


text:
"""
For testing purposes - simply has a single model
"""
type RobustAdapter{S,A} <: RobustMDP{S,A}
    mdp::MDP{S,A}
end

rand_model(rmdp::RobustAdapter, rng::AbstractRNG) = rmdp.mdp

representative_mdp(r::RobustAdapter) = r.mdp


text:
module BenchSpanTree
import Random
import GraphIdx

function grid_tree(
    n1::Integer,
    n2::Integer,
    weights::Vector{Float64} = Float64[],
    root::Integer = 1
)
    if length(weights) < 2 * n1 * n2
        @time weights = randn(2 * n1 * n2)
    end
    graph = GraphIdx.Grid.GridGraph(n1, n2)
    @assert GraphIdx.num_edges(graph) <= length(weights)
    @time mem = GraphIdx.PrimMstMem(graph)
    @time tree = GraphIdx.prim_mst_edges(weights, root, mem)
end

function benchmark(n1)
    n = n1^2
    @info "normal" n
    w = randn(Random.MersenneTwister(42), 2n)
    @time pi = BenchSpanTree.grid_tree(n1, n1, w)
    println(pi[1:5] .- 1, "…", pi[end-4:end] .- 1)
end

end

import Random

@info "precompile"
@time BenchSpanTree.grid_tree(2, 2)

@info "uniform"
w = rand(Random.MersenneTwister(42), 4_000_000)
@time BenchSpanTree.grid_tree(1000, 1000, w)

for n1 in [1000, 2000, 3000]
    BenchSpanTree.benchmark(n1)
end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
module SparseAlgebra

include("basis.jl")

export *

end


text:
"""
        color_graph(g::LightGraphs.AbstractGraphs, ::AcyclicColoring)

Returns a coloring vector following the acyclic coloring rules (1) the coloring
corresponds to a distance-1 coloring, and (2) vertices in every cycle of the
graph are assigned at least three distinct colors. This variant of coloring is
called acyclic since every subgraph induced by vertices assigned any two colors
is a collection of trees—and hence is acyclic.

Reference: Gebremedhin AH, Manne F, Pothen A. **New Acyclic and Star Coloring Algorithms with Application to Computing Hessians**
"""
function color_graph(g::LightGraphs.AbstractGraph, ::AcyclicColoring)

    color = zeros(Int, nv(g))
    forbidden_colors = zeros(Int, nv(g))

    set = DisjointSets{LightGraphs.Edge}([])

    first_visit_to_tree = Array{Tuple{Int, Int}, 1}(undef, ne(g))
    first_neighbor = Array{Tuple{Int, Int}, 1}(undef, nv(g))

    for v in vertices(g)
        #enforces the first condition of acyclic coloring
        for w in outneighbors(g, v)
            if color[w] != 0
                forbidden_colors[color[w]] = v
            end
        end
        #enforces the second condition of acyclic coloring
        for w in outneighbors(g, v)
            if color[w] != 0 #colored neighbor
                for x in outneighbors(g, w)
                    if color[x] != 0 #colored x
                        if forbidden_colors[color[x]] != v
                            prevent_cycle(v, w, x, g, color, forbidden_colors, first_visit_to_tree, set)
                        end
                    end
                end
            end
        end

        color[v] = min_index(forbidden_colors, v)

        #grow star for every edge connecting colored vertices v and w
        for w in outneighbors(g, v)
            if color[w] != 0
                grow_star!(set, v, w, g, first_neighbor, color)
            end
        end

        #merge the newly formed stars into existing trees if possible
        for w in outneighbors(g, v)
            if color[w] != 0
                for x in outneighbors(g, w)
                    if color[x] != 0 && x != v
                        if color[x] == color[v]
                            merge_trees!(set, v, w, x, g)
                        end
                    end
                end
            end
        end
    end

    return color
end

"""
        prevent_cycle(v::Integer,
                    w::Integer,
                    x::Integer,
                    g::LightGraphs.AbstractGraph,
                    color::AbstractVector{<:Integer},
                    forbidden_colors::AbstractVector{<:Integer},
                    first_visit_to_tree::Array{Tuple{Integer, Integer}, 1},
                    set::DisjointSets{LightGraphs.Edge})

Subroutine to avoid generation of 2-colored cycle due to coloring of vertex v,
which is adjacent to vertices w and x in graph g. Disjoint set is used to store
the induced 2-colored subgraphs/trees where the id of set is a key edge of g
"""
function prevent_cycle(v::Integer,
                        w::Integer,
                        x::Integer,
                        g::LightGraphs.AbstractGraph,
                        color::AbstractVector{<:Integer},
                        forbidden_colors::AbstractVector{<:Integer},
                        first_visit_to_tree::AbstractVector{<: Tuple{Integer, Integer}},
                        set::DisjointSets{LightGraphs.Edge})

    edge = find_edge(g, w, x)
    e = find_root(set, edge)
    p, q = first_visit_to_tree[edge_index(g, e)]
    if p != v
        first_visit_to_tree[edge_index(g, e)] = (v, w)
    elseif q != w
        forbidden_colors[color[x]] = v
    end
end

"""
        min_index(forbidden_colors::AbstractVector{<:Integer}, v::Integer)

Returns min{i > 0 such that forbidden_colors[i] != v}
"""
function min_index(forbidden_colors::AbstractVector{<:Integer}, v::Integer)
    return findfirst(!isequal(v), forbidden_colors)
end

"""
        grow_star!(set::DisjointSets{LightGraphs.Edge},
                v::Integer,
                w::Integer,
                g::LightGraphs.AbstractGraph
                first_neighbor::Array{Tuple{Integer, Integer}, 1})

Subroutine to grow a 2-colored star after assigning a new color to the
previously uncolored vertex v, by comparing it with the adjacent vertex w.
Disjoint set is used to store stars in sets, which are identified through key
edges present in g.
"""
function grow_star!(set::DisjointSets{LightGraphs.Edge},
                   v::Integer,
                   w::Integer,
                   g::LightGraphs.AbstractGraph,
                   first_neighbor::AbstractArray{<: Tuple{Integer, Integer}, 1},
                   color::AbstractVector{<: Integer})
    edge = find_edge(g, v, w)
    push!(set, edge)
    p, q = first_neighbor[color[w]]
    if p != v
        first_neighbor[color[w]] = (v, w)
    else
        edge1 = find_edge(g, v, w)
        edge2 = find_edge(g, p, q)
        e1 = find_root(set, edge1)
        e2 = find_root(set, edge2)
        union!(set, e1, e2)
    end
    return nothing
end

"""
        merge_trees!(v::Integer,
                w::Integer,
                x::Integer,
                g::LightGraphs.AbstractGraph,
                set::DisjointSets{LightGraphs.Edge})

Subroutine to merge trees present in the disjoint set which have a
common edge.
"""
function merge_trees!(set::DisjointSets{LightGraphs.Edge},
                    v::Integer,
                    w::Integer,
                    x::Integer,
                    g::LightGraphs.AbstractGraph)
    edge1 = find_edge(g, v, w)
    edge2 = find_edge(g, w, x)
    e1 = find_root(set, edge1)
    e2 = find_root(set, edge2)
    if (e1 != e2)
        union!(set, e1, e2)
    end
end

"""
        find_edge(g::LightGraphs.AbstractGraph, v::Integer, w::Integer)

Returns an edge object of the type LightGraphs.Edge which represents the
edge connecting vertices v and w of the undirected graph g
"""
function find_edge(g::LightGraphs.AbstractGraph,
                   v::Integer,
                   w::Integer)
    for e in edges(g)
        if (src(e) == v && dst(e) == w) || (src(e) == w && dst(e) == v)
            return e
        end
    end
    throw(ArgumentError("$v and $w are not connected in graph g"))
end

"""
        edge_index(g::LightGraphs.AbstractGraph, e::LightGraphs.Edge)

Returns an Integer value which uniquely identifies the edge e in graph
g. Used as an index in main function to avoid custom arrays with non-
numerical indices.
"""
function edge_index(g::LightGraphs.AbstractGraph,
                    e::LightGraphs.Edge)
    for (i, edge) in enumerate(edges(g))
        if edge == e
            return i
        end
    end
    throw(ArgumentError("Edge $e is not present in graph g"))
end


text:
JSONRPC.parse_params(::Type{Val{Symbol("workspace/didChangeWatchedFiles")}}, params) = DidChangeWatchedFilesParams(params)
function process(r::JSONRPC.Request{Val{Symbol("workspace/didChangeWatchedFiles")},DidChangeWatchedFilesParams}, server)
    for change in r.params.changes
        uri = change.uri

        startswith(uri, "file:") || continue

        if change.type == FileChangeTypes.Created || change.type == FileChangeTypes.Changed
            if hasdocument(server, URI2(uri))
                doc = getdocument(server, URI2(uri))

                # Currently managed by the client, we don't do anything
                if get_open_in_editor(doc)
                    continue
                else
                    filepath = uri2filepath(uri)
                    content = try
                        s = read(filepath, String)
                        if !isvalid(s)
                            deletedocument!(server, URI2(uri))
                            continue
                        end
                        s
                    catch err
                        isa(err, Base.IOError) || isa(err, Base.SystemError) || rethrow()
                        deletedocument!(server, URI2(uri))
                        continue
                    end

                    set_text!(doc, content)
                    set_is_workspace_file(doc, true)
                    parse_all(doc, server)
                end
            else
                filepath = uri2filepath(uri)
                content = try
                    s = read(filepath, String)
                    isvalid(s) || continue
                    s
                catch err
                    isa(err, Base.IOError) || isa(err, Base.SystemError) || rethrow()
                    continue
                end

                doc = Document(uri, content, true, server)
                setdocument!(server, URI2(uri), doc)
                parse_all(doc, server)
            end
        elseif change.type == FileChangeTypes.Deleted
            if hasdocument(server, URI2(uri))
                doc = getdocument(server, URI2(uri))

                # We only handle if currently not managed by client
                if !get_open_in_editor(doc)
                    deletedocument!(server, URI2(uri))

                    publishDiagnosticsParams = PublishDiagnosticsParams(uri, missing, Diagnostic[])
                    JSONRPCEndpoints.send_notification(server.jr_endpoint, "textDocument/publishDiagnostics", publishDiagnosticsParams)
                else
                    # TODO replace with accessor function once the other PR
                    # that introduces the accessor is merged
                    doc._workspace_file = false
                end
            end
        else
            error("Unknown change type.")
        end
    end
end

JSONRPC.parse_params(::Type{Val{Symbol("workspace/didChangeConfiguration")}}, params) = params
function process(r::JSONRPC.Request{Val{Symbol("workspace/didChangeConfiguration")},Dict{String,Any}}, server::LanguageServerInstance)
    request_julia_config(server)
end

function request_julia_config(server::LanguageServerInstance)
    server.clientCapabilities.workspace.configuration === false && return # Or !== true?

    response = JSONRPCEndpoints.send_request(server.jr_endpoint, "workspace/configuration", ConfigurationParams([
        ConfigurationItem(missing, "julia.format.indent"), # FormatOptions
        ConfigurationItem(missing, "julia.format.indents"),
        ConfigurationItem(missing, "julia.format.ops"),
        ConfigurationItem(missing, "julia.format.tuples"),
        ConfigurationItem(missing, "julia.format.curly"),
        ConfigurationItem(missing, "julia.format.calls"),
        ConfigurationItem(missing, "julia.format.iterOps"),
        ConfigurationItem(missing, "julia.format.comments"),
        ConfigurationItem(missing, "julia.format.docs"),
        ConfigurationItem(missing, "julia.format.lineends"),
        ConfigurationItem(missing, "julia.format.kw"),
        ConfigurationItem(missing, "julia.lint.call"), # LintOptions
        ConfigurationItem(missing, "julia.lint.iter"),
        ConfigurationItem(missing, "julia.lint.nothingcomp"),
        ConfigurationItem(missing, "julia.lint.constif"),
        ConfigurationItem(missing, "julia.lint.lazy"),
        ConfigurationItem(missing, "julia.lint.datadecl"),
        ConfigurationItem(missing, "julia.lint.typeparam"),
        ConfigurationItem(missing, "julia.lint.modname"),
        ConfigurationItem(missing, "julia.lint.pirates"),
        ConfigurationItem(missing, "julia.lint.useoffuncargs"),
        ConfigurationItem(missing, "julia.lint.run"),
        ConfigurationItem(missing, "julia.lint.missingrefs")
        ]))

    server.format_options = DocumentFormat.FormatOptions(response[1:11]...)
    server.runlinter = something(response[22], true)
    server.lint_missingrefs = Symbol(something(response[23], :all))

    new_SL_opts = StaticLint.LintOptions(response[12:21]...)
    # TODO: implement == for StaticLint.LintOptions
    rerun_lint = any(getproperty(server.lint_options, opt) != getproperty(new_SL_opts, opt) for opt in fieldnames(StaticLint.LintOptions))
    server.lint_options = new_SL_opts

    if rerun_lint
        for doc in getdocuments_value(server)
            lint!(doc, server)
        end
    end

end

JSONRPC.parse_params(::Type{Val{Symbol("workspace/didChangeWorkspaceFolders")}}, params) = DidChangeWorkspaceFoldersParams(params)
function process(r::JSONRPC.Request{Val{Symbol("workspace/didChangeWorkspaceFolders")}}, server)
    for wksp in r.params.event.added
        push!(server.workspaceFolders, uri2filepath(wksp.uri))
        load_folder(wksp, server)
    end
    for wksp in r.params.event.removed
        delete!(server.workspaceFolders, uri2filepath(wksp.uri))
        remove_workspace_files(wksp, server)
    end
end

JSONRPC.parse_params(::Type{Val{Symbol("workspace/symbol")}}, params) = WorkspaceSymbolParams(params)
function process(r::JSONRPC.Request{Val{Symbol("workspace/symbol")},WorkspaceSymbolParams}, server)
    syms = SymbolInformation[]
    for doc in getdocuments_value(server)
        bs = collect_toplevel_bindings_w_loc(getcst(doc), query = r.params.query)
        for x in bs
            p, b = x[1], x[2]
            push!(syms, SymbolInformation(valof(b.name), 1, false, Location(doc._uri, Range(doc, p)), missing))
        end
    end

    return syms
end


text:
module HttpUtil

using ..Coincheck
using ..Methods
using HTTP

function make_http_request(method, url; headers = Nullable(), body = Nullable())
    const _body = (method == Coincheck.Methods.GET) ? Nullable() : body
    const make_request =
        if (method == Coincheck.Methods.GET) HTTP.get
        elseif (method == Coincheck.Methods.POST) HTTP.post
        elseif (method == Coincheck.Methods.DELETE) HTTP.delete
        else nothing
        end
    if isnull(headers) && isnull(body)
        make_request(url)
    elseif isnull(headers)
        make_request(url, body = body)
    elseif isnull(body)
        make_request(url, headers = headers)
    else
        make_request(url, headers = headers, body = body)
    end
end
export make_http_request

function convert_to_url(endpoint, path, args = Nullable())
    _endpoint = isdirpath(endpoint) ? dirname(endpoint) : endpoint
    _path = isdirpath(path) ? dirname(path) : path
    _path = isabspath(_path) ? _path[2:_path.len] : _path
    if isnull(args)
        return "$_endpoint/$_path"
    else
        query = join(map(arg -> "$(arg[1])=$(arg[2])", collect(args)), "&")
        return "$_endpoint/$_path?$query"
    end
end
export convert_to_url

end


text:
mutable struct ArrayType <: SchemaType
    type::String
    items::Schema
    default::Any
end

==(x::ArrayType, y::ArrayType) = x.items == y.items
ArrayType(items::Schema, default=nothing) = ArrayType("array", items, default)
ArrayType() = ArrayType("array", "", nothing)
StructTypes.StructType(::Type{ArrayType}) = StructTypes.Mutable()

struct Array{T} <: AbstractVector{T}
    data::Vector{UInt8}
    values::Vector{Int}
end

Base.IndexStyle(::Type{<:Array}) = Base.IndexLinear()
Base.size(x::Array) = (length(x.values),)
Base.getindex(x::Array{T}, i::Int) where {T} =
    readvalue(Binary(), T, x.data, x.values[i], length(x.data), (;))[1]

function Base.copy(x::Array{T}) where {T}
    len = length(x)
    y = Vector{T}(undef, len)
    for i = 1:len
        @inbounds y[i] = copy(x[i])
    end
    return y
end

juliatype(x::ArrayType) = Array{juliatype(x.items)}
schematype(::StructTypes.ArrayType, ::Type{A}) where {A} = ArrayType(schematype(eltype(A)))

function writevalue(B::Binary, A::ArrayType, x, buf, pos, len, opts)
    xlen = length(x)
    pos = writevalue(B, long, -xlen, buf, pos, len, opts)
    xlen == 0 && return pos
    nb = nbytes(A, x, true)
    pos = writevalue(B, long, nb, buf, pos, len, opts)
    startpos = pos
    for val in x
        pos = writevalue(B, A.items, val, buf, pos, len, opts)
    end
    @assert nb == (pos - startpos)
    # write out final 0 count block
    pos = writevalue(B, long, 0, buf, pos, len, opts)
    return pos
end

function nbytes(A::ArrayType, x, self=false)
    nb = isempty(x) ? 0 : sum(y -> nbytes(A.items, y), x)
    return self ? nb : nbytes(length(x)) + nbytes(nb) + nb + 1
end

function readvalue(B::Binary, AT::ArrayType, ::Type{A}, buf, pos, buflen, opts) where {A <: AbstractVector{T}} where {T}
    sz = startpos = 0
    len, pos = readvalue(B, long, Int, buf, pos, buflen, opts)
    if len == 0
        return StructTypes.construct(A, Array{T}(buf, Int[])), pos
    elseif len < 0
        sz, pos = readvalue(B, long, Int, buf, pos, buflen, opts)
        startpos = pos
        len = -len
    end
    values = Vector{Int}(undef, len)
    for i = 1:len
        @inbounds values[i] = pos
        pos = skipvalue(B, AT.items, T, buf, pos, buflen, opts)
    end
    if sz > 0
        @assert sz == (pos - startpos)
    end
    sz = startpos = 0
    len, pos = readvalue(B, long, Int, buf, pos, buflen, opts)
    len == 0 && return StructTypes.construct(A, Array{T}(buf, values)), pos
    # multiple blocks
    arr = ChainedVector([StructTypes.construct(A, Array{T}(buf, values))])
    while true
        if len < 0
            sz, pos = readvalue(B, long, Int, buf, pos, buflen, opts)
            startpos = pos
            len = -len
        end
        values = Vector{Int}(undef, len)
        for i = 1:len
            @inbounds values[i] = pos
            pos = skipvalue(B, AT.items, T, buf, pos, buflen, opts)
        end
        if sz > 0
            @assert sz == (pos - startpos)
        end
        append!(arr, StructTypes.construct(A, Array{T}(buf, values)))
        len, pos = readvalue(B, long, Int, buf, pos, buflen, opts)
        len == 0 && return arr, pos
    end
end

function skipvalue(B::Binary, AT::ArrayType, ::Type{A}, buf, pos, buflen, opts) where {A <: AbstractVector{T}} where {T}
    while true
        len, pos = readvalue(B, long, Int, buf, pos, buflen, opts)
        if len == 0
            break
        elseif len < 0
            sz, pos = readvalue(B, long, Int, buf, pos, buflen, opts)
            pos += sz
        else
            for i = 1:-len
                pos = skipvalue(B, AT.items, T, buf, pos, buflen, opts)
            end
        end
    end
    return pos
end


text:
import MatFuncs
using Base.Test

function ident(i,j)
    if i == j
        1
    else
        0
    end
end
fm = MatFuncs.MatFunc(Int64, ident, (100,100))
@test fm[1,1] == 1
@test fm[1,2] == 0

function identlong(i,j)
    sleep(2)
    if i == j
        1
    else
        0
    end
end

fm2 = MatFuncs.MatFunc(Int64, identlong, (100000,100000))
@test fm2[1,1] == 1
@test fm2[1,2] == 0


text:
"""
    MenuAdventures.Testing

A sub-module for testing your universe.

```jldoctest
julia> using MenuAdventures

julia> using MenuAdventures.Testing

julia> @universe struct Universe <: AbstractUniverse
        end;

julia> @noun struct Room <: AbstractRoom
            already_lit::Bool = true
        end;

julia> @noun struct Person <: Noun
        end;

julia> function make_universe(interface)
            you = Person(
                "Brandon",
                grammatical_person = second_person,
                indefinite_article = "",
            )
            universe = Universe(you, interface = interface)
            universe[Room("room"), you] = Containing()
            universe
        end;

julia> cd(joinpath(pkgdir(MenuAdventures), "test", "Testing")) do
            check_choices(make_universe, transcript_file = "transcript1.txt")
        end
Conflict in line 6
Existing transcript: " > [control]quit[control]"
New transcript: " > [control]Quit[control]"
false

julia> cd(joinpath(pkgdir(MenuAdventures), "test", "Testing")) do
            check_choices(make_universe, transcript_file = "transcript2.txt")
        end
Conflict in line 6
Existing transcript: ""
New transcript: " > [control]Quit[control]"
false

julia> cd(joinpath(pkgdir(MenuAdventures), "test", "Testing")) do
            check_choices(make_universe, transcript_file = "transcript3.txt")
        end
Conflict in line 7
Existing transcript: " more"
New transcript: ""
false

julia> cd(joinpath(pkgdir(MenuAdventures), "test", "Testing")) do
            check_choices(make_universe, choices_file = "choices3.txt", )
        end
Unchosen choices: [1.0]
false
```
"""
module Testing

using DelimitedFiles: readdlm, writedlm
using MenuAdventures: AbstractUniverse, terminal, TTYTerminal, turn!

const KEY_PRESS = Dict(:up => "\e[A", :down => "\e[B", :enter => "\r")

function run_turn_sequence(io, choices_log)
    for option_number in choices_log
        for _ in 1:(option_number - 1)
            write(io, KEY_PRESS[:down])
        end
        write(io, KEY_PRESS[:enter])
    end
end

"""
    save_choices(make_universe;
        choices_file = "choices.txt",
        transcript_file = "transcript.txt",
        interface = terminal,
        resume = false
    )

Save the choices a user makes as a delimited file, as well as the transcript of the game.

Use to create a transcript to test with [`check_choices`](@ref).
`make_universe` should be a function which takes one argument, an IO interface, and returns an [`AbstractUniverse`](@ref).
`choices_file` will be the delimited file where user choices are saved.
`transcript_file` will be the file where the game transcript will be saved.
If `resume` is true, the game will pick up from the point you left off, based on the existing `choices_file`.
"""
function save_choices(make_universe;
    choices_file = "choices.txt",
    transcript_file = "transcript.txt",
    interface = terminal,
    resume = false
)
    take!(stdin.buffer)
    if resume
        run_turn_sequence(stdin.buffer, readdlm(choices_file))
    end
    universe = make_universe(interface)
    choices_log = universe.choices_log
    turn!(universe)
    open(io -> writedlm(io, choices_log), choices_file, "w")
    open(transcript_file, "w") do transcript
        take!(stdin.buffer)
        run_turn_sequence(stdin.buffer, choices_log)
        universe = make_universe(TTYTerminal("unix", stdin, transcript, stderr))
        turn!(universe)
    end
end

function find_line_difference(string1, string2)
    lines1 = split(string1, '\n')
    lines2 = split(string2, '\n')
    for (line_number, (line1, line2)) in enumerate(zip(lines1, lines2))
        if line1 != line2
            return line_number, line1, line2
        end
    end
    number_of_lines1 = length(lines1)
    number_of_lines2 = length(lines2)
    if number_of_lines2 > number_of_lines1
        next_line = number_of_lines1 + 1
        return next_line, "", lines2[next_line]
    elseif number_of_lines1 > number_of_lines2
        next_line = number_of_lines2 + 1
        return next_line, lines1[next_line], ""
    else
        return nothing
    end
end

export save_choices

function replace_escape_codes(text)
    replace(
        text,
        r"\x1b\[([\x30-\x3F]*)([\x20-\x2F]*)([\x40-\x7E])" =>
        s"[control]"
    )
end

"""
    check_choices(make_universe;
        choices_file = "choices.txt",
        transcript_file = "transcript.txt"
    )

Check the results of choices against a transcript created with [`save_choices`](@ref).

Return `true` if the results match, `false` otherwise.
`make_universe` should be a function which takes one argument, an IO interface, and returns an [`AbstractUniverse`].
`choices_file` should be the delimited file where user choices were saved.
`transcript_file` should be the file where the game transcript was saved.
"""
function check_choices(make_universe;
    choices_file = "choices.txt",
    transcript_file = "transcript.txt"
)
    output = IOBuffer()
    take!(stdin.buffer)
    universe = make_universe(TTYTerminal("unix", stdin, output, stderr))
    choices = readdlm(choices_file)
    run_turn_sequence(stdin.buffer, choices)
    turn!(universe)
    expected_choices = length(choices)
    choices_log = universe.choices_log
    actual_choices = length(choices_log)
    same_turn_count = actual_choices == expected_choices
    if !same_turn_count
        if actual_choices < expected_choices
            print("Unchosen choices: ")
            println(choices[actual_choices+1 : expected_choices])
        elseif expected_choices < actual_choices
            print("Extra choices: ")
            println(choices_log[expected_choices+1 : actual_choices])
        end
    end
    old_result = read(transcript_file, String)
    new_result = String(take!(output))
    difference = find_line_difference(old_result, new_result)
    same_transcript = difference === nothing
    if !same_transcript
        line_number, expected_line, received_line = difference
        print("Conflict in line ")
        println(line_number)
        print("Existing transcript: ")
        show(replace_escape_codes(expected_line))
        println()
        print("New transcript: ")
        show(replace_escape_codes(received_line))
        println()
    end
    same_turn_count && same_transcript
end

export check_choices

end


text:
"""
    {{ name }}({{ params }})

Lorem ipsum dolor sit amet.
"""
function {{ name }}({{ params }})

end


text:
function TC1(A)
    T = mul(A, A', mask=A, desc=S)
    y = reduce(+, T, dims=2)
    k = reduce(+, y)
    return (3 .* mul(A, y) - 2 .* mul(one.(T), y) .+ y) ./ k
end

function TC3(A)
    M = tril(A, -1)
    T = mul(A, A', (+, pair), mask=M, desc=S)
    y = reduce(+, T, dims=2) .+ reduce(+, T', dims=2, desc=S)
    k = reduce(+, y)
    T2 = *(+, second)(T, y) .+ *(+, second)(T', y, desc=S)
    r = (3 .* *(+, second)(A, y)) - (2 .* T2) + y
    return r ./ k
end


text:
"""
Initialize an model with zeros, given mgrid and names.
"""
function Base.zero(::Type{Medium}, mgrid, names=[:vp,:rho])
	nz=length(mgrid[1]); nx=length(mgrid[2])
	return Medium(mgrid,NamedArray([zeros(nz,nx) for i in names], (names,)),
		       NamedArray([fill(0.0,2) for i in names], (names,)),
		       NamedArray([0.0 for i in names],(names,)))
end
Medium(mgrid, names=[:vp,:rho])=zero(Medium, mgrid, names)
Medium(mgrid, m::Medium)=zero(Medium, mgrid, names(m.m)[1])

"""
Return true if a `Seismic` model is just allocated with zeros.
"""
function Base.iszero(mod::Medium)
	if (any(broadcast(iszero,mod.bounds)) && any(broadcast(iszero,mod.ref)))
		return true
	else
		for field in [:vp,:rho,:vs]
			if(field ∈ names(mod.m)[1])
				if(any(iszero.(mod[field])))
					return true
				end
			end
		end
	end
	return false
end

"""
Return true if a `Seismic` model is just allocated with zeros.
"""
function Base.isassigned(mod::Medium, fields=[:vp,:rho,:vs])
	for field in fields
		@assert field ∈ [:vp, :rho, :vs]
	end

	return any(broadcast(iszero,mod.bounds)) && any(broadcast(iszero,mod.ref))
end

"""
Return a similar model (size, bounds and reference values) to the input model, used for allocation.
"""
function Base.similar(mod::Medium)
	mod_new=deepcopy(mod)
	fill!(mod_new)
	return mod_new
end

"Compare if two `Seismic` models are equal"
function Base.isequal(mod1::T, mod2::T) where {T<:Union{Medium}}
	fnames = fieldnames(typeof(mod1))
	vec=[(isequal(getfield(mod1, name),getfield(mod2, name))) for name in fnames]
	return all(vec)
end

"""
Return if two `Seismic` models have same dimensions and bounds.
"""
function Base.size(mod::Medium)
	return length.(mod.mgrid)
end

"""
Copy for `Seismic` models. The models should have same bounds and sizes.
Doesn't allocate any memory.
Only copy primary medium parameters
"""
function Base.copyto!(modo::Medium, mod::Medium)
	for name in [:vp,:vs,:rho]
		((name ∈ names(modo.m)[1]) & (name ∈ names(mod.m)[1])) && (copyto!(modo.m[name], mod.m[name]) )
	end
	return modo
end

function Base.fill!(mod::Medium)
	for s in names(mod.m)[1]
		fill!(mod.m[s], mod.ref[s])
	end
	return mod
end

"""
Print information about `Seismic`
"""
function Base.show(io::Base.IO, mod::Medium)
	println(io, "> Medium")
	println(io, "> number of samples:\t","x\t",length(mod.mgrid[2]),"\tz\t",length(mod.mgrid[1]))
	println(io, "> sampling intervals:\t","x\t",step(mod.mgrid[2]),"\tz\t",step(mod.mgrid[1]))
	println(io, "> vp:\t","min\t",minimum(mod[:vp]),"\tmax\t",maximum(mod[:vp]))
	println(io, "> rho:\t","min\t",minimum(mod[:rho]),"\tmax\t",maximum(mod[:rho]))
	println(io, "> Bounds:")
	println(io, mod.bounds)
end

"""
Return `Seismic` with zeros everywhere;
this method is used for preallocation.

# Arguments
* `mgrid` : used for sizes of χ fields
"""
function Seismic_zeros(mgrid::Vector{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}})
	nz=length(mgrid[1]); nx=length(mgrid[2])
	return Seismic(mgrid,NamedArray([zeros(nz,nx),zeros(nx,nz)], ([:vp,:rho],)),
		       NamedArray([fill(0.0,2),fill(0.0,2)], ([:vp,:rho],)),
		       NamedArray([0.0,0.0],([:vp,:rho],)))
end

"""
Re-parameterization routine
that modifies the fields
`χvp` and `χrho` of an input seismic model
using two input vectors.

# Arguments

* `mod::Seismic` : to be updated
* `x1::Array{Float64,2}` : contrast of inverse bulk modulus
* `x2::Array{Float64,2}` : contrast of inverse density
* `fields:::Vector{Symbol}` : [:χKI, :χrhoI]
"""
function Base.copyto!(
	              mod::Medium,
		      x::AbstractVector,
		      fields::Vector{Symbol}=[:χKI, :χrhoI, :null]
		      )
	iszero(mod) ? error("mod cannot be zero") : nothing
	nznx=prod(length.(mod.mgrid))
	(length(x)≠(count(fields.≠ :null)*nznx)) &&  error("size x")
	if(fields == [:χKI, :χrhoI, :null])
		@inbounds for i in 1:nznx
			K=inv(χ(x[i],mod.ref[:KI],-1))
			rho=inv(χ(x[nznx+i],mod.ref[:rhoI],-1))
			mod.m[:vp][i]=sqrt(K*inv(rho))
			mod.m[:rho][i]=rho
		end
	elseif(fields == [:χKI, :null, :null])
		@inbounds for i in 1:nznx
			K=inv(χ(x[i],mod.ref[:KI],-1))
			rho=mod.m[:rho][i]
			mod.m[:vp][i]=sqrt(K*inv(rho))
			mod.m[:rho][i]=rho
		end
	elseif(fields == [:χvp, :χrho, :null])
		@inbounds for i in 1:nznx; mod.m[:vp][i]=χ(x[i],mod.ref[:vp],-1); end
		@inbounds for i in 1:nznx; mod.m[:rho][i]=χ(x[nznx+i],mod.ref[:rho],-1); end
	elseif(fields == [:null, :χrho, :null])
		@inbounds for i in 1:nznx; mod.m[:rho][i]=χ(x[i],mod.ref[:rho],-1); end
	elseif(fields == [:χvp, :null, :null])
		@inbounds for i in 1:nznx; mod.m[:vp][i]=χ(x[i],mod.ref[:vp],-1); end
	else
		error("invalid fields")
	end
	return mod
end

"""
Input perturbations corresponding to `parameterization`.
Output perturbations of KI and rhoI, i.e., without the contrast.
`mod` is just used for reference values.
"""
function Base.copyto!(δxout::AbstractVector, δx::AbstractVector, mod::Medium, parameterization)
	nznx=prod(length.(mod.mgrid))
	fill!(δxout, 0.0)
	if(parameterization == [:χKI, :χrhoI, :null])
		@inbounds for i in 1:nznx
			δxout[i]=δx[i]*mod.ref[:KI]
			δxout[nznx+i]=δx[nznx+i]*mod.ref[:rhoI]
		end
	elseif(parameterization == [:χKI, :null, :null])
		@inbounds for i in 1:nznx
			δxout[i]=δx[i]*mod.ref[:KI]
		end
	elseif(parameterization == [:χvp, :χrho, :null])
		@inbounds for i in 1:nznx
			δxout[i]=(-1.0*δx[nznx+i]-2.0*δx[i])*mod.ref[:KI]
			δxout[nznx+i]=(-1.0*δx[nznx+i])*mod.ref[:rhoI]
		end
	elseif(parameterization == [:χvp, :null, :null])
		@inbounds for i in 1:nznx
			δxout[i]=-2.0*δx[i]*mod.ref[:KI]
		end
	elseif(parameterization == [:null, :χrho, :null])
		@inbounds for i in 1:nznx
			δxout[i]=-1.0*δx[i]*mod.ref[:rhoI]
		end
	end
end

#pert_of_KI(dvp, drho, vp0, rho0) = - 2. / (rho0 * vp0^3) * dvp - 1 / (rho0^2 * vp0^2) * drho
#pert_of_rhoI(drho, vp0, rho0) = -1. / (rho0^2) * drho

"""
Get other dependent model parameters of a seismic model
that are not present in `Seismic`.

* `:rhoI` : inverse of density
* `:Zp` : P-wave impedance
"""
function Base.getindex(mod::Medium, s::Symbol)
	nznx=length(mod.m[:vp])
	# check if vp, vs, and rho are present
	@assert :vp ∈ names(mod.m)[1]
	@assert :rho ∈ names(mod.m)[1]
	vp=mod.m[:vp]; rho=mod.m[:rho]
	if(s == :vp)
		return vp
	elseif(s == :rho)
		return rho
	elseif(s in [:vs,:χmu])
		@assert :vs ∈ names(mod.m)[1]
		if(s==:vs)
			return mod.m[:vs]
		end
	end

	# add derived fields, if not present
	(s ∉ names(mod.m)[1]) && (mod.m=vcat(mod.m,NamedArray([zero(vp)],[[s],])))
	x=mod.m[s]

	# output derived fields
	if(s == :rhoI)
		@inbounds for i in 1:nznx; x[i]=inv(rho[i]); end
	elseif(s == :χrhoI)
		@inbounds for i in 1:nznx; x[i]=χ(inv(rho[i]),mod.ref[:rhoI],1); end
	elseif(s == :χrho)
		@inbounds for i in 1:nznx; x[i]=χ(rho[i],mod.ref[:rho],1); end
	elseif(s == :χvp)
		@inbounds for i in 1:nznx; x[i]=χ(vp[i],mod.ref[:vp],1); end
	elseif(s == :χK)
		@inbounds for i in 1:nznx; x[i]=χ(vp[i]*vp[i]*rho[i],mod.ref[:K],1); end
	elseif(s == :χmu)
		@inbounds for i in 1:nznx; x[i]=χ(mod.m[:vs][i]*mod.m[:vs][i]*rho[i],
			mod.ref[:mu],1); end
	elseif(s == :χKI)
		@inbounds for i in 1:nznx; x[i]=χ(inv(vp[i]*vp[i]*rho[i]),mod.ref[:KI],1); end
	elseif(s == :KI)
		@inbounds for i in 1:nznx; x[i]=inv(vp[i]*vp[i]*rho[i]); end
	elseif(s == :K)
		@inbounds for i in 1:nznx; x[i]=vp[i]*vp[i]*rho[i]; end
	elseif(s == :Zp)
		@inbounds for i in 1:nznx; x[i]=vp[i]*rho[i]; end
	else
		error(string("invalid attrib\t",s))
	end
	return mod.m[s]
end

"""
Copy medium parameters from Medium `mod`, to vector `x`. Return `x`.
copyto!(x,mod,[:vp])
"""
function Base.copyto!(x::AbstractArray{T}, mod::Medium, fields::Vector{Symbol}) where {T}
	rho=mod.m[:rho];
	vp=mod.m[:vp];
	nznx=length(rho)
	(length(x)≠(count(fields.≠ :null)*nznx)) &&  error("size x")
	i0=0;
	for field in fields
		if(field ≠:null)
			xx=view(x,i0+1:i0+nznx)
			m=Base.getindex(mod, field)
			copyto!(xx,m)
			i0+=nznx
		end
	end
	return x
end

"""
Return a vector of medium parameters `x`, from Medium `mod`.
x=vec(mod)
"""
function Base.vec(mod::Medium, field::Symbol)
	# allocate
	rout=zeros(length(mod.mgrid[1]), length(mod.mgrid[2]))
	copyto!(rout, mod, [field])
	return rout
end


text:
Core.eval(m::Module, x::Sym) = Core.eval(m, expr(x))

@inline expr(x::Sym; kwargs...) = expr(Val(gethead(x)), x; kwargs...)

function expr(::Val{:protected}, x::Sym; kwargs...)
    @assert gethead(x) === :protected
    ex = expr(firstarg(x); kwargs...)
    return :(protect($ex))
end

function expr(::Val{:object}, x::Sym; kwargs...)
    @assert gethead(x) === :object
    return expr(firstarg(x); kwargs...)
end

function expr(::Val{:symbol}, x::Sym{TAG}; annotate=false, type_annotate=false, kwargs...) where {TAG}
    @assert gethead(x) === :symbol
    ex = expr(firstarg(x); kwargs...)
    if annotate || (type_annotate && TAG !== Any)
        return :($ex::$TAG)
    else
        return ex
    end
end

function expr(::Val{:type}, x::Sym; kwargs...)
    @assert gethead(x) === :type
    T = expr(firstarg(x); kwargs...)
    return foldl((a, b) -> Expr(:., a, QuoteNode(b)), map(Symbol, split(repr(T), ".")))
end

function expr(::Val{:fn}, x::Sym; qualified_names=false, kwargs...)
    @assert gethead(x) === :fn
    fn = expr(firstarg(x); kwargs...)
    return foldl((a, b) -> Expr(:., a, QuoteNode(b)), map(Symbol, split(repr(fn), ".")))
end

function expr(::Val{:function}, x::Sym; kwargs...)
    @assert gethead(x) === :function
    fn, body = getargs(x)
    return Expr(:function, expr(fn, type_annotate=true; kwargs...), expr(body; kwargs...))
end

function expr(::Val{:struct}, x::Sym; kwargs...)
    @assert gethead(x) === :struct
    return Expr(
        :struct,
        map(arg -> expr(arg; kwargs...), getargs(x)[1:2])...,
        map(arg -> expr(arg, type_annotate=true; kwargs...), getargs(x)[3:end])...
    )
end

function expr(::Val{:(=)}, x::Sym; kwargs...)
    @assert gethead(x) === :(=)
    lhs, rhs = getargs(x)
    return Expr(:(=), expr(lhs, type_annotate=true; kwargs...), expr(rhs; kwargs...))
end

function expr(::Val{:(->)}, x::Sym; kwargs...)
    @assert gethead(x) === :(->)
    fn, body = getargs(x)
    return Expr(:(->), expr(fn, type_annotate=true; kwargs...), expr(body; kwargs...))
end

function expr(::Val{head}, x::Sym{TAG}; annotate=false, kwargs...) where {TAG,head}
    @assert gethead(x) === head
    ex = Expr(gethead(x), map(arg -> expr(arg; annotate=annotate, kwargs...), getargs(x))...)
    if annotate
        return :($ex::$TAG)
    else
        return ex
    end
end

@inline expr(a::Nothing; kwargs...) = Symbol(repr(a))
@inline expr(a; kwargs...) = a # Catch all

expr(t::Tuple; kwargs...) = Expr(:tuple, map(arg -> expr(arg; kwargs...), t)...)
expr(v::AbstractVector; kwargs...) = Expr(:vect, map(arg -> expr(arg; kwargs...), v)...)
expr(A::AbstractMatrix; kwargs...) = Expr(:vcat, mapslices(x -> Expr(:row, map(arg -> expr(arg; kwargs...), x)...), A, dims=2)...)
expr(x::LinearAlgebra.UniformScaling; kwargs...) = Expr(:call, :*, expr(x.λ; kwargs...), :I)

block(xs...) = Expr(:block, xs...)

getsymbols(x) = map(firstarg, getsyms(x))

getsyms(x) = Sym[]
getsyms(x::Sym) = getsyms(Val(gethead(x)), x)
getsyms(::Val, x::Sym) = mapreduce(getsyms, union, getargs(x))
getsyms(::Val{:symbol}, x::Sym) = Sym[x]
getsyms(x::Tuple) = mapreduce(getsyms, union, x)
getsyms(x::AbstractArray) = mapreduce(getsyms, union, x)

getops(x) = []
getops(x::Sym) = getops(Val(gethead(x)), x)
getops(::Val, x::Sym) = mapreduce(getops, union, getargs(x))
getops(::Val{:call}, x::Sym) = union([firstarg(x)], mapreduce(getops, union, tailargs(x)))
getops(x::Tuple) = mapreduce(getops, union, x)
getops(x::AbstractArray) = mapreduce(getops, union, x)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

# helper function to extract raw data
# from uniform scaling objects
raw(a::UniformScaling) = a.λ
raw(a) = a

"""
    NestedVariogram(cs, γs)

A nested variogram model `γ = c₁γ₁ + c₂γ₂ + ⋯ + cₙγₙ` with
coefficients `cs = (c₁, c₂, ..., cₙ)` and variogram models
`γs = (γ₁, γ₂, ..., γₙ)`.
"""
struct NestedVariogram{CS,GS} <: Variogram{Number,Metric}
  cs::CS
  γs::GS

  function NestedVariogram{CS,GS}(cs, γs) where {CS,GS}
    @assert all(issymmetric.(cs)) "coefficients must be symmetric"
    new(cs, γs)
  end
end

NestedVariogram(cs, γs) = NestedVariogram{typeof(cs),typeof(γs)}(cs, γs)

# variogram interface
(g::NestedVariogram)(h)                  = raw(sum(c*γ(h) for (c, γ) in zip(g.cs, g.γs)))
(g::NestedVariogram)(x::Point, y::Point) = raw(sum(c*γ(x, y) for (c, γ) in zip(g.cs, g.γs)))
sill(g::NestedVariogram)                 = raw(sum(c*sill(γ) for (c, γ) in zip(g.cs, g.γs)))
nugget(g::NestedVariogram)               = raw(sum(c*nugget(γ) for (c, γ) in zip(g.cs, g.γs)))
Base.range(g::NestedVariogram)           = maximum(range(γ) for γ in g.γs)
isstationary(g::NestedVariogram)         = all(isstationary(γ) for γ in g.γs)

# algebraic structure
*(c, γ::Variogram)                          = NestedVariogram((c,), (γ,))
*(c, γ::NestedVariogram)                    = NestedVariogram(map(x->c.*x, γ.cs), γ.γs)
+(γ₁::Variogram, γ₂::Variogram)             = NestedVariogram((1, 1), (γ₁, γ₂))
+(γ₁::NestedVariogram, γ₂::Variogram)       = NestedVariogram((γ₁.cs..., 1), (γ₁.γs..., γ₂))
+(γ₁::Variogram, γ₂::NestedVariogram)       = NestedVariogram((1, γ₂.cs...), (γ₁, γ₂.γs...))
+(γ₁::NestedVariogram, γ₂::NestedVariogram) = NestedVariogram((γ₁.cs..., γ₂.cs...), (γ₁.γs..., γ₂.γs...))

"""
    structures(γ)

Return the individual structures of a (possibly nested)
variogram as a tuple. The structures are the total nugget
`cₒ`, and the coefficients (or contributions) `cs` for the
remaining non-trivial structures `γs` after normalization
(i.e. sill=1, nugget=0).
"""
function structures(γ::Variogram)
  cₒ = nugget(γ)
  c  = sill(γ) - nugget(γ)
  T  = typeof(c)
  γ  = @set γ.sill = one(T)
  γ  = @set γ.nugget = zero(T)
  cₒ, (c,), (γ,)
end

function structures(γ::NestedVariogram)
  ks, gs = γ.cs, γ.γs

  # total nugget and contributions
  cₒ = raw(sum(@. ks * nugget(gs)))
  cs = @. raw(ks * (sill(gs) - nugget(gs)))

  # discard nugget effect terms
  inds = findall(g->!(g isa NuggetEffect), gs)
  cs, γs = cs[inds], gs[inds]

  # adjust sill and nugget
  γs = map(γs) do γ
    T = typeof(sill(γ))
    γ = @set γ.sill = one(T)
    γ = @set γ.nugget = zero(T)
  end

  cₒ, cs, γs
end

# -----------
# IO METHODS
# -----------

function Base.show(io::IO, g::NestedVariogram)
  N = length(g.cs)
  print(io, "NestedVariogram{$N}")
end

function Base.show(io::IO, ::MIME"text/plain", g::NestedVariogram)
  coeffs = 1 .* raw.(g.cs)
  models = [nameof(typeof(γ)) for γ in g.γs]
  params = ["range=$(range(γ)), sill=$(sill(γ)), nugget=$(nugget(γ))" for γ in g.γs]
  println(io, g)
  println(io, "  structures")
  lines  = ["    └─$γ($p)" for (γ, p) in zip(models, params)]
  println(io, join(lines, "\n"))
  println(io, "  coefficients")
  lines  = ["    └─$c" for c in coeffs]
  print(io, join(lines, "\n"))
end


text:
using Test, MPI
include(joinpath("..", "..", "testhelpers.jl"))

include("iterativesolvers.jl")


text:
include("minimal_sky.jl")
using Plots

grad =  ColorGradient([:white, :black])
strip2.Plot_Sky(NN, pix_size, conv_map, cmin=-100, cmax=50, gradient=grad)
m, h= strip2.binned_map(NN, tod, puntamenti)
strip2.Plot_Sky(NN, pix_size, m ./ h, cmin=-1000, cmax=1000, gradient=grad)

noise_ss = strip2.tod_noise(length(tod), 0.0083, 0.1, 3.0, 20.76)
#noise_ss2 = strip2.tod_noise(length(tod), 0.0083, 1.0, 2.0, 20.76)
dataset_baselines_2 = Array{Main.strip2.TOD_fake_pointing, 1}(undef, 1)
dat = strip2.TOD_fake_pointing(tod, puntamenti, noise_ss)
dataset_baselines_2[1] = dat

# Filtro semplice - In questo caso la sorgente cosmologica viene uccisa e portata
# molto vicino al livello di wite noise
@time tod_denoise = strip2.denoise(tod, noise_ss)
# tod_denoise2 = strip2.denoise(tod_denoise, noise_ss2)
mappa_bin, hit = strip2.binned_map(NN, tod_denoise, puntamenti)
#mappa_bin2, hit2 = strip2.binned_map(NN, tod_denoise2, puntamenti)
strip2.Plot_Sky(NN, pix_size, mappa_bin ./ hit, cmin=-0.1, cmax=0.05, gradient=grad)
#strip2.Plot_Sky(NN, pix_size, mappa_bin2 ./ hit2, cmin=-0.001, cmax=0.0005, gradient=grad)

# Ma se risolvo il problema del mapmaker usando la funzione strip2.denoise
# per approssimare la matrice di covarianza allora ottengo i livelli corretti
# Ci mette molto a convergere perche` non ho praticamente nessun cross-angle (worst-case)
map_des, exit_status, iter = strip2.destriper(NN, dataset_baselines_2)
map_d = reshape(map_des, NN, NN)
strip2.Plot_Sky(NN, pix_size, map_d, cmin=-100, cmax=50, gradient=grad)

# Se adotto una scanning strategy con dei cross-angles molto densa:
# Rifare i punti prima con una ss diversa.

tod_d = strip2.get_tod(map_d, puntamenti)
#dataset_baselines_3 = Array{Main.strip2.TOD_fake_pointing, 1}(undef, 1)
#dat2 = strip2.TOD_fake_pointing(tod_d, puntamenti, noise_ss2)
#dataset_baselines_3[1] = dat2

#map_des2, exit_status, iter = strip2.destriper(NN, dataset_baselines_3)
#map_d2 = reshape(map_des2, NN, NN)
#strip2.Plot_Sky(NN, pix_size, map_d2, cmin=-100, cmax=50, gradient=grad)
#tod_d2 = strip2.get_tod(map_d2, puntamenti)

using DSP

pp = DSP.welch_pgram(tod, 5000, 120)
pp2 = DSP.welch_pgram(tod_denoise, 5000, 120)
#pp3 = DSP.welch_pgram(tod_denoise2, 5000, 120)
pp4 = DSP.welch_pgram(tod_d, 5000, 120)
#pp5 = DSP.welch_pgram(tod_d2, 5000, 120)

plot(pp.freq[2:end], pp.power[2:end], xscale=:log, yscale=:log, label="Raw TOD")
plot!(pp2.freq[2:end], pp2.power[2:end], xscale=:log, yscale=:log, label="Filter1")
#plot!(pp3.freq[2:end], pp3.power[2:end], xscale=:log, yscale=:log, label="Filter2")
plot!(pp4.freq[2:end], pp4.power[2:end] * 10, xscale=:log, yscale=:log, label="MaxLike1")
#plot!(pp5.freq[2:end], pp5.power[2:end], xscale=:log, yscale=:log, label="MaxLike2")

using FFTW
using DSP
freq = strip2.fftfreq(length(tod), 0.0083)
fft_tod = fft(DSP.hamming(Int(length(tod))).*tod)
fft_tod_d = fft(DSP.hamming(length(tod)).*tod_d)

# A bare high-pass frequency filter
fft_filtered = fft(DSP.hamming(Int(length(tod_denoise))).*tod_denoise)

powe = real.(fft_tod).^2 .+ imag.(fft_tod).^2
powe_d = real.(fft_tod_d).^2 .+ imag.(fft_tod_d).^2
powe_filtered = real.(fft_filtered).^2 .+ imag.(fft_filtered).^2

# Ricordarsi di moltiplicare per due!
plot(freq[1:100:Int(round(length(tod)/2))], 2*powe[1:100:Int(round(length(tod)/2))]./length(tod), xlims=(0,0.5), ylims=(1e-6,1e10), yscale=:log, label="raw tod")
plot!(freq[1:100:Int(round(length(tod)/2))], noise_ss[1:100:Int(round(length(tod)/2))], xlims=(0,0.5), ylims=(1e-6,1e10), yscale=:log, label="noise model slope=3")
plot!(freq[1:100:Int(round(length(tod)/2))], noise_ss2[1:100:Int(round(length(tod)/2))], xlims=(0,0.5), ylims=(1e-6,1e10), yscale=:log, label="noise model slope=2")
plot!(freq[1:100:Int(round(length(tod)/2))], 2*powe_d[1:100:Int(round(length(tod)/2))]./length(tod), xlims=(0,0.5), yscale=:log, ylims=(1e-6,1e10), label="Demod")
plot!(freq[1:100:Int(round(length(tod)/2))], 2*powe_filtered[1:100:Int(round(length(tod)/2))]./length(tod), xlims=(0,0.5), yscale=:log, label="bare frequency filter")
# plot(pp.freq*1/0.0083) #Attenzione a DSP!!! Fa cose che non capisco bene!

plot!(freq[1:100:Int(round(length(tod)/2))], noise_ss[1:100:Int(round(length(tod)/2))], xlims=(0,0.5), ylims=(1e-6,1e10), yscale=:log, label="test")


text:
module PkgRegistryAnalysis

using DataFrames
using OrderedCollections
using Pkg.TOML: parsefile
using AbstractTrees

#-----------------------------------------------------------------------------# tree stuff
struct Registry
    path::String
    Registry(path = joinpath(homedir(), ".julia", "registries", "General")) = new(path)
end

struct Letter
    path::String
end

struct Package
    path::String
end

AbstractTrees.children(r::Registry) = Letter.(joinpath.([r.path], string.('A':'Z')))
AbstractTrees.children(l::Letter) = Package.(joinpath.([l.path], readdir(l.path)))
AbstractTrees.children(p::Package) = ()

#-----------------------------------------------------------------------------# package_data
function package_data(registrypath = joinpath(homedir(), ".julia", "registries", "General"))
    DataFrame(retrieve_package_data(row.path) for row in AbstractTrees.Leaves(Registry(registrypath)))
end

function retrieve_package_data(path::String)
    package = parsefile(joinpath(path, "Package.toml"))
    versions = sort(parsefile(joinpath(path, "Versions.toml")))
    return (
        name = package["name"],
        repo = package["repo"],
        uuid = package["uuid"],
        latest = last(collect(keys(versions))),
        versions = Dict(k => v["git-tree-sha1"] for (k,v) in pairs(versions)),
        deps = isfile(joinpath(path, "Deps.toml")) ? sort(parsefile(joinpath(path, "Deps.toml"))) : nothing,
        compat = isfile(joinpath(path, "Compat.toml")) ? sort(parsefile(joinpath(path, "Compat.toml"))) : nothing
    )
end

function dependents(df::DataFrame, pkg::String)
    inds = findall(x -> !isnothing(x) && any(y -> pkg in keys(y), values(x)), df.deps)
    df[inds, :]
end

end # module


text:
export
    @op_perf,
    @make_c,
    runtime

# All TODO

macro op_perf(op_perf_body)
    return esc(op_perf_macro(op_perf_body))
end

function op_perf_macro(op_perf_body::Expr)
    return quote
        $op_perf_body
        export $(op_perf_body.args[1].args[1])
    end
end

macro make_c()
    return (:1)
end

function runtime end