text:
# Integer conic problems

function intsoc1test(model::MOI.ModelLike, config::Config)
    atol = config.atol
    rtol = config.rtol
    # Problem SINTSOC1
    # min 0x - 2y - 1z
    #  st  x            == 1
    #      x >= ||(y,z)||
    #      (y,z) binary
    @test MOI.supports_incremental_interface(model)
    @test MOI.supports(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
    )
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{Float64},
        MOI.Zeros,
    )
    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.ZeroOne)
    @test MOI.supports_constraint(
        model,
        MOI.VectorOfVariables,
        MOI.SecondOrderCone,
    )
    MOI.empty!(model)
    @test MOI.is_empty(model)
    x, y, z = MOI.add_variables(model, 3)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([-2.0, -1.0], [y, z]),
            0.0,
        ),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    ceq = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x))],
            [-1.0],
        ),
        MOI.Zeros(1),
    )
    csoc = MOI.add_constraint(
        model,
        MOI.VectorOfVariables([x, y, z]),
        MOI.SecondOrderCone(3),
    )
    if config.query_number_of_constraints
        @test MOI.get(
            model,
            MOI.NumberOfConstraints{
                MOI.VectorAffineFunction{Float64},
                MOI.Zeros,
            }(),
        ) == 1
        @test MOI.get(
            model,
            MOI.NumberOfConstraints{MOI.VectorOfVariables,MOI.SecondOrderCone}(),
        ) == 1
    end
    loc = MOI.get(model, MOI.ListOfConstraintTypesPresent())
    @test length(loc) == 2
    @test (MOI.VectorAffineFunction{Float64}, MOI.Zeros) in loc
    @test (MOI.VectorOfVariables, MOI.SecondOrderCone) in loc
    bin1 = MOI.add_constraint(model, y, MOI.ZeroOne())
    # We test this after the creation of every `VariableIndex` constraint
    # to ensure a good coverage of corner cases.
    @test bin1.value == y.value
    bin2 = MOI.add_constraint(model, z, MOI.ZeroOne())
    @test bin2.value == z.value
    if config.solve
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ -2 atol = atol rtol = rtol
        @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 1 atol = atol rtol =
            rtol
        @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 1 atol = atol rtol =
            rtol
        @test MOI.get(model, MOI.VariablePrimal(), z) ≈ 0 atol = atol rtol =
            rtol
    end
end

const intsoctests = Dict("intsoc1" => intsoc1test)

@moitestset intsoc

const intconictests = Dict("intsoc" => intsoctest)

@moitestset intconic true


text:
using StackOverflow

"""
    getrecentbadgerecipients()

Returns recently awarded badges in the system.
"""
function getrecentbadgerecipients()
    r = HTTP.request("GET", "https://api.stackexchange.com/2.2/badges/recipients?site=stackoverflow")
    #TODO: add fromdate and todate.

    json_obj = StackOverflow.convert_HTTP_Response_To_JSON(r)

    for (k, v) in json_obj
        if occursin("items", k)
            for item in v
                for (key, value) in item
                    println("$(key) $(value)")
                end
                println("\n")
            end
        end
    end

    #TODO: May want to return the data in an array so it can be used.
    # The data can just be printed and/or pushed to a array then returned.

end

ext: jl


text:
using Dolang: sanitize
using Lerche: Tree, Token

import Dolo
using Dolo: ModelCalibration
import Dolang: solve_triangular_system

import YAML
import Dolang
import Dolang: yaml_node_from_string, yaml_node_from_file
using Dolang: SymExpr
import Dolang: LTree
import Dolang

using DataStructures: OrderedDict

# defind language understood in yaml files
language = minilang

abstract type AbstractModel{ID} end

mutable struct Model{ID} <: AbstractModel{ID}

    data::YAML.MappingNode
    symbols::Dict{Symbol, Vector{Symbol}}
    calibration::ModelCalibration
    exogenous
    domain::AbstractDomain
    factories
    definitions

    function Model{ID}(data::YAML.Node) where ID
        model = new{ID}(data)
        model.calibration = get_calibration(model)
        model.symbols = get_symbols(model)
        model.exogenous = get_exogenous(model)
        model.domain = get_domain(model)

        factories = get_factories(model)
        model.factories = factories
        for (k,fact) in factories
            code = Dolang.gen_generated_gufun(fact; dispatch=typeof(model))
            # print_code && println("equation '", eq_type, "'", code)
            Core.eval(Dolo, code)
        end

        # # Create definitions
        defs = get_definitions(model; stringify=true)
        model.definitions = defs
        definitions = OrderedDict{Symbol, SymExpr}( [(stringify(k),v) for (k,v) in defs])
        vars = cat(model.symbols[:exogenous], model.symbols[:states], model.symbols[:controls]; dims=1)
        args = OrderedDict(
            :past => [Dolang.stringify(v,-1) for v in vars],
            :present => [Dolang.stringify(v,0) for v in vars],
            :future => [Dolang.stringify(v,1) for v in vars],
            :params => [Dolang.stringify(e) for e in model.symbols[:parameters]]
        )
        ff = Dolang.FunctionFactory(definitions, args, OrderedDict{Symbol, SymExpr}(), :definitions)
        code = Dolang.gen_generated_gufun(ff;dispatch=typeof(model), funname=:evaluate_definitions)
        Core.eval(Dolo,code)
        model.factories[:definitions] = ff

        return model
    end

end

function Model(url::AbstractString; print_code=false)

    # it looks like it would be cool to use the super constructor ;-)
    if match(r"(http|https):.*", url) != nothing
        res = HTTP.request("GET", url)
        txt = String(res.body)
    else
        txt = read(open(url), String)
    end
    data = Dolang.yaml_node_from_string(txt)
    id = gensym()
    fname = basename(url)
    return Model{id}(data)

end

yaml_import(filename::AbstractString) = Model(filename)

function Base.show(io::IO, model::Model)
    print(io, "Model")
end

function get_symbols(model::Model)
    syms = model.data[:symbols]
    symbols = OrderedDict{Symbol,Vector{Symbol}}()
    for k in keys(syms)
        symbols[Symbol(k)] = [Symbol(e.value) for e in syms[k]]
    end
    return symbols
end

function get_variables(model::Model)
    symbols = get_symbols(model)
    vars = cat(values(symbols)...; dims=1)
    dynvars = setdiff(vars, symbols[:parameters] )
    dynvars = union( dynvars, get_defined_variables(model) )
end

function get_defined_variables(model::Model)
    if !("definitions" in keys(model.data))
        return Symbol[]
    else
        return [Symbol(e) for e in keys(model.data["definitions"])]
    end
end

function get_name(model::Model)
    if "name" in keys(model.data)
        name = model.data[:name].value
    else
        name = "anonymous"
    end
    return name
end

function get_exogenous(model::AModel)

    rdata = model.data
    if !("exogenous" in keys(model.data))
        return nothing
    end
    exo_dict = model.data[:exogenous]
    cond = !(exo_dict.tag=="tag:yaml.org,2002:map")

    syms = cat( [[Symbol(strip(e)) for e in split(k, ",")] for k in keys(exo_dict)]..., dims=1)
    expected = model.symbols[:exogenous]
    if (syms != expected)
        msg = string("In 'exogenous' section, shocks must be declared in the same order as shock symbol. Found: ", syms, ". Expected: ", expected, ".")
        throw(ErrorException(msg))
    end
    calibration = model.calibration.flat
    processes = []
    for k in keys(exo_dict)
        v = exo_dict[k]
        p = Dolang.eval_node(v, calibration, minilang, FromGreek())
        push!(processes, p)
    end
    return ProductProcess(processes...)

end

function get_calibration(model::Model)
    calib = model.data[:calibration]
    # prep calib: parse to Expr, Symbol, or Number
    _calib  = OrderedDict{Symbol,Union{Expr,Symbol,Number}}()
    # add calibration for definitions
    for k in keys(calib)
        x = Dolang.parse_equation(calib[k].value)
        e = Dolang.convert(Expr, x, stringify=false)
        _calib[Symbol(k)] = e
    end
    # so far _calib is a symbolic calibration
    calibration = solve_triangular_system(_calib)
    symbols = get_symbols(model)
    return ModelCalibration(calibration, symbols)
end

function set_calibration!(model::Model, key::Symbol, value)
    # TODO: set proper type for ScalarNode
    # this will fail is parameter wasn't defined before
    model.data[:calibration][key].value = string(value)
    model.calibration = get_calibration(model)
    model.exogenous = get_exogenous(model)
    model.domain = get_domain(model)
end

function set_calibration!(model::Model, values::AbstractDict{Symbol, Any})
    calib = model.data[:calibration]
    for (k,v) in values
        calib[k].value = string(v)
    end
    model.calibration = get_calibration(model)
    model.exogenous = get_exogenous(model)
    model.domain = get_domain(model)
end

function set_calibration!(model::Model; kwargs...)
    set_calibration!(model, Dict(kwargs))
end

function get_equation_block(model, eqname; stringify=true)
    variables = get_variables(model)
    variables_str = String[string(e) for e in variables]

    yml_node = model.data["equations"][eqname]

    if "tag:yaml.org,2002:str" == yml_node.tag
        block = Dolang.parse_equation_block(yml_node.value; variables=variables_str)
        equations = block.children
    else
        equations = [Dolang.parse_equation(line.value; variables=variables_str) for line in yml_node]
    end

    eqs = Union{Expr, Number, Symbol}[]
    eqs_lb = Union{Expr, Number, Symbol}[]
    eqs_ub = Union{Expr, Number, Symbol}[]

    for eq in equations
        if eq.data == "double_complementarity"
            inequality = eq.children[2]
            # @assert inegality.data == "double_inequality"
            eq_lb = Dolang.convert(Expr, inequality.children[1]; stringify=stringify)
            eq_cc = Dolang.convert(Expr, inequality.children[2]; stringify=stringify)
            # TODO: check complementarities order
            eq_ub = Dolang.convert(Expr, inequality.children[3]; stringify=stringify)
            eq = eq.children[1]
            push!(eqs_lb, eq_lb)
            push!(eqs_ub, eq_ub)
        else
            push!(eqs_lb, :(-Inf) )
            push!(eqs_ub, :(Inf) )
        end
        if eq.data == "equality"
            eq = Tree("sub", [eq.children[2], eq.children[1]])
        end
        # eq = Dolang.stringify(eq)
        push!(eqs, Dolang.convert(Expr, eq; stringify=stringify))
    end
    return (eqs, eqs_lb, eqs_ub)
end

function get_assignment_block(model, eqname; stringify=true)
    variables = get_variables(model)
    variables_str = String[string(e) for e in variables]

    yml_node = model.data["equations"][eqname]

    if "tag:yaml.org,2002:str" == yml_node.tag
        block = Dolang.parse_assignment_block(yml_node.value; variables=variables_str)
        equations = block.children
    else
        equations = [Dolang.parse_assignment(line.value; variables=variables_str) for line in yml_node]
    end

    # eqs = OrderedDict{Tuple{Symbol, Int64}, SymExpr}()
    eqs = OrderedDict()
    for eq in equations

        lhs = eq.children[1]
        rhs = eq.children[2]
        dest_s = Symbol(lhs.children[1].children[1].value)
        dest_d = parse(Int, lhs.children[2].children[1].value)

        dest = Dolang.stringify(dest_s, dest_d)

        # TODO : remove type instability
        if stringify
            dest = Dolang.stringify(dest_s, dest_d)
        else
            dest = (dest_s, dest_d)
        end

        expr = Dolang.convert(Expr, rhs; stringify=stringify)

        eqs[dest] = expr

    end
    return eqs
end

function get_domain(model::Model)::AbstractDomain

    states = model.symbols[:states]

    if !("domain" in keys(model.data))
        return EmptyDomain(states)
    end

    domain = model.data["domain"]

    calib = model.calibration.flat

    # TODO deal with empty dict and make robust construction
    min = [Dolang.eval_node(domain[(k)][1], calib) for k in states]
    max = [Dolang.eval_node(domain[(k)][2], calib) for k in states]
    return Domain(states, min, max)

end

function get_definitions(model::Model; tshift=0, stringify=false) # ::OrderedDict{Tuple{Symbol,Int}}
    # parse defs so values are Expr
    if "definitions" in keys(model.data)
        defs = model.data[:definitions]
    else
        defs = Dict()
    end
    dynvars = string.( cat(get_variables(model), keys(defs)...; dims=1) )

    _defs =  OrderedDict{Tuple{Symbol,Int64},Union{Expr, Number, Symbol}}()  # {Symbol,Int}()

    for (kkk,v) in defs
        kk = Dolang.parse_equation(kkk)
        if kk.data == "symbol"
            k = kk.children[1].value
        elseif kk.data=="variable"
            k = kk.children[1].children[1].value
        else
            error("Invalid key.")
        end
        vv = Dolang.parse_equation(v.value; variables=(dynvars))::LTree
        if tshift != 0
            vv = Dolang.time_shift(vv, tshift)::LTree
        end
        _defs[(Symbol(k),tshift)] = Dolang.convert(Expr, vv; stringify=stringify)::SymExpr
    end

    return _defs
end

function get_options(model::Model)
    cc = model.calibration.flat
    if "options" in keys(model.data)
        return Dolang.eval_node(model.data["options"], cc, language)
    else
        return Dict()
    end
end

import Dolang: FunctionFactory
import Dolang: stringify

function get_factories(model::Model)

    facts = Dict()
    for eq_type in ["transition", "arbitrage"]
        if eq_type == "transition"
            factories = (get_factory(model, eq_type), )
        else
            factories = get_factory(model, eq_type)
        end
        for f in factories
            facts[f.funname] = f
        end
    end

    return facts

end

function get_factory(model::Model, eq_type::String)
    if eq_type == "arbitrage"
        defs_0 = get_definitions(model; stringify=true)
        defs_1 = get_definitions(model; tshift=1, stringify=true)
        definitions = OrderedDict{Symbol, SymExpr}([ (Dolang.stringify(k), v) for (k,v) in merge(defs_0, defs_1)])

        eqs, eq_lb, eq_ub = get_equation_block(model, eq_type)

        symbols = get_symbols(model)

        #TODO : fix crazy bug: it doesn't work without the trailing underscore !
        equations = OrderedDict{Symbol, Union{Expr, Number, Symbol}}(
            Symbol(string("out_", i, "_")) => eqs[i]
            for i =1:length(eqs)
        )
        arguments = OrderedDict(
            :m => [stringify(e,0) for e in symbols[:exogenous]],
            :s => [stringify(e,0) for e in symbols[:states]],
            :x => [stringify(e,0) for e in symbols[:controls]],
            :M => [stringify(e,1) for e in symbols[:exogenous]],
            :S => [stringify(e,1) for e in symbols[:states]],
            :X => [stringify(e,1) for e in symbols[:controls]],
            :p => [stringify(e) for e in symbols[:parameters]]
            )

        ff = FunctionFactory(equations, arguments, definitions, Symbol(eq_type))

        equations_lb = OrderedDict{Symbol, Union{Expr, Number, Symbol}}(
                Symbol(string("out_", i, "_")) => eq_lb[i]
                for i =1:length(eqs)
            )
        equations_ub = OrderedDict{Symbol, Union{Expr, Number, Symbol}}(
                Symbol(string("out_", i,"_")) => eq_ub[i]
                for i =1:length(eqs)
            )
        arguments_bounds = OrderedDict(
            :m => [stringify(e,0) for e in symbols[:exogenous]],
            :s => [stringify(e,0) for e in symbols[:states]],
            :p => [stringify(e) for e in symbols[:parameters]]
            )

        # definitions: we should remove definitions depending on controls
        # definitions = OrderedDict{Symbol, SymExpr}([ (Dolang.stringify(k), v) for (k,v) in defs_0] )
        definitions = Dict{Symbol, SymExpr}()
        ff_lb = FunctionFactory(equations_lb, arguments_bounds, definitions, :controls_lb)
        ff_ub = FunctionFactory(equations_ub, arguments_bounds, definitions, :controls_ub)

        return (ff, ff_lb, ff_ub)

    else
        # defs_0  = get_definitions(model; stringify=true)
        defs_m1 = get_definitions(model; tshift=-1, stringify=true)

        definitions = OrderedDict{Symbol, SymExpr}([ (Dolang.stringify(k), v) for (k,v) in  defs_m1 ])

        equations = get_assignment_block(model, eq_type)
        symbols = get_symbols(model)
        arguments = OrderedDict(
            :m => [stringify(e,-1) for e in symbols[:exogenous]],
            :s => [stringify(e,-1) for e in symbols[:states]],
            :x => [stringify(e,-1) for e in symbols[:controls]],
            :M => [stringify(e,0) for e in symbols[:exogenous]],
            :p => [stringify(e) for e in symbols[:parameters]]
        )
        ff = FunctionFactory(equations, arguments, definitions, Symbol(eq_type))
        return ff
    end

end

function get_options(model::AModel; options=Dict())
    if "options" in keys(model.data)
        return model.data["options"]
        # opts = model.data["options"]
        # return rmerge(opts, options)
    else
        return Dict()
    end
end

using Dolo: CartesianGrid

function get_grid(model::AModel; options=Dict())

    domain = get_domain(model)
    d = length(domain.states)

    c = model.calibration.flat

    if ("options" in keys(model.data)) && ( "grid" in keys(model.data["options"]) )

        gopt = model.data["options"]["grid"]
        if gopt.tag == "!Cartesian"
            if "n" in keys(gopt)
                orders = Dolang.eval_node(gopt[:n], c)
            elseif "orders" in keys(gopt)
                orders = Dolang.eval_node(gopt[:orders], c)
            else
                orders = [20 for i=1:d]
            end
            grid = CartesianGrid{d}(domain.min, domain.max, orders)
            if length(orders)!=length(model.calibration[:states])
                msg = string("Check the dimension of the matrix given in the yaml file, section: options-grid-orders. ",
                             "Expected to be of dimension $([length(model.calibration[:states])])")
                error(msg)
            end
        elseif gopt.tag == "Smolyak"
            mu = gopt[:mu]
            grid = SmolyakGrid{d}(domain.min, domain.max, mu)
        elseif gopt.tag == "Random"
            n = gopt[:N]
            grid = RandomGrid{d}(domain.min, domain.max, n)
        else
            error("Unknown grid type.")
        end

    else

        # default grid
        n = [20 for i=1:d]
        grid = CartesianGrid{d}(domain.min, domain.max, n)

    end

    return grid
end


text:
function PrintDivider()
println("____________________________________________________________")
end

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
"Plot the response across time"
@recipe function plot(simulation::TimeSimulation)
    xguide --> "Time (t)"
    yguide --> "Response"
    grid --> false
    title --> "Response from particles of radius $(signif(simulation.frequency_simulation.particles[1].r,2)) contained in a $(lowercase(name(simulation.frequency_simulation.shape)))\n with volfrac=$(signif(calculate_volfrac(simulation.frequency_simulation),2)) measured at ($(simulation.frequency_simulation.listener_positions[1,1]), $(simulation.frequency_simulation.listener_positions[2,1]))"

    (simulation.time_arr, simulation.response)
end

# "Plot the field for an array of time"
# @recipe function plot{T}(TimeSimulation::TimeSimulation{T}; res=res)
#     map(axes(TimeSimulation.time_arr,1)) do i
#         simulation = deepcopy(TimeSimulation)
#         simulation.response = reshape(TimeSimulation.response[i,:],1,:)
#         plot(simulation,TimeSimulation.time_arr[i]; res=res, build_field=false)
#     end
# end

"Plot the field for one time"
@recipe function plot{T}(TimeSimulation::TimeSimulation{T}, t::Union{T,Vector{T}};
                        res=20, xres=res, yres=res, resp_fnc=real,
                        drawshape = false, build_field=true, drawlisteners = build_field)
    simulation = TimeSimulation.frequency_simulation

    if isa(t,T) t_arr = [t]
    else        t_arr = t  end

    @series begin
        # find a box which covers everything
        shape_bounds = bounding_box(simulation.shape)
        listeners_as_particles = map(
            l -> Particle(simulation.listener_positions[:,l],mean_radius(simulation)/2),
            1:size(simulation.listener_positions,2)
        )
        particle_bounds = bounding_box([simulation.particles; listeners_as_particles])
        bounds = bounding_box(shape_bounds, particle_bounds)

        if build_field
          field_simulation = run(simulation, bounds; xres=xres, yres=yres)
          field_TimeSimulation = deepcopy(TimeSimulation) # to use all the same options/fields as TimeSimulation
          field_TimeSimulation.frequency_simulation = field_simulation
          generate_responses!(field_TimeSimulation, t_arr)
        else
          field_TimeSimulation = TimeSimulation
        end

        # For this we sample at the centre of each pixel
        x_pixels = LinRange(bounds.bottomleft[1], bounds.topright[1], xres+1)
        y_pixels = LinRange(bounds.bottomleft[2], bounds.topright[2], yres+1)

        # NOTE only plots the first time plot for now...
        response_mat = transpose(reshape(field_TimeSimulation.response[1,:], (xres+1, yres+1)))
        linetype --> :contour
        fill --> true
        fillcolor --> :balance
        title --> "Field at time=$(t_arr[1])"

        (x_pixels, y_pixels, resp_fnc.(response_mat))
    end
    if drawshape
      @series begin
          simulation.shape
      end
    end
    for i=1:length(simulation.particles) @series simulation.particles[i] end
    if drawlisteners
      @series begin
          line --> 0
          fill --> (0, :lightgreen)
          grid --> false
          colorbar --> true
          aspect_ratio --> 1.0
          legend --> false

          r = mean_radius(simulation.particles)/2
          x(t) = r * cos(t) + simulation.listener_positions[1, 1]
          y(t) = r * sin(t) + simulation.listener_positions[2, 1]

          (x, y, -2π/3, 2π/3)
      end
    end

    # clims --> (-0.65,0.65) # so that white is always = 0, these values are suitable when using default gaussian impulse
end


text:
# Node
abstract AbstractNode{T<:Number,N}
abstract AbstractFourierNode{T<:Number,N} <: AbstractNode{T,N}

AbstractFourierNode{T<:FFTW.fftwReal}(node::AbstractNode{T}, args...) =
    RealFourierNode(node, args...)
AbstractFourierNode{T<:FFTW.fftwComplex}(node::AbstractNode{T}, args...) =
    ComplexFourierNode(node, args...)

immutable RealFourierNode{T<:FFTW.fftwComplex,R<:FFTW.fftwReal,N} <:
        AbstractFourierNode{T,N}
    data::Array{T,N}
    forwardplan::FFTW.rFFTWPlan{R,-1,false,N}
    ranges::NTuple{N,Pair{PathKey,StepRange{Int,Int}}}
end

function RealFourierNode{T<:FFTW.fftwReal}(
        node::AbstractNode{T},
        region::AbstractArray{Int,1},
        flags::UInt32 = FFTW.ESTIMATE,
        timelimit = Inf)
    forwardplan = plan_rfft(
            node.data, region, flags = flags, timelimit = timelimit)
    RealFourierNode(forwardplan * node.data, forwardplan, node.ranges)
end

immutable ComplexFourierNode{T<:FFTW.fftwComplex,N} <:
        AbstractFourierNode{T,N}
    data::Array{T,N}
    forwardplan::FFTW.cFFTWPlan{T,-1,false,N}
    ranges::NTuple{N,Pair{PathKey,StepRange{Int,Int}}}
end

function ComplexFourierNode{T<:FFTW.fftwComplex}(
        node::AbstractNode{T},
        region::AbstractArray{Int,1} ;
        flags = FFTW.ESTIMATE,
        timelimit = Inf)
    forwardplan =
        plan_fft(node.data, region ; flags = flags, timelimit = timelimit)
    ComplexFourierNode(forwardplan * node.data, forwardplan, node.ranges)
end

immutable Node{T<:Number,N} <: AbstractNode{T,N}
    data::Array{T,N}
    ranges::NTuple{N,Pair{PathKey,StepRange{Int,Int}}}
end

immutable InvComplexFourierNode{T<:FFTW.fftwComplex,N,R<:FFTW.fftwReal} <:
        AbstractNode{T,N}
    data::Array{T,N}
    inverseplan::Base.DFT.ScaledPlan{T,FFTW.cFFTWPlan{T,1,false,N},R}
    ranges::NTuple{N,Pair{PathKey,StepRange{Int,Int}}}
end

function InvComplexFourierNode{T<:FFTW.fftwComplex}(
        node::AbstractNode{T},
        region::Vector{Int},
        flags::UInt32 = FFTW.ESTIMATE,
        timelimit = Inf)
    inverseplan =
        plan_ifft(node.data, region ; flags = flags, timelimit = timelimit)
    InvComplexFourierNode(inverseplan * node.data, inverseplan, node.ranges)
end

function pathdepth(node::AbstractNode, refkey::PathKey)
    mapreduce(p -> pathdepth(p.first, refkey), max, 1, node.ranges)
end

function transform!{T}(
        destination::SubArray,
        ψ::AbstractFilter{T,FourierDomain{1}},
        node::AbstractNode,
        dims::Array{Int})
    transform!(destination, ψ, node, dims[1])
end

function transform!(
        destination::SubArray,
        ψ::Analytic1DFilter,
        node::AbstractNode,
        dim::Int)
    inds = Union{Colon,Int}[
        fill(Colon(), dim-1) ; 0 ; fill(Colon(), ndims(destination)-1)]
    @inbounds for ω in (ψ.posfirst):(ψ.posfirst+length(ψ.pos)-1)
        inds[dim] = 1 + ω
        input = view(node.data, inds...)
        output = view(destination, inds...)
        broadcast!(*, output, ψ.pos[1 - ψ.posfirst + ω], input)
    end
end

function transform!(
        destination::SubArray,
        ψ::FullResolution1DFilter,
        node::RealFourierNode,
        dim::Int)
    inds = Union{Colon,Int}[
        fill(Colon(), dim-1) ; 0 ; fill(Colon(), ndims(destination)-1)]
    halfN = size(node.data, dim) - 1
    # positive frequencies (including 0 and π)
    @inbounds for ω in 0:halfN
        inds[dim] = 1 + ω
        input = view(node.data, inds...)
        output = view(destination, inds...)
        broadcast!(*, output, ψ.coeff[1 + ω], input)
    end
    # negative frequencies
    @inbounds for ω in (halfN+1):(2*halfN-1)
        inds[dim] = 1 + 2*halfN - ω
        input = view(node.data, inds...)
        inds[dim] = 1 + ω
        output = view(destination, inds...)
        broadcast!(A_mul_Bc, output, ψ.coeff[1 + ω], input)
    end
end

function transform!(
        destination::SubArray,
        ψ::FullResolution1DFilter,
        node::ComplexFourierNode,
        dim::Int)
    inds = Union{Colon,Int}[
        fill(Colon(), dim-1) ; 0 ; fill(Colon(), ndims(destination)-1)]
    @inbounds for ω in 0:(size(node.data, dim)-1)
        inds[dim] = 1 + ω
        input = view(node.data, inds...)
        output = view(destination, inds...)
        broadcast!(*, output, ψ.coeff[1 + ω], input)
    end
end

function transform!(
        destination::SubArray,
        ψ::Vanishing1DFilter,
        node::ComplexFourierNode,
        dim::Int)
    inds = Union{Colon,Int}[
        fill(Colon(), dim-1) ; 0 ; fill(Colon(), ndims(destination)-1)]
    @inbounds for ω in ψ.an.posfirst+(0:(length(ψ.an.pos)-1))
        inds[dim] = 1 + ω
        input = view(node.data, inds...)
        output = view(destination, inds...)
        broadcast!(*, output, ψ.an.pos[1 - ψ.an.posfirst + ω], input)
    end
    @inbounds for ω in ψ.coan.neglast+(0:-1:(1-length(ψ.coan.neg)))
        inds[dim] = 1 + size(node.data, dim) + ω
        input = view(node.data, inds...)
        inds[dim] = 1 + size(destination, dim) + ω
        output = view(destination, inds...)
        broadcast!(*, output, ψ.coan.neg[1 + end + ω], input)
    end
end

function transform!(
        destination::SubArray,
        ψ::Vanishing1DFilter,
        node::RealFourierNode,
        dim::Int)
    inds = Union{Colon,Int}[
        fill(Colon(), dim-1) ; 0 ; fill(Colon(), ndims(destination)-1)]
    @inbounds for ω in ψ.an.posfirst+(0:(length(ψ.an.pos)-1))
        inds[dim] = 1 + ω
        input = view(node.data, inds...)
        output = view(destination, inds...)
        broadcast!(*, output, ψ.an.pos[1 - ψ.an.posfirst + ω], input)
    end
    inds_in, inds_out = inds, Base.copy(inds)
    @inbounds for ω in ψ.coan.neglast+(0:-1:(1-length(ψ.coan.neg)))
        inds_in[dim] = 1 - ω
        input = view(node.data, inds_in...)
        inds_out[dim] = 1 + size(destination, dim) + ω
        output = view(destination, inds_out...)
        broadcast!(A_mul_Bc, output,
            ψ.coan.neg[end - ψ.coan.neglast + ω], input)
    end
end

function transform!(
        destination::SubArray,
        ψ::VanishingWithMidpoint1DFilter,
        node::RealFourierNode,
        dim::Int)
    inds = Union{Colon,Int}[
        fill(Colon(), dim-1) ; 0 ; fill(Colon(), ndims(destination)-1)]
    @inbounds for ω in ψ.an.posfirst+(0:(length(ψ.an.pos)-1))
        inds[dim] = 1 + ω
        input = view(node.data, inds...)
        output = view(destination, inds...)
        broadcast!(*, output, ψ.an.pos[1 - ψ.an.posfirst + ω], input)
    end
    @inbounds for ω in ψ.coan.neglast+(0:-1:(1-length(ψ.coan.neg)))
        inds[dim] = 1 - ω
        input = view(node.data, inds...)
        inds[dim] = 1 + size(destination, dim) + ω
        output = view(destination, inds...)
        broadcast!(A_mul_Bc, output,
            ψ.coan.neg[end - ψ.coan.neglast + ω], input)
    end
    @inbounds begin
        ωmidpoint = ψ.an.posfirst + length(ψ.an.pos)
        inds[dim] = 1 + ωmidpoint
        input = view(node.data, inds...)
        output = view(destination, inds...)
        broadcast!(*, output, ψ.midpoint, input)
    end
end


text:
module PeriodicBCsTestsSeq
using PartitionedArrays
using Gridap
using GridapDistributed
using Test

nps = [(1,1),(1,2),(2,1),(2,2)]
domain = (0,1,0,1)
cells = (4,4)
isperiodics = [(false,false),(false,true),(true,false),(true,true)]
reffe = ReferenceFE(lagrangian,Float64,1)
ndofss = [25,20,20,16]
for np in nps
  for (ndofs, isperiodic) in zip(ndofss,isperiodics)
    parts = get_part_ids(sequential,np)
    model = CartesianDiscreteModel(parts,domain,cells;isperiodic=isperiodic)
    V = FESpace(model,reffe)
    @test ndofs == num_free_dofs(V)
  end
end

include("../PeriodicBCsTests.jl")
prun(PeriodicBCsTests.main,sequential,(2,2))
prun(PeriodicBCsTests.main,sequential,(2,1))
prun(PeriodicBCsTests.main,sequential,(1,1))
prun(PeriodicBCsTests.main,sequential,(1,2))
prun(PeriodicBCsTests.main,sequential,(2,3))
end # module


text:
#Simple testing of new TCP comms
using Base.Test

function write_memory(sock::TCPSocket, addr::UInt32, data::Vector{UInt32})
	#Have 16 bits for count in datagram so can write up to (1 << 16) - 1 words
	#but need 128bit alignment for SDRAM so 0xfffc
	idx = 1
	max_ct = 0xfffc
	datagrams_written = 0
	init_addr = addr
	while (idx < length(data))
		ct_left = length(data) - idx + 1
		ct = (ct_left < max_ct) ? ct_left : max_ct
		cmd = 0x80000000
		datagram = UInt32[cmd + ct; addr; data[idx:idx+ct-1]]
		write(sock, map(hton, datagram))
		datagrams_written += 1
		idx += ct
		addr += ct*4
	end

	#Check the results
	results = map(ntoh,read(tcp_sock, UInt32, 2*datagrams_written))
	addr = init_addr
	for ct = 1:datagrams_written
		words_written = (ct == datagrams_written) ? length(data)-((datagrams_written-1)*0xfffc) : 0xfffc
		@test results[2*ct-1] == 0x80800000 + UInt32(words_written)
		@test results[2*ct] == UInt32(addr)
		addr += 4*words_written
	end

end

ip_addr = ip"192.168.2.2"

#UDP responder
udp_sock = UDPSocket();
bind(udp_sock, ip"0.0.0.0", 0xbb4f)
send(udp_sock, ip"192.168.2.2", 0xbb4f, [0x01])
resp = bytestring(recv(udp_sock))
@test resp == "I am an APS2"

#TCP reset
# send(udp_sock, ip"192.168.2.2", 0xbb4f, [0x02])

close(udp_sock)

tcp_sock = connect(ip_addr, 0xbb4e)

#Read uptime
datagram = UInt32[0x10000002, 0x44a00050]
write(tcp_sock, map(hton, datagram))
#check response datagram header
resp_header = map(ntoh, read(tcp_sock, UInt32, 2))
@test resp_header == UInt32[0x10000002; 0x44a00050]

uptime_array = map(ntoh,read(tcp_sock, UInt32, 2))
uptime = uptime_array[1] + 1e-9*uptime_array[2]
println("Uptime $uptime")

#memory write
# first BRAMs with short message
write_vals = rand(UInt32, 1 << 10)
write_memory(tcp_sock, 0xc2000000, write_vals)

#Send a read request
datagram = UInt32[0x100003e8, 0xc2000000]
write(tcp_sock, map(hton, datagram))
#check response datagram header
resp_header = map(ntoh, read(tcp_sock, UInt32, 2))
@test resp_header == UInt32[0x100003e8; 0xc2000000]

read_vals = map(ntoh, read(tcp_sock, UInt32, 1000))
@test read_vals == write_vals[1:1000]

#Now SDRAM
write_vals = rand(UInt32, 1 << 18)
tic()
write_memory(tcp_sock, 0x20000000, write_vals)
write_speed = sizeof(eltype(write_vals)) * length(write_vals) / toq()
println("Write speed: $(write_speed/1e6) MB/s")

idx = 0
while ( (idx+1000) < length(write_vals))	#Send a read request
	datagram = UInt32[0x100003e8, 0x20000000+(idx*4)]
	write(tcp_sock, map(hton, datagram))
	#check response datagram header
	resp_header = map(ntoh, read(tcp_sock, UInt32, 2))
	@test resp_header == UInt32[0x100003e8; 0x20000000+(idx*4)]
	read_vals = map(ntoh, read(tcp_sock, UInt32, 1000))
	@test read_vals == write_vals[idx+1:idx+1000]
	idx += 1000
end

#Test status register request to wrapped ApsMsgProc
read_req = UInt32[0x37000010, 0x00000000]
write(tcp_sock, map(hton, read_req))
result = map(ntoh, read(tcp_sock, UInt32, 17))
if result[4] == 0xdddddddd
	@test result[1:5] == UInt32[0x97000010; 0x00000a20 ; 0xbadda555; 0xdddddddd; 0x0ddba111]
else
	@test result[1:5] == UInt32[0x97000010; 0x00000a20 ; 0xbadda555; 0xeeeeeeee; 0x0ddba111]
end

#Test configuration DRAM read/write
num_words = 1 << 8
write_vals = rand(UInt32, num_words)
cmd = 0x25000000
datagram = UInt32[cmd+num_words; 0; write_vals]
write(tcp_sock, map(hton, datagram))
write_resp = ntoh(read(tcp_sock, UInt32))
@test write_resp == 0x85000000
read_req = UInt32[0x35000000 + num_words, 0x0]
write(tcp_sock, map(hton, read_req))
resp_header = ntoh(read(tcp_sock, UInt32))
@test resp_header == UInt32(0x95000000 + num_words)
read_vals = map(ntoh, read(tcp_sock, UInt32, num_words))
@test read_vals == write_vals

#
# #Test flash read
# read_req = UInt32[0x32000004, 0x00ff0000]
# write(tcp_sock, map(hton, read_req))
# result = map(ntoh, read(tcp_sock, UInt32, 5))
# @test result == UInt32[0x92000004; 0x4651db00; 0x002e0000; 0xc0a80202; 0x00000000]

close(tcp_sock)


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

name = "Serd"
# <-- this is a lie, we're building v0.30.10, but we need to bump version to build for julia v1.6
version_fake = v"0.30.11"
version = v"0.30.10"

# Collection of sources required to complete build
sources = [
    ArchiveSource("http://download.drobilla.net/serd-$(version).tar.bz2", "affa80deec78921f86335e6fc3f18b80aefecf424f6a5755e9f2fa0eb0710edf")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/serd*/

install_license COPYING
./waf configure --prefix=$prefix
./waf
./waf install
exit
"""

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

# The products that we will ensure are always built
products = [
    LibraryProduct(["libserd-0", "serd"], :libserd)
]

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

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version_fake, sources, script, platforms, products, dependencies, julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
@testset "Extract free parameters" begin
  pmodel = Model.m
  D = @struct_2_dict(pmodel)

  @free D["pmodel.x"] Distributions.Normal
  @free D["pmodel.a.x"]

  fp = extract_free_parameters(D)
  @test fp[1].val ≈ 1.0
  @test fp[2].val ≈ 3.0
  @test fp[2].prior == Normal

end


text:
"""
ReadableNumStyle field naming
    _integral digits_ preceed the fraction_marker (decimal point)
    _fractional digits_ follow  the fraction_marker (decimal point)
"""
struct ReadableNumStyle
    integral_digits_spanned::Int32
    fractional_digits_spanned::Int32
    between_integral_spans::Char
    between_fractional_spans::Char
    between_parts::Char
end

# default values

const IDIGS = 3%Int32                      # integral_digits_spanned
const FDIGS = 5%Int32                      # fractional_digits_spanned
const IBTWN = FRACPOINT != ',' ? ',' : '.' # between_integral_digits
const FBTWN = '_'                          # between_fractional_digits

# constructors cover likely argument orderings and omisions

ReadableNumStyle() = ReadableNumStyle(IDIGS, FDIGS, IBTWN, FBTWN, FRACPOINT)

function ReadableNumStyle(
             idigs::I, rdigs::I=FDIGS%I, ibtwn::Char=IBTWN, rbtwn::Char=FBTWN, fracpt::Char=FRACPOINT) where {I <: Integer}
    pns = ReadableNumStyle( idigs%Int32, rdigs%Int32, ibtwn, rbtwn, fracpt )
    set_pretty_number_style!(pns)
    return pns
end

ReadableNumStyle(
    ibtwn::Char, rbtwn::Char=FBTWN, idigs::I=IDIGS%I, rdigs::I=FDIGS%I, fracpt::Char=FRACPOINT
    ) where {I <: Integer} =
    ReadableNumStyle(idigs, rdigs, ibtwn, rbtwn, fracpt)
ReadableNumStyle(
    idigs::S, ibtwn::Char, rdigs::S=FDIGS%S, rbtwn::Char=FBTWN, fracpt::Char=FRACPOINT
    ) where {S <: Signed} =
    ReadableNumStyle( idigs%Int32, rdigs%Int32, ibtwn, rbtwn, fracpt )
ReadableNumStyle(
    ibtwn::Char, idigs::S, rbtwn::Char=FBTWN, rdigs::S=FDIGS%S, fracpt::Char=FRACPOINT
    ) where {S <: Signed} =
    ReadableNumStyle(idigs, rdigs, ibtwn, rbtwn, fracpt)

ReadableNumStyle(digs::I, btwn::Char, fracpt::Char=FRACPOINT) where {I <: Integer} =
    ReadableNumStyle( digs%Int32, digs%Int32, btwn, btwn, fracpt )
ReadableNumStyle(btwn::Char, digs::I, fracpt::Char=FRACPOINT) where {I <: Integer} =
    ReadableNumStyle( digs, btwn, fracpt )

# remember the most recent pretty number style

const PRETTY_NUMBER_STYLE_HOLDER = [ ReadableNumStyle() ]
function get_pretty_number_style()
    return PRETTY_NUMBER_STYLE_HOLDER[1]
end
function set_pretty_number_style!(pns::ReadableNumStyle)
    PRETTY_NUMBER_STYLE_HOLDER[1] = pns
    return nothing
end

# accept extended precision numbers and make them become readable

function readable(x::T, pns::ReadableNumStyle=get_pretty_number_style()) where {T <: Real}
    numstr = string(string(x),FRACPOINT)
    return a_readable_number(numstr, pns)
end

function readable(x::String, pns::ReadableNumStyle=get_pretty_number_style())
    numstr = string(x, FRACPOINT)
    return a_readable_number(numstr, pns)
end


text:
using MathOptInterface
const MOI = MathOptInterface
const MOIB = MOI.Bridges
const MOIT = MOI.Test

const BRIDGED = MOIB.full_bridge_optimizer(SCIP.Optimizer(display_verblevel=0), Float64)
const CONFIG_BRIDGED = MOIT.Config(atol=1e-5, rtol=1e-5, exclude=Any[
    MOI.ConstraintDual, MOI.ConstraintName, MOI.DualObjectiveValue, MOI.VariableBasisStatus, MOI.ConstraintBasisStatus,
])

@testset "MOI unit tests" begin
    excluded = copy(MOI_BASE_EXCLUDED)
    append!(excluded, [
        "test_linear_Interval_inactive",
        "test_linear_integration", # Can not delete variable while model contains constraints
        "test_model_ordered_indices", # TODO should fix? ListOf in order of creation
        "test_quadratic_duplicate_terms", # Can not delete variable while model contains constraints
        "test_quadratic_nonhomogeneous", # unsupported by bridge
        "ScalarAffineFunction_ZeroOne",
        "ScalarQuadraticFunction_ZeroOne",
        "VectorAffineFunction_NormOneCone",
        "VectorQuadraticFunction_NormOneCone",
        "test_basic_VectorOfVariables_SecondOrderCone",
        "test_basic_VectorAffineFunction_SecondOrderCone", # SOC tests fail because of lower bound requirement of RHS var.
        "test_basic_VectorQuadraticFunction_SecondOrderCone",
        "test_basic_VectorAffineFunction_RotatedSecondOrderCone",
        "test_basic_VectorQuadraticFunction_RotatedSecondOrderCone",
        "test_basic_VectorOfVariables_RotatedSecondOrderCone",
        "VectorAffineFunction_SOS",
        "VectorQuadraticFunction_SOS",
        "test_conic_GeometricMeanCone_Vector",
        "test_basic_VectorOfVariables_GeometricMeanCone",
        "test_basic_VectorQuadraticFunction_GeometricMeanCone",
        "test_basic_VectorAffineFunction_GeometricMeanCone",
        "test_conic_SecondOrderCone_nonnegative_post_bound",
        "test_variable_delete_SecondOrderCone",
        "test_modification_func_scalaraffine_",
        "test_modification_func_vectoraffine_",
    ])
    MOIT.runtests(
        BRIDGED,
        CONFIG_BRIDGED,
        warn_unsupported=false,
        exclude = excluded,
    )
end


text:
using Unitful
import Unitful: nm, μm, mm, cm, m, km, inch, ft, mi, rad, °
# include("prepare.jl")
mutable struct RI
    water
    cornea
    lens
    lens2distal_retina
    distal_retina2proximal_retina
    gap
    behind_mirror
end

mutable struct Thick
    water
    cornea
    lens
    lens2distal_retina
    distal_retina2proximal_retina
    gap
    behind_mirror
end

mutable struct Radii
    cornea_distal
    lens_distal
    lens_proximal
    distal_retina
    proximal_retina
    mirror
end

ri = RI(zeros(7)...)
thick = Thick(μm*zeros(7)...)
radii = Radii(μm*zeros(6)...)

##################################################
############## EDIT BELOW HERE ###################
##################################################
include("variables.jl")
#= Light properties
source_distance = 1m
source_angle    = 0°

# Eye morphology
aperture        = 251μm
morphing_factor = 1.12

## Volumes
### Refractive indices

ri.water                         = 1.334
ri.cornea                        = 1.37
ri.lens                          = 1.42
ri.lens2distal_retina            = 1.35
ri.distal_retina2proximal_retina = 1.35
ri.gap                           = 1.34

### Thicknesses

thick.cornea                        = 23μm
thick.lens                          = 215μm
thick.lens2distal_retina            = 6μm
thick.distal_retina2proximal_retina = 81μm
thick.gap                           = 49μm

## Interfaces
### Radii of the ellipsoids that form the interface

radii.cornea_distal   = (244μm, 227μm)
radii.lens_distal     = (195μm, 213μm)
radii.lens_proximal   = 337μm
radii.distal_retina   = 337μm
radii.proximal_retina = 337μm
radii.mirror          = 417μm

# Figure ranges
min_aperture = 100μm
max_aperture = 300μm
min_distance = 1cm
max_distance = 1m=#

##################################################
############## EDIT ABOVE HERE ###################
##################################################

# include("plotit.jl")

using Rayden, Scallop

thick.behind_mirror = Inf*μm
ri.behind_mirror = 0.0
thick.water = 0μm

function format_radii{T <: Number}(x::Tuple{T})
    xx1, = x
    x1 = format(xx1)
    return Vec(x1, x1, x1)
end
function format_radii{T <: Number}(x::Tuple{T, T})
    xx1, xx2 = x
    x1 = format(xx1)
    x2 = format(xx2)
    return Vec(x1, x1, x2)
end
function format_radii{T <: Number}(xx::T)
    x = format(xx)
    return Vec(x, x, x)
end

_format(x::T) where T <: Number = uconvert(rad, x)
_format(x::T) where T <: Unitful.Length = uconvert(μm, x)
format(x::T) where T <: Number = Float64(ustrip(_format(x)))
format(x::Int) = Float64(x)
format(x::Float64) = x

morphing_factor = format(morphing_factor)
source_distance = format(source_distance)

volumes = Dict(String(k) => Volume(Float64(getfield(ri, k)), format(getfield(thick, k))) for k in fieldnames(ri))
planes = Dict(String(k) => Plane(format_radii(getfield(radii, k)), !r"distal$"(String(k))) for k in fieldnames(radii))
aperture = format(aperture)

ellipsoids = getmembranes(volumes, planes, morphing_factor, aperture)
opticunits = scallop(ellipsoids, volumes)
l = getLight(source_distance, aperture, ellipsoids["cornea_distal"], format(source_angle))

min_aperture = format(min_aperture)
max_aperture = format(max_aperture)
min_distance = format(min_distance)
max_distance = format(max_distance)

# include("plottingFunctions.jl")

using Plots
plotlyjs()

# function ploteye(ellipsoids::Dict{String,Ellipsoid}, l::Light, opticunits::Vector{OpticUnit}; n_rays::Int = 15)
function ploteye(;n_rays::Int = 15)
    include("main.jl")
    x, y = coordinates2d(ellipsoids)
    plot(x, y, aspect_ratio = :equal, leg=false, color = :black, reuse=false)
    r = Ray()
    for b in linspace(0, 1, n_rays), b2 in [0, 0.5]
        getRay!(r, l, b, b2)
        x = Float64[]
        y = Float64[]
        push!(x, r.orig[1])
        push!(y, r.orig[3])
        complete = true
        for i in opticunits
            if raytrace!(r, i)
                complete = false
                break
            end
            push!(x, r.orig[1])
            push!(y, r.orig[3])
        end
        complete && plot!(x, y, color = :grey)
    end
    plot!(ylims = (-500,100), xlims = (-300,300), xlabel = "μm", ylabel = "μm")
    gui()
end

# function plotpixels(opticunits::Vector{OpticUnit}, l::Light; n_rays::Int = 100000)
function plotpixels(;n_rays::Int = 100000)
    include("main.jl")
    psf = getpsf(opticunits, l, n_rays)
    ph = []
    for (k, v) in psf
        push!(ph, heatmap(collect(v), title = k, aspect_ratio = :equal, reuse = false))
    end
    plot(reuse = false, ph...)
    gui()
end

# function plotfwhm_aperture(opticunits::Vector{OpticUnit}, source_distance::Number, min_aperture::Number, max_aperture::Number; n_data::Int = 10, n_rays::Int = 10000, angle_step = 0.25°)
function plotfwhm_aperture(;n_data::Int = 10, n_rays::Int = 10000, angle_step = 0.25°)
    include("main.jl")
    apertures = linspace(min_aperture, max_aperture, n_data)
    fwhms = Dict(k => zeros(n_data) for k in retinas)
    for (i, aperture) in enumerate(apertures), (k, v) in getfwhm(opticunits, source_distance, aperture, n_rays, format(angle_step))
        fwhms[k][i] = get(v)
    end
    plot(reuse=false)
    [plot!(apertures, v, label = k) for (k, v) in fwhms]
    plot!(xlabel = "Aperture (μm)", ylabel = "FWHM (°)")
    gui()
end

# function plotfwhm_distance(opticunits::Vector{OpticUnit}, aperture::Number, min_distance::Number, max_distance::Number; n_data::Int = 10, n_rays::Int = 10000, angle_step = 0.25°, u = cm)
function plotfwhm_distance(;n_data::Int = 10, n_rays::Int = 10000, angle_step = 0.25°)
    include("main.jl")
    viewdistances = logspace(log10(min_distance), log10(max_distance), n_data)
    fwhms = Dict(k => zeros(n_data) for k in retinas)
    for (i, viewdistance) in enumerate(viewdistances), (k, v) in getfwhm(opticunits, viewdistance, aperture, n_rays, format(angle_step))
        fwhms[k][i] = get(v)
    end
    plot(xscale = :log10, xlabel = "Source distance (cm)", ylabel = "FWHM (°)", reuse=false)
    [plot!(ustrip.(uconvert.(cm, μm.*viewdistances)), v, label = k) for (k, v) in fwhms]
    gui()
end

# function plotfwhm_aperture_distance(opticunits::Vector{OpticUnit}, min_aperture::Number, max_aperture::Number, min_distance::Number, max_distance::Number; n_data::Int = 10, n_rays::Int = 10000, angle_step = 0.25°)
function plotfwhm_aperture_distance(;n_data::Int = 10, n_rays::Int = 10000, angle_step = 0.25°)
    include("main.jl")
    viewdistances = logspace(log10(min_distance), log10(max_distance), n_data)
    apertures = linspace(min_aperture, max_aperture, n_data)
    fwhms = Dict(k => zeros(n_data,n_data) for k in retinas)
    for (i, viewdistance) in enumerate(viewdistances), (j, aperture) in enumerate(apertures), (k, v) in getfwhm(opticunits, viewdistance, aperture, n_rays, format(angle_step))
        fwhms[k][i,j] = get(v)
    end
    ph = []
    for (k, v) in fwhms
        push!(ph, heatmap(apertures, ustrip.(uconvert.(cm, μm.*viewdistances)), v, title = k, xlabel = "Aperture (μm)", ylabel = "Source distance (cm)", reuse = false))
    end
    plot(ph..., yscale = :log10, zlims = (0,20), reuse=false)
    gui()
end

#=l = getLight(1e9, aperture, ellipsoids["cornea_distal"], .05)
psf = get2dpsf(opticunits, l, 100000)
plot()
for (k, v) in psf
    plot!(collect(v), label = k)
end
gui()=#

#=nt = 20
thetas = linspace(0, .5, nt)
ratios = Dict(k => zeros(nt) for k in retinas)
for (i, θ) in enumerate(thetas)
    l = getLight(1e9, aperture, ellipsoids["cornea_distal"], θ)
    psf = get2dpsf(opticunits, l, 100000)
    for k in keys(ratios)
        ratios[k][i] = maximum(psf[k])/psf[k][0]
    end
end
plot(ylims = (0,5))
for (k, v) in ratios
    plot!(thetas, v, label = k)
end
gui()=#


text:
using Flux
_convshift(n) = (i = div(n, 2); mod(n, 2) == 0 ? (1 - i:i) : (-i : i) )
"""
	_addmatvec!(o, i, W, x, j)

	add a product of matrix `W` with a j-th column of `x` to i-th columns of `o`

"""
function _addmatvec!(o::Matrix, i, W::Matrix, x::Matrix, j)
	@inbounds for s in 1:size(W, 2)
		for r in 1:size(W, 1)
			o[r, i] += W[r, s] * x[s, j]
		end
	end
end

function _addmatvec!(o::Matrix, i, W::Matrix, x::SparseMatrixCSC, j)
	rn = x.colptr[j]:x.colptr[j+1] - 1
	isempty(rn) && return
	rowptr = x.rowval[rn]
	vals = x.nzval[rn]
	@inbounds for (s, v) in zip(rowptr, vals)
		for r in 1:size(W, 1)
			o[r, i] += W[r, s] * v
		end
	end
end

"""
	_addmattvec!(o, i, W, x, j)

	add a product of a transposed matrix `W` with a j-th column of `x` to i-th columns of `o`

"""
function _addmattvec!(o::Matrix, i, W::Matrix, x::AbstractMatrix, j)
	@inbounds for s in 1:size(W, 1)
		for r in 1:size(W, 2)
			o[r, i] += W[s, r] * x[s, j]
		end
	end
end

"""
	_outeradd!(W, Δ, i, x, j)

	add an outer product of i-th column of `Δ` and transposed `j`-th columns of `x` to `W`

"""
function _addvecvect!(W::Matrix, Δ::Matrix, i, x::AbstractMatrix, j)
	@inbounds for r in 1:size(W, 2)
		for s in 1:size(W, 1)
			W[s, r] += Δ[s, i] * x[r, j]
		end
	end
end

function _addvecvect!(W::Matrix, Δ::Matrix, i, x::SparseMatrixCSC, j)
	rn = x.colptr[j]:x.colptr[j+1] - 1
	isempty(rn) && return
	rowptr = x.rowval[rn]
	vals = x.nzval[rn]
	@inbounds for (r, v) in zip(rowptr, vals)
		for s in 1:size(W, 1)
			W[s, r] += Δ[s, i] * v
		end
	end
end

function bagconv!(o, x, bags::T, W...) where {T<:Union{AlignedBags, Vector{UnitRange{Int64}}}}
	offsets = _convshift(length(W))
	for b in bags
		for ri in b
			o[:, ri] .= 0
			for (i, k) in enumerate(offsets)
				if first(b) <= k + ri  <= last(b)
					_addmatvec!(o, ri, W[i], x, k + ri)
				end
			end
		end
	end
	o
end

function bagconv!(o, x, bags::T, W...) where {T<:Union{ScatteredBags,Vector{Vector{Int64}}}}
	offsets = _convshift(length(W))
	for b in bags
		for (bi, ri) in enumerate(b)
			o[:, ri] .= 0
			for (i, k) in enumerate(offsets)
				if 0 < k + bi  <= length(b)
					_addmatvec!(o, ri, W[i], x, b[k+bi])
				end
			end
		end
	end
	o
end

function bagconv(x, bags, W...)
	o = similar(W[1], size(W[1], 1), size(x, 2))
	o .= 0
	for i in 1:Threads.nthreads()
		bagconv!(o, x, subsetof(bags), W...)
	end
	o
end

function ∇wbagconv!(∇W, Δ, x, bags::T, W...) where {T<:Union{AlignedBags, Vector{UnitRange{Int64}}}}
	offsets = _convshift(length(W))
	foreach(w -> w .= 0, ∇W)
	for b in bags
		for ri in b
			for (i, k) in enumerate(offsets)
				if first(b) <= k + ri  <= last(b)
					_addvecvect!(∇W[i], Δ, ri, x, k + ri)
				end
			end
		end
	end
end

function ∇wbagconv!(∇W, Δ, x, bags::T, W...) where {T<:Union{ScatteredBags,Vector{Vector{Int64}}}}
	offsets = _convshift(length(W))
	foreach(w -> w .= 0, ∇W)
	for b in bags
		for (bi, ri) in enumerate(b)
			for (i, k) in enumerate(offsets)
				if 0 < k + bi  <= length(b)
					_addvecvect!(∇W[i], Δ, ri, x, b[k+bi])
				end
			end
		end
	end
end

function ∇wbagconv(Δ, x, bags, W...)
	∇Ws = [[zero(w) for w in W] for i in 1:Threads.nthreads()]
	for i in 1:Threads.nthreads()
		∇wbagconv!(∇Ws[i], Δ, x, subsetof(bags), W...)
	end
	foreach(i -> foreach(x -> x[1] .+= x[2], zip(∇Ws[1], ∇Ws[i])), 2:Threads.nthreads())
	return(tuple(∇Ws[1]...))
end

function ∇xbagconv(Δ, x, bags::T, W...) where {T<:Union{AlignedBags, Vector{UnitRange{Int64}}}}
	offsets = _convshift(length(W))
	∇x = fill(zero(eltype(x)), size(x))
	for b in bags
		for ri in b
			for (i, k) in enumerate(offsets)
				if first(b) <= k + ri  <= last(b)
					_addmattvec!(∇x, k + ri, W[i], Δ, ri)
				end
			end
		end
	end
	∇x
end

function ∇xbagconv(Δ, x, bags::T, W...) where {T<:Union{ScatteredBags,Vector{Vector{Int64}}}}
	offsets = _convshift(length(W))
	∇x = fill(zero(eltype(x)), size(x))
	for b in bags
		for (bi, ri) in enumerate(b)
			for (i, k) in enumerate(offsets)
				if 0 < k + bi  <= length(b)
					_addmattvec!(∇x, b[k+bi], W[i], Δ, ri)
				end
			end
		end
	end
	∇x
end

function ∇xwbagconv(Δ, x, bags::T, W...) where {T<:Union{AlignedBags, Vector{UnitRange{Int64}}}}
	offsets = _convshift(length(W))
	∇x = fill(zero(eltype(x)), size(x))
	∇W = [zero(w) for w in W]
	for b in bags
		for ri in b
			for (i, k) in enumerate(offsets)
				if first(b) <= k + ri  <= last(b)
					_addmattvec!(∇x, k + ri, W[i], Δ, ri)
					_addvecvect!(∇W[i], Δ, ri, x, k + ri)
				end
			end
		end
	end
	∇x, tuple(∇W...)
end

function ∇xwbagconv(Δ, x, bags::T, W...) where {T<:Union{ScatteredBags,Vector{Vector{Int64}}}}
	offsets = _convshift(length(W))
	∇x = zero(W[1])
	∇W = [zero(w) for w in W]
	for b in bags
		for (bi, ri) in enumerate(b)
			for (i, k) in enumerate(offsets)
				if first(b) <= k + ri  <= last(b)
					_addmattvec!(∇x, b[k+bi], W[i], Δ, ri)
					_addvecvect!(∇W[i], Δ, ri, x, b[k+bi])
				end
			end
		end
	end
	∇x, tuple(∇W...)
end

Zygote.@adjoint function bagconv(x, bags, fs::Matrix...)
  bagconv(x, bags, fs...), Δ -> (∇xbagconv(Δ, x, bags,  fs...), nothing,  ∇wbagconv(Δ, x, bags,  fs...)...)
end

"""
	struct BagConv{T, F}
		W::T
		σ::F
	end

	Convolution over a matrix `X` correctly handing borders between bags. The convolution is little bit special, as it
	assumes that input is always a matrix (never a tensor) and the kernel always spans the full dimension of the vector.

	BagConv(d::Int, o::Int, n::Int, σ = identity)
	`d` --- input dimension
	`o` --- output dimension (number of channels)
	`n` --- size of convolution
	`σ` --- transfer function used after the convolution

	note that of `n` equals one, then the convolution boils down to multiplication of input data `x` with a matrix `W`
"""
struct BagConv{T, F}
	W::T
	σ::F
end

Flux.@functor BagConv

function BagConv(d::Int, o::Int, n::Int, σ = identity)
	W = (n > 1) ? tuple([randn(o, d) .* sqrt(2.0/(o + d)) for _ in 1:n]...) : randn(o, d) .* sqrt(2.0/(o + d))
	BagConv(W, σ)
end

(m::BagConv{T,F} where {T<:Tuple,F})(x, bags) = m.σ.(bagconv(x, bags, m.W...))
(m::BagConv{T,F} where {T<:AbstractMatrix,F})(x, bags) = m.σ.(m.W * x)
(m::BagConv)(x::ArrayNode, bags::AbstractBags) = ArrayNode(bagconv(x.data, bags, m.W...))
(m::BagConv)(x::BagNode) = ArrayNode(bagconv(x.data.data, x.bags, m.W...))

convsum(bags, xs::AbstractMatrix) = xs
function convsum(bags, xs...)
	offsets = _convshift(length(xs))
	o = similar(xs[1]) .= 0
	for b in bags
		for ri in b
			for (i, k) in enumerate(offsets)
				if first(b) <= k + ri  <= last(b)
					o[:, ri] .+= view(xs[i], :, k + ri)
				end
			end
		end
	end
	o
end

function ∇convsum(Δ, bags, n)
	offsets = _convshift(n)
	o = [zero(Δ) for i in 1:n]
	for b in bags
		for ri in b
			for (i, k) in enumerate(offsets)
				if first(b) <= k + ri  <= last(b)
					o[i][:, k + ri] .+= view(Δ, :, ri)
				end
			end
		end
	end
	tuple(o...)
end

Zygote.@adjoint function convsum(bags, xs...)
  convsum(bags, xs...), Δ -> (nothing,  ∇convsum(Δ, bags, length(xs))...)
end

legacy_bagconv(x, bags, f::AbstractArray{T, 3}) where {T} = convsum(bags, [f[:, :, i]*x for i in 1:size(f, 3)]...)


text:
for N in [Float64, Rational{Int}, Float32]
    # random ball
    rand(BallInf)

    # 1D BallInf
    b = BallInf(N[0], N(1))
    # Test Dimension
    @test dim(b) == 1
    # Test Support Vector
    d = N[1]
    @test σ(d, b) == N[1]
    d = N[-1]
    @test σ(d, b) == N[-1]

    # 2D BallInf
    b = BallInf(N[0, 0], N(1))
    # Test Dimension
    @test dim(b) == 2
    # Test Support Vector
    d = N[1, 1]
    @test σ(d, b) == N[1, 1]
    d = N[-1, 1]
    @test σ(d, b) == N[-1, 1]
    d = N[-1, -1]
    @test σ(d, b) == N[-1, -1]
    d = N[1, -1]
    @test σ(d, b) == N[1, -1]

    # 2D BallInf not 0-centered
    b = BallInf(N[1, 2], N(1))
    # Test Dimension
    @test dim(b) == 2
    # Test Support Vector
    d = N[1, 1]
    @test σ(d, b) == N[2, 3]
    d = N[-1, 1]
    @test σ(d, b) == N[0, 3]
    d = N[-1, -1]
    @test σ(d, b) == N[0, 1]
    d = N[0, -1]
    @test σ(d, b) == N[2, 1]

    # 2D BallInf radius =/= 1
    b = BallInf(N[0, 0], N(2))
    # Test Dimension
    @test dim(b) == 2
    # Test Support Vector
    d = N[1, 1]
    @test σ(d, b) == N[2, 2]
    d = N[-1, 1]
    @test σ(d, b) == N[-2, 2]
    d = N[-1, -1]
    @test σ(d, b) == N[-2, -2]
    d = N[1, -1]
    @test σ(d, b) == N[2, -2]

    # center
    c = N[1, 2]
    b = BallInf(c, N(2))
    @test center(b) == c && center(b, 1) == N(1) && center(b, 2) == N(2)

    # support vector for single entry vector
    svec = σ(SingleEntryVector(2, 3, N(2)), BallInf(zeros(N, 3), N(2)))
    @test svec[1] ∈ Interval(N[-2, 2]) && svec[2] == N(2) && svec[3] ∈ Interval(N[-2, 2])

    # support function
    B = BallInf(N[1, 2], N(1))
    @test ρ(N[1, 1], B) == N(5)
    @test ρ(N[1, 0], B) == N(2)
    @test ρ(N[0, 1], B) == N(3)
    @test ρ(N[-1, -1], B) == N(-1)
    @test ρ(N[-1, 0], B) == N(0)
    @test ρ(N[-1, 1], B) == N(3)

    # boundedness
    @test isbounded(b)

    # isempty
    @test !isempty(b)

    # isuniversal
    answer, w = isuniversal(b, true)
    @test !isuniversal(b) && !answer && w ∉ b

    # membership
    b = BallInf(N[1, 1], N(1))
    @test N[0.5, -0.5] ∉ b
    @test N[0.5, 1.5] ∈ b

    # an_element function
    b = BallInf(N[1, 2], N(3))
    @test an_element(b) ∈ b

    # check that vertices_list for zero radius doesn't repeat vertices
    b = BallInf(N[1, 2], N(0))
    vl = vertices_list(b)
    @test vl == [center(b)]

    # high and low
    b = BallInf(N[1, 2], N(1))
    @test high(b) == N[2, 3]
    @test low(b) == N[0, 1]

    # isflat
    @test !isflat(BallInf(N[1, 1], N(1))) && isflat(BallInf(N[1, 1], N(0)))

    # split
    b = BallInf(N[3, 3], N(1))
    @test split(b, [1, 1]) == [Hyperrectangle(N[3, 3], N[1, 1])]
    @test ispermutation(split(b, [2, 1]),
        [Hyperrectangle(N[2.5, 3], N[0.5, 1]),
         Hyperrectangle(N[3.5, 3], N[0.5, 1])])

    # split along the second generator if we see the BallInf as a zonotope
    B2 = BallInf(N[1, 1], N(1))
    Z2 = split(B2, 2)
    @test Z2[1] == Zonotope(N[1, 0.5], N[1 0; 0 0.5])
    @test Z2[2] == Zonotope(N[1, 1.5], N[1 0; 0 0.5])

    # repaeated split method if we see the BallInf as a zonotope
    Z4 = split(B2, [1, 2], [1, 1])
    @test Z4[1] == Zonotope(N[0.5, 0.5], N[0.5 0; 0 0.5])
    @test Z4[2] == Zonotope(N[0.5, 1.5], N[0.5 0; 0 0.5])
    @test Z4[3] == Zonotope(N[1.5, 0.5], N[0.5 0; 0 0.5])
    @test Z4[4] == Zonotope(N[1.5, 1.5], N[0.5 0; 0 0.5])

    # translation
    b = BallInf(N[1, 2], N(1))
    @test translate(b, N[1, 2]) == BallInf(N[2, 4], N(1))

    # affine map
    M = N[2 1; 0 2]
    B = BallInf(N[0, 0], N(1))
    v = N[1, -1]
    am = affine_map(M, B, v)
    @test ispermutation(vertices_list(am),
                        [N[4, 1], N[0, 1], N[-2, -3], N[2, -3]])
    amv = affine_map(M, B, v, algorithm="vrep") # pass a custom algorithm
    @test amv isa VPolygon && isequivalent(am, amv)

    # volume
    B = BallInf(N[0, 0], N(1))
    @test volume(B) ≈ N(4)
    if N <: AbstractFloat
        B = BallInf(zeros(N, 100), N(1/2 + 1e-5))
        @test volume(B) ≈ N(1.0020019812942185)
    end

    # area
    B = BallInf(N[0, 0], N(1))
    @test area(B) ≈ N(4)

    # concretize
    B = BallInf(N[0, 0], N(1))
    @test concretize(B) === B

    # constraints iterator
    @test ispermutation(collect(constraints(B)), constraints_list(B))

    # vertices iterator
    @test ispermutation(collect(vertices(B)), vertices_list(B))

    # conversion to a zonoope
    B = BallInf(N[0, 0], N(1))
    ZB = convert(Zonotope, B)
    @test ZB == Zonotope(N[0, 0], N[1 0; 0 1])
    B = BallInf(N[0, 0], N(0))  # flat case
    ZB = convert(Zonotope, B)
    @test ZB == Zonotope(N[0, 0], Matrix{N}(undef, 2, 0))

   # conversion to a zonotope, static arrays
   B = BallInf(SA[N(0), N(0)], N(1))
   ZB = convert(Zonotope, B)
   @test ZB == Zonotope(SVector{2}(N[0, 0]), SMatrix{2, 2}(N[1 0; 0 1]))
   B = BallInf(SA[N(0), N(0)], N(0))  # flat case
   ZB = convert(Zonotope, B)
   @test ZB == Zonotope(SVector{2}(N[0, 0]), SMatrix{2, 0, N, 0}())
   # specialized method (no prunning)
   B = BallInf(SA[1.0, 2.0], 1.0)
   ZB = LazySets._convert_2D_static(Zonotope, B)
   @test ZB == Zonotope(SA[1.0, 2.0], SA[1.0 0.0; 0.0 1.0])

   # internal function
   B = BallInf(SA[N(0), N(0)], N(1))
   Zs = LazySets._convert_2D_static(Zonotope, B)
   @test Zs == Zonotope(SVector{2}(N[0, 0]), SMatrix{2, 2}(N[1 0; 0 1]))
end

# tests that only work with Float64
for N in [Float64]
    # concrete Minkowski sum
    b = BallInf(N[1, 2], N(1))
    p = minkowski_sum(b, N[2 0; 0 1] * b)
    @test p isa HPolytope{N} && ispermutation(constraints_list(p),
        [HalfSpace(N[0, -1], N(-2)), HalfSpace(N[0, 1], N(6)),
         HalfSpace(N[-1, 0], N(0)), HalfSpace(N[1, 0], N(6))])
end


text:
using NavierStokes2D
using Test

#=
grid = Grid(1., 2., 8, 16)
params = Parameters(0.1, 0.2, (0.0, 1.0))
include("analytic_functions.jl")
u₀ = make_gaussian(grid)
u_analytic(t) = solution_gaussian(grid, t, params.ν)
model = DiffusionModel(zeros(8, 16), grid, params)
model_plus_ana = DiffusionModel(u₀, u_analytic, grid, params)
prob2 = DiffusionProblem(model_plus_ana, CrankNickolsonDiffusion())
sol = solve(prob2)
=#

@testset "Grid tests" begin

    grid = Grid(1., 2., 16, 32)
    @test grid.xmax == 1.0
    @test grid.ymax == 2.0
    @test grid.nx == 16
    @test grid.ny == 32
    @test size(grid.xRange) == (16,)
    @test size(grid.yRange) == (32,)
    @test isa(grid.Δx, Float64)
    @test isa(grid.Δy, Float64)

end

@testset "Parameters tests" begin

    params = Parameters(0.1, 0.2, (0.0, 1.0))
    @test params.ν == 0.1
    @test params.σ == 0.2
    @test length(params.tspan) == 2

end

@testset "DiffusionModel tests" begin

    grid = Grid(1., 2., 16, 32)
    params = Parameters(0.1, 0.2, (0.0, 1.0))
    model = DiffusionModel(zeros(32, 16), grid, params)
    @test model.u_analytic == nothing
    @test isa(model.Δt, Float64)
    @test isa(model.t, Vector{Float64})

    include("analytic_functions.jl")
    u₀ = make_gaussian(grid)
    u_analytic(t) = solution_gaussian(grid, t, params.ν)
    model_plus_ana = DiffusionModel(u₀, u_analytic, grid, params)
    @test isa(model_plus_ana.u_analytic, Function)
    @test isa(model_plus_ana.Δt, Float64)
    @test isa(model_plus_ana.t, Vector{Float64})

end

@testset "Cache tests" begin

    u = zeros(16, 32)
    explicitDiffusionCache = ExplicitDiffusionCache(u, copy(u), 0.0)
    @test explicitDiffusionCache.u == u
    @test explicitDiffusionCache.uprev == u
    @test explicitDiffusionCache.t == 0.0

    grid = Grid(1., 2., 16, 32)
    params = Parameters(0.1, 0.2, (0.0, 1.0))
    model = DiffusionModel(zeros(32, 16), grid, params)
    cache = init_cache(model, ExplicitDiffusion())
    @test isa(cache, ExplicitDiffusionCache)

    # NOTE: Assume other caches test similarly

end

@testset "Diffusion methods tests" begin

    grid = Grid(1., 2., 16, 32)
    params = Parameters(0.1, 0.2, (0.0, 1.0))
    model = DiffusionModel(zeros(32, 16), grid, params)

    set_method = ExplicitDiffusion()
    cache = init_cache(model, ExplicitDiffusion())

    perform_step!(model, set_method, cache)

    @test isa(model, DiffusionModel)

    # TODO: Add more explicit tests here
    # NOTE: Assume other methods test similarly

end

@testset "Diffusion problems tests" begin

    grid = Grid(1., 2., 16, 32)
    params = Parameters(0.1, 0.2, (0.0, 1.0))
    model = DiffusionModel(zeros(32, 16), grid, params)

    set_method = ExplicitDiffusion()
    cache = init_cache(model, ExplicitDiffusion())

    prob = DiffusionProblem(model, set_method, cache)

    @test isa(prob, DiffusionProblem)
    # NOTE: Assume other methods test similarly

    prob1 = DiffusionProblem(model, ExplicitDiffusion())
    prob2 = DiffusionProblem(model, CrankNickolsonDiffusion())
    prob3 = DiffusionProblem(model, ADIDiffusion())

    @test typeof(prob1).parameters[1] == ExplicitDiffusion
    @test typeof(prob2).parameters[1] == CrankNickolsonDiffusion
    @test typeof(prob3).parameters[1] == ADIDiffusion
    @test typeof(prob1).parameters[2] == ExplicitDiffusionCache
    @test typeof(prob2).parameters[2] == CNDiffusionCache
    @test typeof(prob3).parameters[2] == ADIDiffusionCache

end

@testset "Diffusion solution tests" begin

    grid = Grid(1., 2., 8, 16)
    params = Parameters(0.1, 0.2, (0.0, 1.0))
    model = DiffusionModel(zeros(16, 8), grid, params)
    prob = DiffusionProblem(model, ExplicitDiffusion())

    sol = DiffusionSolution(prob, true)

    @test sol.u_analytic == nothing
    @test sol.errors == nothing

    @test length(sol.prob.model.t) == 7
    sol2 = solve(prob)
    @test isa(sol2, DiffusionSolution{Nothing})
    @test sol2.u_analytic == nothing
    @test sol2.errors == nothing
    # NOTE: Assume other methods test similarly

    include("analytic_functions.jl")
    u₀ = make_gaussian(grid)
    u_analytic(t) = solution_gaussian(grid, t, params.ν)
    model_plus_ana = DiffusionModel(u₀, u_analytic, grid, params)
    prob_plus_ana = DiffusionProblem(model_plus_ana, ExplicitDiffusion())
    sol_plus_ana = DiffusionSolution(prob_plus_ana, true)

    @test isa(sol_plus_ana.u_analytic, Array{Float64, 2})
    @test isa(sol_plus_ana.errors, Array{Float64, 2})

    sol2_plus_ana = solve(prob_plus_ana)
    @test isa(sol2_plus_ana, DiffusionSolution{Array{Float64, 2}})
    @test isa(sol2_plus_ana.u_analytic, Array{Float64, 2})
    @test isa(sol2_plus_ana.errors, Array{Float64, 2})
    # NOTE: Assume other methods test similarly

end


text:
function full_pressure()

    # Create material with flat white spectrum
    whitespec = SampledSpectrum(300, 800, [1.0, 1.0, 1.0])
    mat = Lambert(whitespec)

    # Create sphere shape
    sph = Sphere(1, 1.0, translation(0.0, 1.0, 0.0))
    disk = Disk(translation(2, 0, 0))

    # Create primitive combining shape and material
    sph_primitive = GeometricPrimitive(sph, mat, nothing, 1)
    dsk_primitive = GeometricPrimitive(disk, mat, nothing, 1)

    # Create list of primitives
    primitives = GeometricPrimitive[sph_primitive, dsk_primitive]

    # Generate bounding box hierarchy from primitives
    stuff = BVHAccelerator(primitives)

    # Create a distant light source with a single-line 532 nm spectrum
    p_light = DistantLight(SingleLine(532.0, 1.0), rotation(X_AXIS, 0.0))
    p_light2 = DistantLight(SingleLine(532.0, 1.0), rotation(X_AXIS, pi))
    lights = LightSource[p_light]

    # Create scene with bounding box hierarchy and light source
    scene = Scene(stuff, lights)

    # Set up pixel sampler with 9 rays per pixel (3 rays in both X and Y directions)
    res = 256
    sampler = StratifiedSampler(res, res, 5)

    # Declare an integrator with maximum ray depth of 1, using 100000 rays per light source
    # and initialize empty arrays to store the force and torque from each light source.
    forces = [Vector3(0) for i in lights]
    torques = [Vector3(0) for i in lights]
    counts = [0 for i in lights]
    integrator = PressureIntegrator(1, forces, torques, counts)

    # Run the renderer
    @time render(scene, integrator, sampler ; debug=true)

    forces = forces ./ counts
    torques = torques ./ counts

    return forces, torques
end

@testset "Full pressure computation" begin
    f, t = full_pressure()
    @test isapprox(f[1].z, -13/9*π  - 5/3 * π ; atol=5)
end

function setup_disk_pressure(angle::Real)
    whitespec = SampledSpectrum(300, 800, [1.0, 1.0, 1.0])
    mat = Lambert(whitespec)

    # Create sphere shape
    disk = Disk(translation(2, 0, 0))

    # Create primitive combining shape and material
    dsk_primitive = GeometricPrimitive(disk, mat, nothing, 1)

    # Create list of primitives
    primitives = GeometricPrimitive[dsk_primitive]

    # Generate bounding box hierarchy from primitives
    stuff = BVHAccelerator(primitives)

    # Create a distant light source with a single-line 532 nm spectrum
    # and intensity of 2.0 Watts / m^2 / sr
    rot = rotation(X_AXIS, angle)
    p_light = DiskLight(rot(Z_AXIS), 0.01, SingleLine(532.0, 2.0))
    lights = LightSource[p_light]

    # Create scene with bounding box hierarchy and light source
    scene = Scene(stuff, lights)

    # Set up pixel sampler with 9 rays per pixel (5 rays in both X and Y directions)
    res = 256
    sampler = StratifiedSampler(res, res, 5)

    # Declare an integrator with maximum ray depth of 1, using 100000 rays per light source
    # and initialize empty arrays to store the force and torque from each light source.
    forces = [Vector3(0) for i in lights]
    torques = [Vector3(0) for i in lights]
    counts = [0 for i in lights]
    integrator = PressureIntegrator(1, forces, torques, counts)

    return scene, integrator, sampler
end

@testset "Disk pressure with varying incidence" begin

    # Solid angle of the light source times intensity
    I0 = 2.0
    F0 = 2π * (1 - cos(0.01)) * I0

    # Test against analytical solution for various angles
    angles = [0, π/4, π/3, π/2]
    results = [-F0 * cos(angle) * (cos(angle) + 2//3) * π for angle in angles]

    for (angle, expected) in zip(angles, results)
        scene, integrator, sampler = setup_disk_pressure(angle)
        render(scene, integrator, sampler ; debug=true)
        force = integrator.force[1] / integrator.counts[1]

        @test isapprox(force.z, expected ; atol=1e-5)
    end

end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
# ---
# title: 1390. Four Divisors
# id: problem1390
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/four-divisors/description/>
# hidden: true
# ---
#
# Given an integer array `nums`, return the sum of divisors of the integers in
# that array that have exactly four divisors.
#
# If there is no such integer in the array, return `0`.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [21,4,7]
#     Output: 32
#     Explanation:
#     21 has 4 divisors: 1, 3, 7, 21
#     4 has 3 divisors: 1, 2, 4
#     7 has 2 divisors: 1, 7
#     The answer is the sum of divisors of 21 only.
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 10^4`
#   * `1 <= nums[i] <= 10^5`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module Collections

using Artifacts: @artifact_str, load_artifacts_toml
using DocStringExtensions: METHODLIST, TYPEDEF, TYPEDFIELDS
using LazyArtifacts

# export public types and constructors
export GriddedCollection
export canopy_height_collection, clumping_index_collection, elevation_collection, gpp_collection, lai_collection, land_mask_collection, leaf_chlorophyll_collection, leaf_nitrogen_collection,
       leaf_phosphorus_collection, pft_collection, sif_collection, soil_color_collection, soil_hydraulics_collection, sla_collection, surface_area_collection, tree_density_collection,
       vcmax_collection, wood_density_collection

# export public functions
export clean_collections!, query_collection

# collection types
"""
$(TYPEDEF)

Structure for general gridded dataset collection.

# Fields
$(TYPEDFIELDS)

---
# Examples
```julia
vcmax_collection = GriddedCollection("VCMAX", ["2X_1Y_V1", "2X_1Y_V2"], "2X_1Y_V2");
```
"""
struct GriddedCollection
    "Artifact label name"
    LABEL::String
    "Supported combinations"
    SUPPORTED_COMBOS::Vector{String}
    "Default combination"
    DEFAULT_COMBO::String
end

# constructors for GriddedCollection
"""
    canopy_height_collection()

<details>
<summary>
Method to create a general dataset collection for canopy height. Supported datasets are (click to view bibtex items)
- `20X_1Y_V1` [(Simard et al., 2011)](https://doi.org/10.1029/2011JG001708)
- `2X_1Y_V2` [(Boonman et al., 2020)](https://doi.org/10.1111/geb.13086)
</summary>

```
@article{simard2011mapping,
    author = {Simard, Marc and Pinto, Naiara and Fisher, Joshua B and Baccini, Alessandro},
    year = {2011},
    title = {Mapping forest canopy height globally with spaceborne lidar},
    journal = {Journal of Geophysical Research: Biogeosciences},
    volume = {116},
    number = {G4021}
}
@article{boonman2020assessing,
    author = {Boonman, Coline CF and Ben{\\'i}tez-L{\\'o}pez, Ana and Schipper, Aafke M and Thuiller, Wilfried and Anand, Madhur and Cerabolini, Bruno EL and Cornelissen, Johannes HC and
              Gonzalez-Melo, Andres and Hattingh, Wesley N and Higuchi, Pedro and others},
    year = {2020},
    title = {Assessing the reliability of predicted plant trait distributions at the global scale},
    journal = {Global Ecology and Biogeography},
    volume = {29},
    number = {6},
    pages = {1034--1051}
}
```
</details>
"""
canopy_height_collection() = GriddedCollection("CH", ["20X_1Y_V1", "2X_1Y_V2"], "20X_1Y_V1");

"""
    clumping_index_collection()

<details>
<summary>
Method to create a general dataset collection for clumping index. Supported datasets are (click to view bibtex items)
- `240X_1Y_V1` [(He et al., 2012)](https://doi.org/10.1016/j.rse.2011.12.008)
- `2X_1Y_V1` [(regridded; He et al., 2012)](https://doi.org/10.1016/j.rse.2011.12.008)
- `2X_1Y_V2` [(Braghiere et al., 2019)](https://doi.org/10.1029/2018GB006135)

V2 dataset are classified for different plant functional types. The indices are Broadleaf, Needleleaf, C3 grasses, C4 grasses, and shrubland.
</summary>

```
@article{he2012global,
    author={He, Liming and Chen, Jing M and Pisek, Jan and Schaaf, Crystal B and Strahler, Alan H},
    year={2012},
    title={Global clumping index map derived from the MODIS BRDF product},
    journal={Remote Sensing of Environment},
    volume={119},
    pages={118--130}
}
@article{braghiere2019underestimation,
    author = {Braghiere, R{\\'e}nato Kerches and Quaife, T and Black, E and He, L and Chen, JM},
    year = {2019},
    title = {Underestimation of global photosynthesis in Earth System Models due to representation of vegetation structure},
    journal = {Global Biogeochemical Cycles},
    volume = {33},
    number = {11},
    pages = {1358--1369}
}
```
</details>
"""
clumping_index_collection() = GriddedCollection("CI", ["240X_1Y_V1", "2X_1Y_V1", "2X_1Y_V2"], "2X_1Y_V1");

"""
    elevation_collection()

<details>
<summary>
Method to create a general dataset collection for surface elevation. Supported datasets are (click to view bibtex items)
- `2X_1Y_V1` [(Yamazaki et al., 2017)](https://doi.org/10.1002/2017GL072874)
</summary>

```
@article{yamazaki2017high,
    author = {Yamazaki, Dai and Ikeshima, Daiki and Tawatari, Ryunosuke and Yamaguchi, Tomohiro and O'Loughlin, Fiachra and Neal, Jeffery C and Sampson, Christopher C and Kanae, Shinjiro and
              Bates, Paul D},
    year = {2017},
    title = {A high-accuracy map of global terrain elevations},
    journal = {Geophysical Research Letters},
    volume = {44},
    number = {11},
    pages = {5844--5853}
}
```
</details>
"""
elevation_collection() = GriddedCollection("ELEV", ["4X_1Y_V1"], "4X_1Y_V1");

"""
    gpp_collection()

<details>
<summary>
Method to create a general dataset collection for gross primary productivity. Supported datasets are (click to view bibtex items)
- `MPI_RS_2X_1M_YYYY_V1` [(YYYY from 2001 to 2019; Tramontana et al., 2016)](https://doi.org/10.5194/bg-13-4291-2016)
- `MPI_RS_2X_8D_YYYY_V1` [(YYYY from 2001 to 2019; Tramontana et al., 2016)](https://doi.org/10.5194/bg-13-4291-2016)
- `VPM_5X_8D_YYYY_V2` [(YYYY from 2000 to 2019; Zhang et al., 2017)](https://doi.org/10.1038/sdata.2017.165)
- `VPM_12X_8D_YYYY_V2` [(YYYY from 2000 to 2019; Zhang et al., 2017)](https://doi.org/10.1038/sdata.2017.165)
</summary>

```
@article{tramontana2016predicting,
    author = {Tramontana, Gianluca and Jung, Martin and Schwalm, Christopher R and Ichii, Kazuhito and Camps-Valls, Gustau and R{\'a}duly, Botond and Reichstein, Markus and Arain, M Altaf and
              Cescatti, Alessandro and Kiely, Gerard and others},
    year = {2016},
    title = {Predicting carbon dioxide and energy fluxes across global FLUXNET sites with regression algorithms},
    journal = {Biogeosciences},
    volume = {13},
    number = {14},
    pages = {4291--4313}
}
@article{zhang2017global,
    author = {Zhang, Yao and Xiao, Xiangming and Wu, Xiaocui and Zhou, Sha and Zhang, Geli and Qin, Yuanwei and Dong, Jinwei},
    year = {2017},
    title = {A global moderate resolution dataset of gross primary production of vegetation for 2000--2016},
    journal = {Scientific data},
    volume = {4},
    pages = {170165}
}
```
</details>
"""
gpp_collection() = (
    _supported = [];
    for _year in 2001:2019
        push!(_supported, "MPI_RS_2X_1M_$(_year)_V1");
        push!(_supported, "MPI_RS_2X_8D_$(_year)_V1");
    end;
    for _year in 2000:2019
        push!(_supported, "VPM_5X_8D_$(_year)_V2");
        push!(_supported, "VPM_12X_8D_$(_year)_V2");
    end;

    return GriddedCollection("GPP", _supported, "MPI_RS_2X_1M_2019_V1")
);

"""
    lai_collection()

<details>
<summary>
Method to create a general dataset collection for leaf area index. Supported datasets are (click to view bibtex items)
- `MODIS_2X_1M_YYYY_V1` [(YYYY from 2000 to 2020; Yuan et al., 2011)](https://doi.org/10.1016/j.rse.2011.01.001)
- `MODIS_2X_8D_YYYY_V1` [(YYYY from 2000 to 2020; Yuan et al., 2011)](https://doi.org/10.1016/j.rse.2011.01.001)
- `MODIS_10X_1M_YYYY_V1` [(YYYY from 2000 to 2020; Yuan et al., 2011)](https://doi.org/10.1016/j.rse.2011.01.001)
- `MODIS_10X_8D_YYYY_V1` [(YYYY from 2000 to 2020; Yuan et al., 2011)](https://doi.org/10.1016/j.rse.2011.01.001)
- `MODIS_20X_1M_YYYY_V1` [(YYYY from 2000 to 2020; Yuan et al., 2011)](https://doi.org/10.1016/j.rse.2011.01.001)
- `MODIS_20X_8D_YYYY_V1` [(YYYY from 2000 to 2020; Yuan et al., 2011)](https://doi.org/10.1016/j.rse.2011.01.001)
</summary>

```
@article{yuan2011reprocessing,
	author = {Yuan, Hua and Dai, Yongjiu and Xiao, Zhiqiang and Ji, Duoying and Shangguan, Wei},
	year = {2011},
	title = {Reprocessing the MODIS Leaf Area Index products for land surface and climate modelling},
	journal = {Remote Sensing of Environment},
	volume = {115},
	number = {5},
	pages = {1171--1187}
}
```
</details>
"""
lai_collection() = (
    _supported = [];
    for _year in 2000:2020
        push!(_supported, "MODIS_2X_1M_$(_year)_V1");
        push!(_supported, "MODIS_2X_8D_$(_year)_V1");
        push!(_supported, "MODIS_10X_1M_$(_year)_V1");
        push!(_supported, "MODIS_10X_8D_$(_year)_V1");
        push!(_supported, "MODIS_20X_1M_$(_year)_V1");
        push!(_supported, "MODIS_20X_8D_$(_year)_V1");
    end;

    return GriddedCollection("LAI", _supported, "MODIS_2X_8D_2020_V1")
);

"""
    land_mask_collection()

Method to create a general dataset collection for land mask. Supported datasets are (click to view bibtex items)
- `4X_1Y_V1` [(ERA5)]
"""
land_mask_collection() = GriddedCollection("LM", ["4X_1Y_V1"], "4X_1Y_V1");

"""
    leaf_chlorophyll_collection()

<details>
<summary>
Method to create a general dataset collection for leaf chlorophyll content. Supported datasets are (click to view bibtex items)
- `2X_7D_V1` [(Croft et al., 2017)](https://doi.org/10.1016/j.rse.2019.111479)
</summary>

```
@article{croft2020global,
    author = {Croft, H and Chen, JM and Wang, R and Mo, G and Luo, S and Luo, X and He, L and Gonsamo, A and Arabian, J and Zhang, Y and others},
    year = {2020},
    title = {The global distribution of leaf chlorophyll content},
    journal = {Remote Sensing of Environment},
    volume = {236},
    pages = {111479},
}
```
</details>
"""
leaf_chlorophyll_collection() = (
    @warn "This dataset is only meant for those who has reached to the authors (Croft et al) for permissions. We (developers of GriddingMachine) are not responsible for unauthorized usage";

    return GriddedCollection("CHL", ["2X_7D_V1"], "2X_7D_V1")
);

"""
    leaf_nitrogen_collection()

<details>
<summary>
Method to create a general dataset collection for leaf nitrogen content. Supported datasets are (click to view bibtex items)
- `2X_1Y_V1` [(Butler et al., 2017)](https://doi.org/10.1073/pnas.1708984114)
- `2X_1Y_V2` [(Boonman et al., 2020)](https://doi.org/10.1111/geb.13086)
</summary>

```
@article{butler2017mapping,
    author = {Butler, Ethan E and Datta, Abhirup and Flores-Moreno, Habacuc and Chen, Ming and Wythers, Kirk R and Fazayeli, Farideh and Banerjee, Arindam and Atkin, Owen K and Kattge, Jens and
              Amiaud, Bernard and others},
    year = {2017},
    title = {Mapping local and global variability in plant trait distributions},
    journal = {Proceedings of the National Academy of Sciences},
    volume = {114},
    number = {51},
    pages = {E10937--E10946}
}
@article{boonman2020assessing,
    author = {Boonman, Coline CF and Ben{\\'i}tez-L{\\'o}pez, Ana and Schipper, Aafke M and Thuiller, Wilfried and Anand, Madhur and Cerabolini, Bruno EL and Cornelissen, Johannes HC and
              Gonzalez-Melo, Andres and Hattingh, Wesley N and Higuchi, Pedro and others},
    year = {2020},
    title = {Assessing the reliability of predicted plant trait distributions at the global scale},
    journal = {Global Ecology and Biogeography},
    volume = {29},
    number = {6},
    pages = {1034--1051}
}
```
</details>
"""
leaf_nitrogen_collection() = GriddedCollection("LNC", ["2X_1Y_V1", "2X_1Y_V2"], "2X_1Y_V1");

"""
    leaf_phosphorus_collection()

<details>
<summary>
Method to create a general dataset collection for leaf phosphorus content. Supported datasets are (click to view bibtex items)
- `2X_1Y_V1` [(Butler et al., 2017)](https://doi.org/10.1073/pnas.1708984114)
</summary>

```
@article{butler2017mapping,
    author = {Butler, Ethan E and Datta, Abhirup and Flores-Moreno, Habacuc and Chen, Ming and Wythers, Kirk R and Fazayeli, Farideh and Banerjee, Arindam and Atkin, Owen K and Kattge, Jens and
              Amiaud, Bernard and others},
    year = {2017},
    title = {Mapping local and global variability in plant trait distributions},
    journal = {Proceedings of the National Academy of Sciences},
    volume = {114},
    number = {51},
    pages = {E10937--E10946}
}
```
</details>
"""
leaf_phosphorus_collection() = GriddedCollection("LPC", ["2X_1Y_V1"], "2X_1Y_V1");

"""
    pft_collection()

<details>
<summary>
Method to create a general dataset collection for plant function type ratio. Supported datasets are (click to view bibtex items)
- `2X_1Y_V1` [(Lawrence and Chase, 2007)](https://doi.org/10.1029/2006JG000168)
</summary>

```
@article{lawrence2007representing,
	author = {Lawrence, Peter J and Chase, Thomas N},
	year = {2007},
	title = {Representing a new MODIS consistent land surface in the Community Land Model (CLM 3.0)},
	journal = {Journal of Geophysical Research: Biogeosciences},
	volume = {112},
	pages = {G01023}
}
```
</details>
"""
pft_collection() = GriddedCollection("PFT", ["2X_1Y_V1"], "2X_1Y_V1");

"""
    sif_collection()

<details>
<summary>
Method to create a general dataset collection for solar-induced chlorophyll fluorescence. Supported datasets are (click to view bibtex items)
- `TROPOMI_740_1X_1M_YYYY_V1` [(YYYY from 2018 to 2019; Köhler et al., 2018)](https://doi.org/10.1029/2018GL079031)
- `TROPOMI_740_12X_8D_YYYY_V1` [(YYYY from 2018 to 2019; Köhler et al., 2018)](https://doi.org/10.1029/2018GL079031)
- `TROPOMI_740DC_1X_1M_YYYY_V1` [(YYYY from 2018 to 2019; Köhler et al., 2018)](https://doi.org/10.1029/2018GL079031)
- `TROPOMI_740DC_12X_8D_YYYY_V1` [(YYYY from 2018 to 2019; Köhler et al., 2018)](https://doi.org/10.1029/2018GL079031)
</summary>

```
@article{kohler2018global,
    author = {K{\\"o}hler, Philipp and Frankenberg, Christian and Magney, Troy S and Guanter, Luis and Joiner, Joanna and Landgraf, Jochen},
    year = {2018},
    title = {Global retrievals of solar-induced chlorophyll fluorescence with {TROPOMI}: {F}irst results and intersensor comparison to {OCO-2}},
    journal = {Geophysical Research Letters},
    volume = {45},
    number = {19},
    pages = {10,456--10,463}
}
```
</details>
"""
sif_collection() = (
    _supported = [];
    for _year in 2018:2019
        push!(_supported, "TROPOMI_740_1X_1M_$(_year)_V1");
        push!(_supported, "TROPOMI_740_12X_8D_$(_year)_V1");
        push!(_supported, "TROPOMI_740DC_1X_1M_$(_year)_V1");
        push!(_supported, "TROPOMI_740DC_12X_8D_$(_year)_V1");
    end;

    return GriddedCollection("SIF", _supported, "TROPOMI_740_1X_1M_2019_V1")
);

"""
    sla_collection()

<details>
<summary>
Method to create a general dataset collection for SLA (specific leaf area). Supported datasets are (click to view bibtex items)
- `2X_1Y_V1` [(Butler et al., 2017)](https://doi.org/10.1073/pnas.1708984114)
- `2X_1Y_V2` [(Boonman et al., 2020)](https://doi.org/10.1111/geb.13086)
</summary>

```
@article{butler2017mapping,
    author = {Butler, Ethan E and Datta, Abhirup and Flores-Moreno, Habacuc and Chen, Ming and Wythers, Kirk R and Fazayeli, Farideh and Banerjee, Arindam and Atkin, Owen K and Kattge, Jens and
              Amiaud, Bernard and others},
    year = {2017},
    title = {Mapping local and global variability in plant trait distributions},
    journal = {Proceedings of the National Academy of Sciences},
    volume = {114},
    number = {51},
    pages = {E10937--E10946}
}
@article{boonman2020assessing,
    author = {Boonman, Coline CF and Ben{\\'i}tez-L{\\'o}pez, Ana and Schipper, Aafke M and Thuiller, Wilfried and Anand, Madhur and Cerabolini, Bruno EL and Cornelissen, Johannes HC and
              Gonzalez-Melo, Andres and Hattingh, Wesley N and Higuchi, Pedro and others},
    year = {2020},
    title = {Assessing the reliability of predicted plant trait distributions at the global scale},
    journal = {Global Ecology and Biogeography},
    volume = {29},
    number = {6},
    pages = {1034--1051}
}
```
</details>
"""
sla_collection() = GriddedCollection("SLA", ["2X_1Y_V1", "2X_1Y_V2"], "2X_1Y_V1");

"""
    soil_color_collection()

<details>
<summary>
Method to create a general dataset collection for soil color class to use with soil albedo. Supported datasets are (click to view bibtex items)
- `2X_1Y_V1` [(Lawrence and Chase, 2007)](https://doi.org/10.1029/2006JG000168)
</summary>

```
@article{lawrence2007representing,
    author = {Lawrence, Peter J and Chase, Thomas N},
    year = {2007},
    title = {Representing a new MODIS consistent land surface in the Community Land Model (CLM 3.0)},
    journal = {Journal of Geophysical Research: Biogeosciences},
    volume = {112},
    pages = {G01023}
}
```
</details>
"""
soil_color_collection() = GriddedCollection("SC", ["2X_1Y_V1"], "2X_1Y_V1");

"""
    soil_hydraulics_collection()

<details>
<summary>
Method to create a general dataset collection for soil hydraulic parameters (residual soil water content - SWCR, saturated soil water content - SWCS, van Genuchten α - VGA, van Genuchten n - VGN).
    Supported datasets are (click to view bibtex items)
- `SWCR_120X_1Y_V1` [(Dai et al., 2019)](https://doi.org/10.1029/2019MS001784)
- `SWCR_12X_1Y_V1` [(regridded; Dai et al., 2019)](https://doi.org/10.1029/2019MS001784)
- `SWCS_120X_1Y_V1` [(Dai et al., 2019)](https://doi.org/10.1029/2019MS001784)
- `SWCS_12X_1Y_V1` [(regridded; Dai et al., 2019)](https://doi.org/10.1029/2019MS001784)
- `VGA_120X_1Y_V1` [(Dai et al., 2019)](https://doi.org/10.1029/2019MS001784)
- `VGA_12X_1Y_V1` [(regridded; Dai et al., 2019)](https://doi.org/10.1029/2019MS001784)
- `VGN_120X_1Y_V1` [(Dai et al., 2019)](https://doi.org/10.1029/2019MS001784)
- `VGN_12X_1Y_V1` [(regridded; Dai et al., 2019)](https://doi.org/10.1029/2019MS001784)
</summary>

```
@article{dai2019global,
    author = {Dai, Yongjiu and Xin, Qinchuan and Wei, Nan and Zhang, Yonggen and Shangguan, Wei and Yuan, Hua and Zhang, Shupeng and Liu, Shaofeng and Lu, Xingjie},
	year = {2019},
    title = {A global high-resolution data set of soil hydraulic and thermal properties for land surface modeling},
	journal = {Journal of Advances in Modeling Earth Systems},
	volume = {11},
	number = {9},
	pages = {2996--3023}
}
```
</details>
"""
soil_hydraulics_collection() = (
    _supported = ["SWCR_120X_1Y_V1", "SWCR_12X_1Y_V1", "SWCS_120X_1Y_V1", "SWCS_12X_1Y_V1", "VGA_120X_1Y_V1", "VGA_12X_1Y_V1", "VGN_120X_1Y_V1", "VGN_12X_1Y_V1"];

    return GriddedCollection("SOIL", _supported, "SWCS_12X_1Y_V1")
);

"""
    surface_area_collection()

<details>
<summary>
Method to create a general dataset collection for earth surface area. Supported datasets are (click to view bibtex items)
- `2X_1Y_V1` [(Lawrence and Chase, 2007)](https://doi.org/10.1029/2006JG000168)
- `2X_1Y_V1` [(regridded; Lawrence and Chase, 2007)](https://doi.org/10.1029/2006JG000168)
</summary>

```
@article{lawrence2007representing,
author = {Lawrence, Peter J and Chase, Thomas N},
year = {2007},
title = {Representing a new MODIS consistent land surface in the Community Land Model (CLM 3.0)},
journal = {Journal of Geophysical Research: Biogeosciences},
volume = {112},
pages = {G01023}
}
```
</details>
"""
surface_area_collection() = GriddedCollection("SA", ["2X_1Y_V1", "1X_1Y_V1"], "2X_1Y_V1");

"""
    tree_density_collection()

<details>
<summary>
Method to create a general dataset collection for tree density (number of trees per area). Supported datasets are (click to view bibtex items)
- `120X_1Y_V1` [(Crowther et al., 2017)](https://doi.org/10.1038/nature14967)
- `2X_1Y_V1` [(regridded; Crowther et al., 2020)](https://doi.org/10.1038/nature14967)
</summary>

```
@article{crowther2015mapping,
    author = {Crowther, Thomas W and Glick, Henry B and Covey, Kristofer R and Bettigole, Charlie and Maynard, Daniel S and Thomas, Stephen M and Smith, Jeffrey R and Hintler, Gregor and
              Duguid, Marlyse C and Amatulli, Giuseppe and others},
    year = {2015},
    title = {Mapping tree density at a global scale},
    journal = {Nature},
    volume = {525},
    number = {7568},
    pages = {201--205}
}
```
</details>
"""
tree_density_collection() = GriddedCollection("TD", ["120X_1Y_V1", "2X_1Y_V1"], "2X_1Y_V1");

"""
    vcmax_collection()

<details>
<summary>
Method to create a general dataset collection for Vcmax. Supported datasets are (click to view bibtex items)
- `2X_1Y_V1` [(Smith et al., 2019)](https://doi.org/10.1111/ele.13210)
- `2X_1Y_V2` [(Luo et al., 2019)](https://doi.org/10.1038/s41467-021-25163-9)
</summary>

```
@article{smith2019global,
    author = {Smith, Nicholas G. and Keenan, Trevor F. and Prentice, I. Colin and Wang, Han and Wright, Ian J. and Niinemets, Ülo and Crous, Kristine Y. and Domingues, Tomas F. and
              Guerrieri, Rossella and {Yoko Ishida}, F. and Zhou, Shuangxi},
    year = {2019},
    title = {Global photosynthetic capacity is optimized to the environment},
    journal = {Ecology Letters},
    volume = {22},
    number = {3},
    pages = {506–517}
}
@article{luo2021global,
	author = {Luo, Xiangzhong and Keenan, Trevor F. and Chen, Jing M. and Croft, Holly and {Colin Prentice}, I. and Smith, Nicholas G. and Walker, Anthony P. and Wang, Han and Wang, Rong and
              Xu, Chonggang and Zhang, Yao},
	year = {2021},
	title = {Global variation in the fraction of leaf nitrogen allocated to photosynthesis},
	journal = {Nature Communications},
	volume = {12},
	number = {1},
	pages = {4866}
}
```
</details>
"""
vcmax_collection() = GriddedCollection("VCMAX", ["2X_1Y_V1", "2X_1Y_V2"], "2X_1Y_V2");

"""
    wood_density_collection()

<details>
<summary>
Method to create a general dataset collection for wood density. Supported datasets are (click to view bibtex items)
- `2X_1Y_V1` [(Boonman et al., 2020)](https://doi.org/10.1111/geb.13086)
</summary>

```
@article{boonman2020assessing,
    author = {Boonman, Coline CF and Ben{\\'i}tez-L{\\'o}pez, Ana and Schipper, Aafke M and Thuiller, Wilfried and Anand, Madhur and Cerabolini, Bruno EL and Cornelissen, Johannes HC and
              Gonzalez-Melo, Andres and Hattingh, Wesley N and Higuchi, Pedro and others},
    year = {2020},
    title = {Assessing the reliability of predicted plant trait distributions at the global scale},
    journal = {Global Ecology and Biogeography},
    volume = {29},
    number = {6},
    pages = {1034--1051}
}
```
</details>
"""
wood_density_collection() = GriddedCollection("WD", ["2X_1Y_V1"], "2X_1Y_V1");

# query file from gridded collections
"""
This function queries data path for a dataset, supported methods are

$(METHODLIST)

"""
function query_collection end

"""
    query_collection(ds::GriddedCollection)

This method queries the local data path from collection for the default data, given
- `ds` [`GriddedCollection`](@ref) type collection

---
# Examples
```julia
dat_file = query_collection(canopy_height_collection());
```
"""
query_collection(ds::GriddedCollection) = query_collection(ds, ds.DEFAULT_COMBO);

"""
    query_collection(ds::GriddedCollection, version::String)

This method queries the local data path from collection, given
- `ds` [`GriddedCollection`](@ref) type collection
- `version` Queried dataset version (must be in `ds.SUPPORTED_COMBOS`)

---
# Examples
```julia
dat_file = query_collection(canopy_height_collection(), "20X_1Y_V1");
```
"""
query_collection(ds::GriddedCollection, version::String) = (
    # make sure requested version is in the
    @assert version in ds.SUPPORTED_COMBOS;

    # determine file name from label and supported version
    _fn = "$(ds.LABEL)_$(version)";

    return query_collection(_fn)
);

"""
    query_collection(artname::String)

This method queries the local data path from given artifact name
- `artname` Artifact name
"""
query_collection(artname::String) = (
    _metas = load_artifacts_toml(joinpath(@__DIR__, "../Artifacts.toml"));
    _artns = [_name for (_name,_) in _metas];
    @assert artname in _artns;

    return @artifact_str(artname) * "/$(artname).nc"
);

"""
This function cleans up the collections, supported methods are

    $(METHODLIST)
"""
function clean_collections! end

"""
    clean_collections!(selection::String="old")

This method cleans up all selected artifacts of GriddingMachine.jl (through identify the `GRIDDINGMACHINE` file in the artifacts), given
- `selection` A string indicating which artifacts to clean up
    - `old` Artifacts from an old version of GriddingMachine.jl (default)
    - `all` All Artifacts from GriddingMachine.jl

---
# Examples
```julia
clean_collections!();
clean_collections!("old");
clean_collections!("all");
```
"""
clean_collections!(selection::String="old") = (
    # read the SHA1 identifications in Artifacts.toml
    _metas = load_artifacts_toml(joinpath(@__DIR__, "../Artifacts.toml"));
    _hashs = [_meta["git-tree-sha1"] for (_,_meta) in _metas];

    # iterate through the artifacts and remove the old one that is not in current Artifacts.toml or remove all artifacts within GriddingMachine.jl
    _artifact_dirs = readdir("$(homedir())/.julia/artifacts");
    for _dir in _artifact_dirs
        if isdir("$(homedir())/.julia/artifacts/$(_dir)")
            if isfile("$(homedir())/.julia/artifacts/$(_dir)/GRIDDINGMACHINE")
                if selection == "all"
                    rm("$(homedir())/.julia/artifacts/$(_dir)"; recursive=true, force=true);
                else
                    if !(_dir in _hashs)
                        rm("$(homedir())/.julia/artifacts/$(_dir)"; recursive=true, force=true);
                    end;
                end;
            end;
        end;
    end;

    return nothing
);

"""
    clean_collections!(selection::Vector{String})

This method cleans up all selected artifacts in GriddingMachine.jl, given
- `selection` A vector of artifact names

---
# Examples
```julia
clean_collections!(["PFT_2X_1Y_V1"]);
```
"""
clean_collections!(selection::Vector{String}) = (
    # read the SHA1 identifications in Artifacts.toml
    _metas = load_artifacts_toml(joinpath(@__DIR__, "../Artifacts.toml"));
    _hashs = [_metas[_artn]["git-tree-sha1"] for _artn in selection];

    # iterate the artifact hashs to remove corresponding folder
    for _dir in _hashs
        rm("$(homedir())/.julia/artifacts/$(_dir)"; recursive=true, force=true);
    end;

    return nothing
);

"""
    clean_collections!(selection::GriddedCollection)

This method cleans up all selected artifacts in GriddingMachine.jl, given
- `selection` A [`GriddedCollection`](@ref) type collection

---
# Examples
```julia
clean_collections!(pft_collection());
```
"""
clean_collections!(selection::GriddedCollection) = (
    clean_collections!(["$(selection.LABEL)_$(_ver)" for _ver in selection.SUPPORTED_COMBOS]);

    return nothing
);

end


text:
using Test
using Foo

@test 1 == 1