text:
"""
    initialize_flextopo_model(config::Config)

Initial part of the FlexTopo model concept. Reads the input settings and data as defined in the
Config object. Will return a Model that is ready to run.
"""
function initialize_flextopo_model(config::Config)
    # unpack the paths to the NetCDF files
    static_path = input_path(config, config.input.path_static)

    reader = prepare_reader(config)
    clock = Clock(config, reader)
    Δt = clock.Δt

    reinit = get(config.model, "reinit", true)::Bool
    do_reservoirs = get(config.model, "reservoirs", false)::Bool
    do_lakes = get(config.model, "lakes", false)::Bool
    do_pits = get(config.model, "pits", false)::Bool

    kw_river_tstep = get(config.model, "kw_river_tstep", 0)
    kw_land_tstep = get(config.model, "kw_land_tstep", 0)
    kinwave_it = get(config.model, "kin_wave_iteration", false)::Bool

    classes = get(config.model, "classes", "")
    nclass = length(classes)
    @info "Classes are set to names `$(join(classes,", "))` and have size `$nclass`."
    kclass = [1]    #needed to initialize

    # dictionary of available functions for each store
    dic_function = Dict{String,Function}(
        "common_snow_hbv" => common_snow_hbv,
        "common_snow_no_storage" => common_snow_no_storage,
        "common_glaciers" => common_glaciers,
        "interception_overflow" => interception_overflow,
        "interception_no_storage" => interception_no_storage,
        "hortonponding" => hortonponding,
        "hortonponding_no_storage" => hortonponding_no_storage,
        "hortonrunoff" => hortonrunoff,
        "hortonrunoff_no_storage" => hortonrunoff_no_storage,
        "rootzone_storage" => rootzone_storage,
        "rootzone_no_storage" => rootzone_no_storage,
        "fast_no_storage" => fast_no_storage,
        "fast_storage" => fast_storage,
        "slow_no_storage" => slow_no_storage,
        "common_slow_storage" => common_slow_storage,
    )

    select_snow = get(config.model, "select_snow", ["common_snow_hbv"])
    select_interception = get(config.model, "select_interception", ["interception_overflow"])
    select_hortonponding = get(config.model, "select_hortonponding", ["hortonponding_no_storage"])
    select_hortonrunoff = get(config.model, "select_hortonrunoff", ["hortonrunoff_no_storage"])
    select_rootzone = get(config.model, "select_rootzone", ["rootzone_storage"])
    select_fast = get(config.model, "select_fast", ["fast_storage"])
    select_slow = get(config.model, "select_slow", ["common_slow_storage"])

    nc = NCDataset(static_path)

    subcatch_2d =
        ncread(nc, config, "subcatchment"; optional = false, allow_missing = true)
    # indices based on catchment
    inds, rev_inds = active_indices(subcatch_2d, missing)
    n = length(inds)

    # glacier parameters
    g_tt = ncread(
        nc,
        config,
        "vertical.g_tt";
        sel = inds,
        defaults = 0.0,
        type = Float,
        fill = 0.0,
    )
    g_cfmax =
        ncread(
            nc,
            config,
            "vertical.g_cfmax";
            sel = inds,
            defaults = 3.0,
            type = Float,
            fill = 0.0,
        ) .* (Δt / basetimestep)

    g_sifrac = ncread(
        nc,
        config,
        "vertical.g_sifrac";
        sel = inds,
        defaults = 0.001,
        type = Float,
        fill = 0.0,
    )
    glacierfrac = ncread(
        nc,
        config,
        "vertical.glacierfrac";
        sel = inds,
        defaults = 0.0,
        type = Float,
        fill = 0.0,
    )
    glacierstore = ncread(
        nc,
        config,
        "vertical.glacierstore";
        sel = inds,
        defaults = 5500.0,
        type = Float,
        fill = 0.0,
    )

    #snow param (single class)
    cfmax =
        ncread(
            nc,
            config,
            "vertical.cfmax";
            sel = inds,
            defaults = 3.75653,
            type = Float,
        ) .* (Δt / basetimestep)

    tt = ncread(
        nc,
        config,
        "vertical.tt";
        sel = inds,
        defaults = -1.41934,
        type = Float,
    )

    tti = ncread(nc, config, "vertical.tti"; sel = inds, defaults = 1.0, type = Float)

    ttm = ncread(
        nc,
        config,
        "vertical.ttm";
        sel = inds,
        defaults = -1.41934,
        type = Float,
    )

    whc = ncread(nc, config, "vertical.whc"; sel = inds, defaults = 0.1, type = Float)

    cfr =
        ncread(nc, config, "vertical.cfr"; sel = inds, defaults = 0.05, type = Float)

    #parameters which are not class specific
    pcorr =
        ncread(nc, config, "vertical.pcorr"; sel = inds, defaults = 1.0, type = Float)
    ecorr =
        ncread(nc, config, "vertical.ecorr"; sel = inds, defaults = 1.0, type = Float)
    rfcf =
        ncread(nc, config, "vertical.rfcf"; sel = inds, defaults = 1.0, type = Float)
    sfcf =
        ncread(nc, config, "vertical.sfcf"; sel = inds, defaults = 1.0, type = Float)
    ks =
        ncread(
            nc,
            config,
            "vertical.ks";
            sel = inds,
            defaults = 0.006,
            type = Float,
        ) .* (Δt / basetimestep)

    #initialize parameters that differ per class
    hrufrac = ncread(
        nc,
        config,
        "vertical.hrufrac";
        sel = inds,
        defaults = 1.0 / length(classes),
        type = Float,
        dimname = :classes,
    )
    imax = ncread(
        nc,
        config,
        "vertical.imax";
        sel = inds,
        defaults = 3.0,
        type = Float,
        dimname = :classes,
    )
    shmax = ncread(
        nc,
        config,
        "vertical.shmax";
        sel = inds,
        defaults = 30.0,
        type = Float,
        dimname = :classes,
    )
    khf =
        ncread(
            nc,
            config,
            "vertical.khf";
            sel = inds,
            defaults = 0.5,
            type = Float,
            dimname = :classes,
        ) .* (Δt / basetimestep)
    facc0 = ncread(
        nc,
        config,
        "vertical.facc0";
        sel = inds,
        defaults = -3.0,
        type = Float,
        dimname = :classes,
    )
    facc1 = ncread(
        nc,
        config,
        "vertical.facc1";
        sel = inds,
        defaults = 0.0,
        type = Float,
        dimname = :classes,
    )
    fdec = ncread(
        nc,
        config,
        "vertical.fdec";
        sel = inds,
        defaults = 0.2,
        type = Float,
        dimname = :classes,
    )
    fmax =
        ncread(
            nc,
            config,
            "vertical.fmax";
            sel = inds,
            defaults = 2.0,
            type = Float,
            dimname = :classes,
        ) .* (Δt / basetimestep)
    shmin = ncread(
        nc,
        config,
        "vertical.shmin";
        sel = inds,
        defaults = 0.2,
        type = Float,
        dimname = :classes,
    )
    kmf = ncread(
        nc,
        config,
        "vertical.kmf";
        sel = inds,
        defaults = 1.0,
        type = Float,
        dimname = :classes,
    )
    srmax = ncread(
        nc,
        config,
        "vertical.srmax";
        sel = inds,
        defaults = 260.0,
        type = Float,
        dimname = :classes,
    )
    beta = ncread(
        nc,
        config,
        "vertical.beta";
        sel = inds,
        defaults = 0.3,
        type = Float,
        dimname = :classes,
    )
    lp = ncread(
        nc,
        config,
        "vertical.lp";
        sel = inds,
        defaults = 0.3,
        type = Float,
        dimname = :classes,
    )
    perc =
        ncread(
            nc,
            config,
            "vertical.perc";
            sel = inds,
            defaults = 0.30,
            type = Float,
            dimname = :classes,
        ) .* (Δt / basetimestep)
    cap =
        ncread(
            nc,
            config,
            "vertical.cap";
            sel = inds,
            defaults = 0.20,
            type = Float,
            dimname = :classes,
        ) .* (Δt / basetimestep)
    kf =
        ncread(
            nc,
            config,
            "vertical.kf";
            sel = inds,
            defaults = 0.1,
            type = Float,
            dimname = :classes,
        ) .* (Δt / basetimestep)
    alfa = ncread(
        nc,
        config,
        "vertical.alfa";
        sel = inds,
        defaults = 1.0,
        type = Float,
        dimname = :classes,
    )
    ds = ncread(
        nc,
        config,
        "vertical.ds";
        sel = inds,
        defaults = 0.2,
        type = Float,
        dimname = :classes,
    )

    # read x, y coordinates and calculate cell length [m]
    y_nc = read_y_axis(nc)
    x_nc = read_x_axis(nc)
    y = permutedims(repeat(y_nc, outer = (1, length(x_nc))))[inds]
    cellength = abs(mean(diff(x_nc)))

    sizeinmetres = get(config.model, "sizeinmetres", false)::Bool
    xl, yl = cell_lengths(y, cellength, sizeinmetres)

    interceptionstorage = zeros(Float, nclass, n)
    hortonpondingstorage = zeros(Float, nclass, n)
    hortonrunoffstorage = zeros(Float, nclass, n)
    srootzone_over_srmax = zeros(Float, nclass, n) .+ 1.0

    states_ = fill(mv, nclass, n)

    potsoilevap = fill(mv, nclass, n)
    intevap = fill(mv, nclass, n)
    precipeffective = fill(mv, nclass, n)
    hortonevap = fill(mv, nclass, n)
    rootevap = fill(mv, nclass, n)
    qhortonpond = fill(mv, nclass, n)
    qhortonrootzone = fill(mv, nclass, n)
    facc = zeros(Float, nclass, n)
    qhortonrun = fill(mv, nclass, n)
    qrootzone = fill(mv, nclass, n)
    qrootzonefast = fill(mv, nclass, n)
    qfast = fill(mv, nclass, n)
    actevap = fill(mv, nclass, n)
    qpercolation = fill(mv, nclass, n)
    qcapillary = fill(mv, nclass, n)

    wb_interception = fill(mv, nclass, n)
    wb_hortonponding = fill(mv, nclass, n)
    wb_hortonrunoff = fill(mv, nclass, n)
    wb_rootzone = fill(mv, nclass, n)
    wb_fast = fill(mv, nclass, n)

    flextopo = FLEXTOPO{Float,nclass}(
        Δt = Float(tosecond(Δt)),
        nclass = nclass,
        n = n,
        dic_function = dic_function,
        kclass = kclass,
        classes = classes,
        select_snow = select_snow,
        select_interception = select_interception,
        select_hortonponding = select_hortonponding,
        select_hortonrunoff = select_hortonrunoff,
        select_rootzone = select_rootzone,
        select_fast = select_fast,
        select_slow = select_slow,
        hrufrac = svectorscopy(hrufrac, Val{nclass}()),
        pcorr = pcorr,
        ecorr = ecorr,
        #glaciers
        g_tt = g_tt,
        g_cfmax = g_cfmax,
        g_sifrac = g_sifrac,
        glacierfrac = glacierfrac,
        glacierstore = glacierstore,
        # snow
        tti = tti,
        tt = tt,
        ttm = ttm,
        cfmax = cfmax,
        whc = whc,
        cfr = cfr,
        rfcf = rfcf,
        sfcf = sfcf,
        #interception
        imax = svectorscopy(imax, Val{nclass}()),
        #horton
        shmax = svectorscopy(shmax, Val{nclass}()),
        khf = svectorscopy(khf, Val{nclass}()),
        facc0 = svectorscopy(facc0, Val{nclass}()),
        facc1 = svectorscopy(facc1, Val{nclass}()),
        fdec = svectorscopy(fdec, Val{nclass}()),
        fmax = svectorscopy(fmax, Val{nclass}()),
        shmin = svectorscopy(shmin, Val{nclass}()),
        kmf = svectorscopy(kmf, Val{nclass}()),
        #root zone
        srmax = svectorscopy(srmax, Val{nclass}()),
        lp = svectorscopy(lp, Val{nclass}()),
        beta = svectorscopy(beta, Val{nclass}()),
        perc = svectorscopy(perc, Val{nclass}()),
        cap = svectorscopy(cap, Val{nclass}()),
        #fast
        ds = svectorscopy(ds, Val{nclass}()),
        alfa = svectorscopy(alfa, Val{nclass}()),
        kf = svectorscopy(kf, Val{nclass}()),
        #slow
        ks = ks,

        # # default (cold) states:
        snow = zeros(Float, n),
        snowwater = zeros(Float, n),
        interceptionstorage = svectorscopy(interceptionstorage, Val{nclass}()),
        hortonpondingstorage = svectorscopy(hortonpondingstorage, Val{nclass}()),
        hortonrunoffstorage = svectorscopy(hortonrunoffstorage, Val{nclass}()),
        rootzonestorage = svectorscopy(srmax, Val{nclass}()),
        faststorage = 0.0 .* svectorscopy(srmax, Val{nclass}()),
        srootzone_over_srmax = svectorscopy(srootzone_over_srmax, Val{nclass}()),
        slowstorage = zeros(Float, n) .+ 30.0,
        #states previous time step
        states_m = fill(mv, n),
        states_ = svectorscopy(states_, Val{nclass}()),
        #states averaged over all classes
        interceptionstorage_m = zeros(Float, n),
        hortonpondingstorage_m = zeros(Float, n),
        hortonrunoffstorage_m = zeros(Float, n),
        srootzone_m = zeros(Float, n),
        faststorage_m = zeros(Float, n),
        srootzone_over_srmax_m = zeros(Float, n),

        # variables:
        precipitation = fill(mv, n),
        temperature = fill(mv, n),
        potential_evaporation = fill(mv, n),
        epotcorr = fill(mv, n),
        precipcorr = fill(mv, n),
        rainfallplusmelt = fill(mv, n),
        snowfall = fill(mv, n),
        snowmelt = fill(mv, n),
        potsoilevap = svectorscopy(potsoilevap, Val{nclass}()),
        precipeffective = svectorscopy(precipeffective, Val{nclass}()),
        intevap = svectorscopy(intevap, Val{nclass}()),
        hortonevap = svectorscopy(hortonevap, Val{nclass}()),
        rootevap = svectorscopy(rootevap, Val{nclass}()),
        qhortonpond = svectorscopy(qhortonpond, Val{nclass}()),
        qhortonrootzone = svectorscopy(qhortonrootzone, Val{nclass}()),
        facc = svectorscopy(facc, Val{nclass}()),
        qhortonrun = svectorscopy(qhortonrun, Val{nclass}()),
        qrootzone = svectorscopy(qrootzone, Val{nclass}()),
        qrootzonefast = svectorscopy(qrootzonefast, Val{nclass}()),
        qrootzoneslow_m = fill(mv, n),
        qfast = svectorscopy(qfast, Val{nclass}()),
        actevap = svectorscopy(actevap, Val{nclass}()),
        qpercolation = svectorscopy(qpercolation, Val{nclass}()),
        qcapillary = svectorscopy(qcapillary, Val{nclass}()),
        qpercolation_m = fill(mv, n),
        qcapillary_m = fill(mv, n),
        # combined for the classes
        actevap_m = fill(mv, n),
        intevap_m = fill(mv, n),
        hortonevap_m = fill(mv, n),
        rootevap_m = fill(mv, n),
        qslow = fill(mv, n),
        qfast_tot = fill(mv, n),
        runoff = fill(mv, n),
        wb_snow = fill(mv, n),
        wb_interception = svectorscopy(wb_interception, Val{nclass}()),
        wb_hortonponding = svectorscopy(wb_hortonponding, Val{nclass}()),
        wb_hortonrunoff = svectorscopy(wb_hortonrunoff, Val{nclass}()),
        wb_rootzone = svectorscopy(wb_rootzone, Val{nclass}()),
        wb_fast = svectorscopy(wb_fast, Val{nclass}()),
        wb_slow = fill(mv, n),
        wb_tot = fill(mv, n),
    )

    modelsize_2d = size(subcatch_2d)
    river_2d = ncread(
        nc,
        config,
        "river_location";
        optional = false,
        type = Bool,
        fill = false,
    )
    river = river_2d[inds]
    riverwidth_2d = ncread(
        nc,
        config,
        "lateral.river.width";
        optional = false,
        type = Float,
        fill = 0,
    )
    riverwidth = riverwidth_2d[inds]
    riverlength_2d = ncread(
        nc,
        config,
        "lateral.river.length";
        optional = false,
        type = Float,
        fill = 0,
    )
    riverlength = riverlength_2d[inds]

    inds_riv, rev_inds_riv = active_indices(river_2d, 0)
    nriv = length(inds_riv)

    # reservoirs
    pits = zeros(Bool, modelsize_2d)
    if do_reservoirs
        reservoirs, resindex, reservoir, pits =
            initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, tosecond(Δt))
    else
        reservoir = ()
        reservoirs = nothing
        resindex = fill(0, nriv)
    end

    # lakes
    if do_lakes
        lakes, lakeindex, lake, pits =
            initialize_natural_lake(config, nc, inds_riv, nriv, pits, tosecond(Δt))
    else
        lake = ()
        lakes = nothing
        lakeindex = fill(0, nriv)
    end

    ldd_2d = ncread(nc, config, "ldd"; optional = false, allow_missing = true)
    ldd = ldd_2d[inds]
    if do_pits
        pits_2d =
            ncread(nc, config, "pits"; optional = false, type = Bool, fill = false)
        ldd = set_pit_ldd(pits_2d, ldd, inds)
    end

    βₗ = ncread(
        nc,
        config,
        "lateral.land.slope";
        optional = false,
        sel = inds,
        type = Float,
    )
    clamp!(βₗ, 0.00001, Inf)

    dl = map(detdrainlength, ldd, xl, yl)
    dw = (xl .* yl) ./ dl
    olf = initialize_surfaceflow_land(
        nc,
        config,
        inds;
        sl = βₗ,
        dl = dl,
        width = map(det_surfacewidth, dw, riverwidth, river),
        wb_pit = pits[inds],
        iterate = kinwave_it,
        tstep = kw_land_tstep,
        Δt = Δt,
    )

    graph = flowgraph(ldd, inds, pcr_dir)

    riverlength = riverlength_2d[inds_riv]
    riverwidth = riverwidth_2d[inds_riv]

    ldd_riv = ldd_2d[inds_riv]
    if do_pits
        ldd_riv = set_pit_ldd(pits_2d, ldd_riv, inds_riv)
    end
    graph_riv = flowgraph(ldd_riv, inds_riv, pcr_dir)

    # the indices of the river cells in the land(+river) cell vector
    index_river = filter(i -> !isequal(river[i], 0), 1:n)
    frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, βₗ, n)

    rf = initialize_surfaceflow_river(
        nc,
        config,
        inds_riv;
        dl = riverlength,
        width = riverwidth,
        wb_pit = pits[inds_riv],
        reservoir_index = resindex,
        reservoir = reservoirs,
        lake_index = lakeindex,
        lake = lakes,
        river = river,
        iterate = kinwave_it,
        tstep = kw_river_tstep,
        Δt = Δt,
    )

    # setup subdomains for the land and river kinematic wave domain, if nthreads = 1
    # subdomain is equal to the complete domain
    toposort = topological_sort_by_dfs(graph)
    toposort_riv = topological_sort_by_dfs(graph_riv)
    index_pit_land = findall(x -> x == 5, ldd)
    index_pit_river = findall(x -> x == 5, ldd_riv)
    subbas_order, indices_subbas, topo_subbas =
        kinwave_set_subdomains(config, graph, toposort, index_pit_land)
    subriv_order, indices_subriv, topo_subriv =
        kinwave_set_subdomains(config, graph_riv, toposort_riv, index_pit_river)

    modelmap = (vertical = flextopo, lateral = (land = olf, river = rf))
    indices_reverse = (
        land = rev_inds,
        river = rev_inds_riv,
        reservoir = isempty(reservoir) ? nothing : reservoir.reverse_indices,
        lake = isempty(lake) ? nothing : lake.reverse_indices,
    )

    # if fews_run = true, then classes are specified as Float64 index so FEWS can import
    # NetCDF output (3D data/SVector).
    fews_run = get(config, "fews_run", false)::Bool
    classes = fews_run ? Float64.(collect(1:length(flextopo.classes))) : flextopo.classes

    writer = prepare_writer(
        config,
        modelmap,
        indices_reverse,
        x_nc,
        y_nc,
        nc,
        extra_dim = (name = "classes", value = classes),
    )
    close(nc)

    # for each domain save:
    # - the directed acyclic graph (graph),
    # - the traversion order (order),
    # - upstream_nodes,
    # - subdomains for the kinematic wave domains for parallel execution (execution order of
    #   subbasins (subdomain_order), traversion order per subbasin (topo_subdomain) and
    #   Vector indices per subbasin matching the traversion order of the complete domain
    #   (indices_subdomain))
    # - the indices that map it back to the two dimensional grid (indices)

    # for the land domain the x and y length [m] of the grid cells are stored
    # for reservoirs and lakes indices information is available from the initialization
    # functions
    land = (
        graph = graph,
        upstream_nodes = filter_upsteam_nodes(graph, olf.wb_pit),
        subdomain_order = subbas_order,
        topo_subdomain = topo_subbas,
        indices_subdomain = indices_subbas,
        order = toposort,
        indices = inds,
        reverse_indices = rev_inds,
        xl = xl,
        yl = yl,
    )
    river = (
        graph = graph_riv,
        upstream_nodes = filter_upsteam_nodes(graph_riv, rf.wb_pit),
        subdomain_order = subriv_order,
        topo_subdomain = topo_subriv,
        indices_subdomain = indices_subriv,
        order = toposort_riv,
        indices = inds_riv,
        reverse_indices = rev_inds_riv,
    )

    model = Model(
        config,
        (; land, river, reservoir, lake, index_river, frac_toriver),
        (subsurface = nothing, land = olf, river = rf),
        flextopo,
        clock,
        reader,
        writer,
        FlextopoModel(),
    )

    # read and set states in model object if reinit=true
    if reinit == false
        instate_path = input_path(config, config.state.path_input)
        state_ncnames = ncnames(config.state)
        set_states(instate_path, model, state_ncnames; type = Float, dimname = :classes)
        # update kinematic wave volume for river and land domain
        @unpack lateral = model
        lateral.land.volume .= lateral.land.h .* lateral.land.width .* lateral.land.dl
        lateral.river.volume .= lateral.river.h .* lateral.river.width .* lateral.river.dl

        if do_lakes
            # storage must be re-initialized after loading the state with the current
            # waterlevel otherwise the storage will be based on the initial water level
            lakes.storage .=
                initialize_storage(lakes.storfunc, lakes.area, lakes.waterlevel, lakes.sh)
        end
    end

    return model
end

function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:FlextopoModel}
    @unpack lateral, vertical, network, clock, config = model

    inds_riv = network.index_river

    #COMMON SNOW
    vertical.dic_function[vertical.select_snow[1]](vertical)

    #lateral snow transport
    if get(config.model, "masswasting", false)::Bool
        lateral_snow_transport!(vertical.snow, vertical.snowwater, lateral.land.sl, network.land)
    end

    if get(config.model, "glacier", false)::Bool
        common_glaciers(vertical, config)
    end

    for (k, class) in enumerate(vertical.classes)
        vertical.kclass[1] = k

        #INTERCEPTION
        vertical.dic_function[vertical.select_interception[k]](vertical)

        #HORTON
        vertical.dic_function[vertical.select_hortonponding[k]](vertical)
        vertical.dic_function[vertical.select_hortonrunoff[k]](vertical)

        #ROOT-ZONE
        vertical.dic_function[vertical.select_rootzone[k]](vertical)

        #FAST
        vertical.dic_function[vertical.select_fast[k]](vertical)

    end

    #COMMON SLOW
    vertical.dic_function[vertical.select_slow[1]](vertical)

    # WAT BAL
    watbal(vertical)

    surface_routing(model)

    write_output(model)

    # update the clock
    advance!(clock)

    return model
end


text:
using MacroTools
using Compat, Compat.Test

@testset "MacroTools" begin

include("match.jl")
include("split.jl")
include("destruct.jl")

end

ext: jl


text:
using KnowedMimeTypes
using Base.Test

# write your own tests here
@test get_mime(".pdf") != Void

mt = MimeTypes(Dict("a" => "b"))

@test get_mime(mt, "a") == "b"


text:
function Base.summary(tech::TechnicalParams)
    return "$(typeof(tech))"
end

function Base.summary(sys::System)
    try
        return "System with $(sys.junctions) junctions and $(sys.arcs)"
    catch
        return "System (empty)"
    end
end

# function Base.show(io::IO, sys::System)
#     println(io, "basepower=$(sys.basepower)")
#     show(io, sys.data)
# end

# function Base.show(io::IO, ::MIME"text/plain", sys::System)
#     println(io, "System")
#     println(io, "======")
#     println(io, "Base Power: $(sys.basepower)\n")
#     show(io, MIME"text/plain"(), sys.data)
# end

# function Base.show(io::IO, ::MIME"text/html", sys::System)
#     println(io, "
System
")
#     println(io, "
Base Power: $(sys.basepower)
")
#     show(io, MIME"text/html"(), sys.data)
# end


text:
using IntelVectorMath
using ForwardDiff
using BenchmarkTools

function IVM.atan(xd::Vector{ForwardDiff.Dual{T, V, 1}},
                  yd::Vector{ForwardDiff.Dual{T, V, 1}}) where {T, V}

    x = ForwardDiff.value.(xd)
    y = ForwardDiff.value.(yd)
    z = IVM.atan(x, y)

    δx = ForwardDiff.partials.(xd, 1)
    δy = ForwardDiff.partials.(yd, 1)
    δz = @. (δy*x - δx*y)/(x^2 + y^2)

    ForwardDiff.Dual{T}.(z, δz)
end

struct fast_atan_dual_cache{T}
    x::Vector{T}
    y::Vector{T}
    z::Vector{T}
    δ::Matrix{T}
end

function allocate_fast_atan_dual_cache(T::Type, N::Int)
    x = Vector{T}(undef, N)
    y = Vector{T}(undef, N)
    z = Vector{T}(undef, N)
    δ = Matrix{T}(undef, N, 3)
    fast_atan_dual_cache(x, y, z, δ)
end

function IVM.atan!(zd::Vector{ForwardDiff.Dual{T, V, 1}},
                   xd::Vector{ForwardDiff.Dual{T, V, 1}},
                   yd::Vector{ForwardDiff.Dual{T, V, 1}},
                   cache::fast_atan_dual_cache{V}) where {T, V}

   x, y, z = cache.x, cache.y, cache.z
   δx, δy, δz = view(cache.δ,:,1), view(cache.δ,:,2), view(cache.δ,:,3)

    x .= ForwardDiff.value.(xd)
    y .= ForwardDiff.value.(yd)
    IVM.atan!(z, x, y)

    δx .= ForwardDiff.partials.(xd, 1)
    δy .= ForwardDiff.partials.(yd, 1)
    δz .= @. (δy*x - δx*y)/(x^2 + y^2)

    zd .= ForwardDiff.Dual{T}.(z, δz)
end

N = 500
x = rand(N)
y = rand(N)
δx = rand(N)
δy = rand(N)
xd = ForwardDiff.Dual{DERTag}.(x, δx)
yd = ForwardDiff.Dual{DERTag}.(y, δy)
zd = IVM.atan(xd,yd)

cache_atan = allocate_fast_atan_dual_cache(Float64, N)
@btime IVM.atan!(zd, xd, yd, cache_atan)
@btime $zd .= atan.($xd, $yd)

z = rand(N)
@btime IVM.atan!(z, x, y)
@btime z .= atan.(x, y)

struct DERTag end
v = rand(3)
p1 = rand(3)
p2 = rand(3)
p = [p1, p2]
d = ForwardDiff.Dual{DERTag}.(v, p...)

ForwardDiff.partials.(d, 1)
ForwardDiff.partials.(d, 2)
ForwardDiff.value.(d)


text:
@inline function calc_A_second_order(x::VInt16, ::Val{N}) where N
    p = 15; r = 1; A = Int16(23170); C = Int16(-54292 >> N)
    rounding = one(x) << (p - N - 1)
    (rounding + A + (x * (x * C) >> N) >> r) >> (p - N)
end

@inline function calc_B_second_order(x::VInt16, ::Val{N}) where N
    r = 14; B = Int16(-16384 >> N);
    rounding = one(x) << (r - N - 1)
    (rounding + x * B) >> (r - N)
end

@inline function calc_A(x::VInt16, bits::B) where B <: Union{Val{3}, Val{4}, Val{5}}
    calc_A_second_order(x, bits)
end
@inline function calc_B(x::VInt16, bits::B) where B <: Union{Val{3}, Val{4}, Val{5}}
    calc_B_second_order(x, bits)
end


text:
"""
Dynamic table data structure to store and update priorities

Implementation
Stores ids that identify the priorities. Each group stores a range of priority ids.
The basic design assumes a lower-most group of "zero" priorities, and a second
group storing all non-zero priorities that are < a `minpriority`. All other groups
store priorities within consecutive ranges. The `minpriority` is taken fixed at creation,
but the `maxpriority` will increase dynamically based on inserted / updated priorities.
The ranges are assumed to be powers of two:
   bin 1 = {0},
   bin 2 = (0,`minpriority`),
   bin 3 = [`minpriority`,`2*minpriority`)...
   bin N = [`.5*maxpriority`,`maxpriority`)
*Assumes* the `priortogid` function that maps priorities to groups maps the upper end of
the interval to the next group. i.e. maxpriority -> N+1
"""

"""
One group (i.e. bin) of priority ids within the table
"""
mutable struct PriorityGroup{T,W <: AbstractVector}
    "(strict) upper bound for priorities in this group"
    maxpriority::T

    "number of priority ids"
    numpids::Int

    "priority ids associated with group"
    pids::W
end
PriorityGroup{U}(maxpriority::T) where {T,U} = PriorityGroup(maxpriority, 0, Vector{U}())

@inline function insert!(pg::PriorityGroup, pid)
    @unpack numpids, pids = pg

    if numpids == length(pids)
        push!(pids, pid)
    else
        @inbounds pids[numpids+1] = pid
    end

    # return the position of insertion
    pg.numpids += 1
end

@inline function remove!(pg::PriorityGroup, pididx)
    @unpack numpids, pids = pg

    @inbounds lastpid = pids[numpids]

    # simply swap the last id with the one to remove
    @inbounds pids[pididx] = lastpid
    pg.numpids -= 1

    # return the pid that was swapped to pididx
    lastpid
end

function Base.show(io::IO, pg::PriorityGroup)
    println("  ", summary(pg))
    println("  maxpriority = ", pg.maxpriority)
    println("  numpids = ", pg.numpids)
    println("  pids = ", pg.pids[1:pg.numpids])
end

"""
Table to store the groups.
"""
mutable struct PriorityTable{F,S,T,U<:Function}
    "non-zero values below this are binned together, static"
    minpriority::F

    "values above this cause new groups to be added"
    maxpriority::F

    "bins storing priority ids within a given range"
    groups::Vector{PriorityGroup{F,Vector{S}}}

    "stores the sums of the priorities within each group"
    gsums::Vector{F}

    "stores the sum of the group sums"
    gsum::F

    "maps priority id to group and idx within the group"
    pidtogroup::Vector{Tuple{T,T}}

    "mapping from priority value to group id that stores it"
    priortogid::U
end

"""
Setup table from a vector of priorities. The id
of a priority is its position within this vector.
"""
function PriorityTable(priortogid::Function, priorities::AbstractVector, minpriority, maxpriority)

    numgroups  = priortogid(maxpriority)
    numgroups -= one(typeof(numgroups))
    pidtype    = typeof(numgroups)
    ptype      = eltype(priorities)
    groups     = Vector{PriorityGroup{ptype,Vector{pidtype}}}()
    pidtogroup = Vector{Tuple{Int,Int}}(undef, length(priorities))
    gsum       = zero(ptype)
    gsums      = zeros(ptype, numgroups)

    # create the groups, {0}, (0,minpriority), [minpriority,2*minpriority)...
    push!(groups, PriorityGroup{pidtype}(zero(ptype)))
    gmaxprior = minpriority
    for i = 2:numgroups
        push!(groups, PriorityGroup{pidtype}(gmaxprior))
        gmaxprior *= 2
    end
    #@assert 2*maxpriority <= gmaxprior "gmaxprior = $gmaxprior, maxpriority=$maxpriority"
    gmaxprior = groups[end].maxpriority

    pt = PriorityTable(minpriority, gmaxprior, groups, gsums, gsum, pidtogroup, priortogid)

    # insert priority ids into the groups
    for (pid,priority) in enumerate(priorities)
        insert!(pt, pid, priority)
    end

    pt
end

########################## ACCESSORS ##########################
@inline function numgroups(pt::PriorityTable)
    length(pt.groups)
end

@inline function numpriorities(pt::PriorityTable)
    length(pt.pidtogroup)
end

@inline function groupsum(pt::PriorityTable)
    pt.gsum
end

"""
Adds extra groups to the table to accomodate a new maxpriority.
"""
function padtable!(pt::PriorityTable, pid, priority)
    @unpack maxpriority, groups, gsums = pt
    pidtype = typeof(pid)

    while priority >= maxpriority
        maxpriority *= 2
        push!(groups, PriorityGroup{pidtype}(maxpriority))
        push!(gsums, zero(eltype(gsums)))
    end
    pt.maxpriority = maxpriority
    nothing
end

# assumes pid is at most 1 greater than last priority (id) currently in table
# i.e. pid = length(pidtogroup) + 1
function insert!(pt::PriorityTable, pid, priority)
    @unpack maxpriority, groups, gsums, pidtogroup, priortogid = pt
    pidtype = typeof(pid)

    # find group for new priority
    gid = priortogid(priority)

    # add new (empty) groups if priority is too big
    if priority >= maxpriority
        padtable!(pt, pid, priority)
        #@assert (gid == length(groups))
    end

    # update table and priority's group
    pt.gsum += priority
    @inbounds gsums[gid] += priority
    @inbounds pididx = insert!(groups[gid], pid)

    if pid <= length(pidtogroup)
        @inbounds pidtogroup[pid] = (gid,pididx)
    else
        #@assert pid == length(pidtogroup) + 1
        push!(pidtogroup, (gid,pididx))
    end

    nothing
end

function update!(pt::PriorityTable, pid, oldpriority, newpriority)
    @unpack maxpriority, groups, gsums, pidtogroup, priortogid = pt

    oldgid = priortogid(oldpriority)
    newgid = priortogid(newpriority)

    # expand the table if necessary
    if newpriority >= maxpriority
        padtable!(pt, pid, newpriority)
    end

    # update the global priority sum
    pdiff    = newpriority - oldpriority
    pt.gsum += pdiff

    if oldgid == newgid
        # update the group priority too
        @inbounds gsums[newgid] += pdiff
    else
        @inbounds begin
            # location in oldgid of pid to move
            pididx = pidtogroup[pid][2]

            # remove pid from old group
            movedpid = remove!(groups[oldgid], pididx)

            # update position in group of the pid that swapped places with the removed one
            pidtogroup[movedpid] = (oldgid,pididx)

            # insert the updated pid back and store it's new position in the group
            pididx = insert!(groups[newgid], pid)
            pidtogroup[pid] = (newgid,pididx)

            # update sums, special case if group empty to avoid FP error in running sums
            grpsz = groups[oldgid].numpids
            gsums[oldgid]  = (grpsz == zero(grpsz)) ? zero(oldpriority) : gsums[oldgid] - oldpriority
            gsums[newgid] += newpriority
        end
    end
    nothing
end

function Base.show(io::IO, pt::PriorityTable)
    println(summary(pt))
    println("(minpriority,maxpriority) = ", (pt.minpriority, pt.maxpriority))
    println("sum of priorities = ", pt.gsum)
    println("num of groups = ", length(pt.groups))
    println("pidtogroup = ", pt.pidtogroup)
    for (i,group) in enumerate(pt.groups)
        if group.numpids > 0
            println("group = ",i,", group sum = ", pt.gsums[i])
            Base.show(io,group)
        end
    end
end

#######################################################
# routines for DirectCR
#######################################################

# map priority (i.e. jump rate) to integer
# add two as 0. -> 1 and  priority < minpriority ==> pid -> 2
@inline function priortogid(priority, mingid)
    (priority <= eps(typeof(priority))) && return 1
    gid = exponent(priority) + 1
    (gid <= mingid) && return 2
    return gid - mingid + 2
end

@inline function sample(pg::PriorityGroup, priorities, rng=Random.GLOBAL_RNG)
    @unpack maxpriority, numpids, pids = pg

    pididx = 0
    pid    = zero(eltype(pids))

    @inbounds while true

        # pick a random element
        r      = rand(rng) * numpids
        pididx = trunc(Int, r)
        pid    = pids[pididx+1]

        # acceptance test
        ( (r - pididx)*maxpriority < priorities[pid] ) && break
    end

    pid
end

function sample(pt::PriorityTable, priorities, rng=Random.GLOBAL_RNG)
    @unpack groups, gsum, gsums = pt

    # return id zero if total priority is zero
    (gsum < eps(typeof(gsum))) && return zero(eltype(groups[1].pids))

    # sample a group, search from end (largest priorities)
    # NOTE, THIS ASSUMES THE FIRST PRIORITY IS ZERO!!!
    # gid = length(gsums)
    # @inbounds r = rand(rng) * gsum - gsums[gid]
    # while r > zero(r)
    #     gid -= one(gid)
    #     iszero(gid) && return gid   # if no result found return zero
    #     @inbounds r -= gsums[gid]
    # end
    gid = 0
    r   = rand(rng) * gsum
    @inbounds for i = length(gsums):-1:1
        r -= gsums[i]
        if r <= zero(r)
            gid = i
            break
        end
    end
    iszero(gid) && return gid

    # sample element within the group
    @inbounds sample(groups[gid], priorities, rng)
end


text:
@testset "SVC" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x <= 5, Int)
    @variable(m, -5 <= y <= 5, Int)
    @constraint(m, x <= y)
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]
    @test constraint isa CS.SingleVariableConstraint

    # doesn't check the length
    @test !CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [3, 2])
    @test CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [2, 2])
    @test CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 2])

    constr_indices = constraint.indices
    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[2],
        -5,
    )
    @test CS.fix!(com, com.search_space[constr_indices[2]], 3)
    @test !CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[1],
        4,
    )
    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[1],
        3,
    )

    # need to create a backtrack_vec to reverse pruning
    dummy_backtrack_obj = CS.BacktrackObj(com)
    push!(com.backtrack_vec, dummy_backtrack_obj)
    # reverse previous fix
    CS.reverse_pruning!(com, 1)
    com.c_backtrack_idx = 1

    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[1],
        4,
    )

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x <= 5, Int)
    @variable(m, -5 <= y <= 5, Int)
    @constraint(m, x <= y)
    optimize!(m)
    com = CS.get_inner_model(m)
    constraint = com.constraints[1]

    # feasible and no changes
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    for ind in constr_indices
        @test sort(CS.values(com.search_space[ind])) == -5:5
    end
    @test CS.fix!(com, com.search_space[constr_indices[2]], 4)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    @test sort(CS.values(com.search_space[1])) == -5:4

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x <= 5, Int)
    @variable(m, -5 <= y <= 5, Int)
    @constraint(m, x <= y)
    optimize!(m)
    com = CS.get_inner_model(m)
    constraint = com.constraints[1]

    # Should be synced to the other variables
    @test CS.remove_above!(com, com.search_space[constr_indices[2]], 1)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    for ind in constr_indices
        @test sort(CS.values(com.search_space[ind])) == -5:1
    end
end

@testset "svc is_constraint_violated test" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, 0 <= x <= 5, Int)
    @variable(m, 0 <= y <= 3, Int)
    @constraint(m, x <= y)
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]

    variables = com.search_space
    @test CS.fix!(com, variables[constraint.indices[1]], 5; check_feasibility = false)
    @test CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, 0 <= x <= 5, Int)
    @variable(m, 0 <= y <= 3, Int)
    @constraint(m, x <= y)
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]

    variables = com.search_space
    @test CS.fix!(com, variables[constraint.indices[1]], 3; check_feasibility = false)
    @test !CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)
end


text:
using FrankWolfe
using LinearAlgebra
using Test

# necessary for away step
function Base.convert(::Type{GenericArray{Float64, 1}}, v::FrankWolfe.ScaledHotVector{Float64})
    return GenericArray(collect(v))
end

@testset "GenericArray is maintained" begin
    n = Int(1e3)
    k = n

    f(x::GenericArray) = dot(x, x)
    function grad!(storage, x)
        @. storage = 2 * x
    end

    lmo = FrankWolfe.ProbabilitySimplexOracle(1)

    x0 = GenericArray(collect(FrankWolfe.compute_extreme_point(lmo, zeros(n))))

    x, v, primal, dual_gap0, trajectory = FrankWolfe.frank_wolfe(
        f,
        grad!,
        lmo,
        x0,
        max_iteration=k,
        line_search=FrankWolfe.Agnostic(),
        print_iter=k / 10,
        verbose=true,
        emphasis=FrankWolfe.memory,
    )

    @test f(x) < f(x0)
    @test x isa GenericArray

    x, v, primal, dual_gap0, trajectory = FrankWolfe.lazified_conditional_gradient(
        f,
        grad!,
        lmo,
        x0,
        max_iteration=k,
        line_search=FrankWolfe.Agnostic(),
        print_iter=k / 10,
        verbose=true,
        emphasis=FrankWolfe.memory,
        VType=FrankWolfe.ScaledHotVector{Float64},
    )

    @test f(x) < f(x0)
    @test x isa GenericArray

    @test_broken x, v, primal, dual_gap0, trajectory = FrankWolfe.away_frank_wolfe(
        f,
        grad!,
        lmo,
        x0,
        max_iteration=k,
        line_search=FrankWolfe.Agnostic(),
        print_iter=k / 10,
        verbose=true,
        emphasis=FrankWolfe.memory,
    )
end


text:
@testset "_merge - AbstractDict" begin
    @testset "Simple" begin
        expected = Dict("a" => 1, "b" => 2)

        a = Dict("a"=>1)
        b = Dict("b"=>2)

        @test AWS._merge(a, b) == expected

    end

    @testset "Complex Different Nested Dict Keys" begin
        expected = Dict("common" => Dict("b" => 2))

        a = Dict("common" => Dict("a" => 1))
        b = Dict("common" => Dict("b" => 2))

        @test AWS._merge(a, b) == expected
    end

    @testset "Complex Same Nested Dict Keys" begin
        expected = Dict("common" => Dict("a" => 2))

        a1 = Dict("common" => Dict("a" => 1))
        a2 = Dict("common" => Dict("a" => 2))

        @test AWS._merge(a1, a2) == expected
    end
end

@testset "mergewith(_merge) - AbstractDict" begin
    @testset "Simple" begin
        expected = Dict("a" => 1, "b" => 2)

        a = Dict("a"=>1)
        b = Dict("b"=>2)

        @test mergewith(AWS._merge, a, b) == expected
    end

    @testset "Complex Differed Nested Dict Keys" begin
        expected = Dict("common" => Dict("a" => 1, "b" => 2))

        a = Dict("common" => Dict("a" => 1))
        b = Dict("common" => Dict("b" => 2))

        @test mergewith(AWS._merge, a, b) == expected
    end

    @testset "Complex Same Nested Dict Keys" begin
        expected = Dict("common" => Dict("a" => 2))

        a1 = Dict("common" => Dict("a" => 1))
        a2 = Dict("common" => Dict("a" => 2))

        @test mergewith(AWS._merge, a1, a2) == expected
    end
end

@testset "_merge - AbstractString" begin
    expected = "b"

    @test AWS._merge("a", "b") == expected
end


text:
module IterLSQ

include("diff_variable.jl")
include("fit.jl")
include("plot_fit.jl")

end # module


text:
type FullModel <: GP
	d::Int64
	k::Kernel
	lik::Likelihood
	X::Array{Float64, 2}
	y::Array{Float64, 1}
	y_raw::Array{Float64, 1}
	K::Array{Float64, 2} #Covariance matrix
	Cfull::Array{Float64, 2}
    C::Union(Array{Float64, 2}, Base.LinAlg.Cholesky{Float64}) #Cholesky decomp of C
    CinvY::Array{Float64, 1} # C^-1  * Y
    mu::Float64
    std::Float64
end

function FullModel(k::Kernel, lik::Likelihood)
	FullModel(0, k, lik, Array(Float64, (0, 0)), Array(Float64, (0,)),
		Array(Float64, (0,)),
		Array(Float64, (0,0)), Array(Float64, (0, 0)), Array(Float64, (0,0)),
		Array(Float64, (0,)), 0.0, 1.0)
end

function observe!{T <: FloatingPoint}(m::FullModel, X::Array{T}, y::Array{T})
	if size(m.X) == (0,0)
		m.d = size(X, 2)
		m.X = X
		m.y_raw = y
	else
		@assert m.d == size(X, 2) "Ensure X has the same dimension as the model!"
		m.X = vcat(m.X, X)
		m.y_raw = vcat(m.y_raw, y)
	end
	m.mu = mean(m.y_raw)
	m.std = std(m.y_raw)
	m.y = (m.y_raw .- m.mu) / (m.std + eps())
	refresh_cache!(m)
end

function predict(m::FullModel, Xp)
	Kpd = kernel(m.k, Xp, m.X)
	d = size(Xp, 2)
	z = zeros(1, d)
	c = kernel(m.k, z, z)[1]
	mn = (Kpd * m.CinvY)
	variance = (c .- sum(Kpd .* (m.C \ Kpd')', 2))
	mn * (m.std + eps()) .+ m.mu, variance * (m.std + eps())^2
end

function refresh_cache!(m::FullModel)
	m.K = kernel(m.k, m.X)
	sigma2 = exp(2m.lik.log_sigma)
	m.Cfull = m.K + eye(m.K) * sigma2
	try
		m.C = cholfact(m.Cfull)
	catch e
		m.C = m.Cfull
	end
	m.CinvY = vec(m.C \ m.y)
end


text:
abstract type CompositeFieldVector{T,L,N} <: AbstractVector{T} end

" determines is a field is included in the composite vector "
composite(typ, field) = true

(::Type{F})(typ, args...) where F <: CompositeFieldVector = begin
    lookup = buildlookup(typ, filter(F, args), 1)
    tlookup = tuplelookup(lookup)
    F{typ, tlookup, length(lookup), typeof.(args)...}(args...)
end

filter(F, args) = begin
    include = composite.(F, fieldnames(F))
    x = [args...]
    x[.!include] .= nothing
    x
end

buildlookup(t, fs::Array, i, l...) = begin
    if length(fs) > 1
        (buildlookup(t, fs[1], i, l...)..., buildlookup(t, fs[2:end], i + 1, l...)...)
    elseif length(fs) == 1
        (buildlookup(t, fs[1], i, l...)...,)
    else
        tuple()
    end
end
buildlookup(t, fs::Tuple, i, l...) = buildlookup(t, [fs...], 1, l..., i)
buildlookup(t, fs::Tuple{}, i, l...) = tuple()
buildlookup(t, f, i, l...) = begin
    typeof(f) <: Type(t) && return ((l..., i),)
    fns = fieldnames(f)
    if length(fns) > 0
        fns[composite.(typeof(f), fns)]
        buildlookup(t, Any[getfield.(f, fns)...], 1, l..., i)
    else
        tuple()
    end
end
@require Unitful begin
    buildlookup(t, f::T, i, l...) where T <: Unitful.Quantity =
        typeof(f.val) <: Type(t) ? ((l..., i),) : nothing
end

function tuplelookup(lookup)
    typelookup = []
    for l in lookup
        push!(typelookup, Tuple{l...})
    end
    Tuple{typelookup...}
end

body(x::UnionAll) = body(x.body)
body(x::DataType) = x


text:
using AbstractPlotting.MakieLayout

function app(
    points,
    F::Type=Alpha;
    palette=DEFAULT_PALETTE,
    resolution=(1600, 900),
    backgroundcolor=:white,
    time=Observable(0.0),
    slider_values=1000,
    dim_max=1,
    ripserer_kwargs...,
)
    @info "computing diagrams"
    diagram = ripserer(F, points; verbose=true, ripserer_kwargs..., dim_max=dim_max)
    filtration = F(points)
    return app(
        points,
        filtration,
        diagram;
        palette=palette,
        resolution=resolution,
        backgroundcolor=backgroundcolor,
        time=time,
        slider_values=slider_values,
        dim_max=dim_max,
    )
end

function app(
    points,
    filtration,
    diagram;
    palette=DEFAULT_PALETTE,
    resolution=(1600, 900),
    backgroundcolor=:white,
    time=Observable(0.0),
    slider_values=1000,
    dim_max=1,
)
    outer_padding = 20
    scene, layout = layoutscene(
        outer_padding; resolution=resolution, backgroundcolor=backgroundcolor
    )

    filtration_axis = layout[1:2, 1:2] = LScene(scene; title="Filtration")
    diagram_axis = layout[1, 3] = LAxis(scene; title="Diagram")
    bcd_axis = layout[2, 3] = LAxis(scene; title="Barcode")

    min_time = minimum(Ripserer.births(filtration))
    max_time = threshold(filtration)

    time_slider, t = time_slider!(scene, min_time, max_time, slider_values, time)
    toggles, triangles, edges, criticals = toggles!(scene, dim_max)
    triangles = @lift $triangles && $edges
    layout[3, 1:2] = time_slider
    layout[3, 3] = toggles

    markersize = length(points[1]) == 2 ? 5 : 10

    plot!(
        filtration_axis, filtration, points; time=t, palette, triangles, edges, markersize
    )
    plot_critical_chains!(
        filtration_axis, diagram, points, criticals; time=t, palette, markersize
    )
    if cameracontrols(filtration_axis.scene) isa Camera3D
        cameracontrols(filtration_axis.scene).rotationspeed[] = 0.01f0
    end

    plot_diagram!(diagram_axis, diagram; time=t, palette=palette)
    plot_barcode!(bcd_axis, diagram; time=t, palette=palette)

    linkxaxes!(diagram_axis, bcd_axis)
    tightlimits!(diagram_axis)
    tightlimits!(bcd_axis)

    return scene
end

function time_slider!(scene, min_time, max_time, slider_values, time)
    time_slider = LSlider(scene; range=range(min_time, max_time; length=slider_values))
    time_str = @lift string("t=", round($(time_slider.value); digits=3))
    time_label = LText(scene, time_str; tellwidth=false)
    set_close_to!(time_slider, to_value(time))
    on(time) do t
        set_close_to!(time_slider, t)
    end
    return (
        grid!(reshape([time_label, time_slider], (2, 1)); tellwidth=false),
        time_slider.value,
    )
end

function toggles!(scene, dim_max)
    triangle_toggle = LToggle(scene; active=true, tellwidth=false)
    triangle_label = LText(scene, "triangles"; tellwidth=false)
    edge_toggle = LToggle(scene; active=true, tellwidth=false)
    edge_label = LText(scene, "edges"; tellwidth=false)
    critical_toggles = [LToggle(scene; active=false, tellwidth=false) for d in 0:dim_max]
    critical_labels = [
        LText(scene, "$d-critical simplices"; tellwidth=false) for d in 0:dim_max
    ]

    layout = [
        triangle_toggle triangle_label
        edge_toggle edge_label
        critical_toggles critical_labels
    ]
    return (
        grid!(layout),
        triangle_toggle.active,
        edge_toggle.active,
        [t.active for t in critical_toggles],
    )
end

function plot_critical_chains!(scene, diagram, points, criticals; time, kwargs...)
    death_chains = [
        FilteredChain(filter(!isnothing, death_simplex.(d)), points) for d in diagram
    ]
    birth_chains = [FilteredChain(birth_simplex.(d), points) for d in diagram]
    crit_ts = [@lift $(criticals[i]) ? $time : -Inf for i in eachindex(diagram)]
    for (birth_chain, death_chain, time, i) in
        zip(birth_chains, death_chains, crit_ts, eachindex(crit_ts))
        for chain in (birth_chain, death_chain)
            if length(chain.triangles) == 0
                chainplot!(scene, chain; edgecolor=2 + i, linewidth=3, time, kwargs...)
            else
                chainplot!(scene, chain; trianglecolor=2 + i, time, kwargs...)
            end
        end
    end
end


text:
#module types

#export Projectile

mutable struct RectErrorB
    σ_x::Float64
    σ_y::Float64
    μ_x::Float64
    μ_y::Float64
end

mutable struct Hull
    Φ::Union{Float64,Nothing}
    ψ::Union{Float64,Nothing}
    χ::Union{Float64,Nothing}
end
"""
    The turret is defined by its azimuth angle
"""
mutable struct Turret
     ξ::Union{Float64,Nothing}#azimuth
end

mutable struct Canon
     θ::Union{Float64,Nothing} #elevation
     lw::Float64
     u₀::Float64
end

 mutable struct Tourelle
     θ::Union{Float64,Nothing} #elevation
     ξ::Union{Float64,Nothing}#azimuth
     Φ::Union{Float64,Nothing}
     ψ::Union{Float64,Nothing}
     χ::Union{Float64,Nothing}
     rWY::Union{Float64,Nothing}
 end

 mutable struct Sight
     θ::Float64
     ξ::Float64
 end

mutable struct Tank
    hull::Hull
    turret::Turret
    canon::Canon
    sight::Sight
    rWY::Float64
end

 #mutable struct Target
#     ρ::Float64
# end


text:
#TODO Rewrite battery problem
"""
    run_battery_opt(data::ClustData)

operational battery storage optimization problem
runs every day seperately and adds results in the end
"""
function run_battery_opt(data::ClustData)
  prnt=false
  num_periods = data.K # number of periods, 1day, one week, etc.
  num_hours = data.T # hours per period (24 per day, 48 per 2days)
  el_price = data.data["el_price-$(data.region)"]
  weight = data.weights
  # time steps
  del_t = 1; # hour

  # example battery Southern California Edison
  P_battery = 100; # MW
  E_battery = 400; # MWh
  eff_Storage_in = 0.95;
  eff_Storage_out = 0.95;
  #Stor_init = 0.5;

  # optimization
  # Sets
  # time
  t_max = num_hours;

  E_in_arr = zeros(num_hours,num_periods)
  E_out_arr = zeros(num_hours,num_periods)
  stor = zeros(num_hours +1,num_periods)

  obj = zeros(num_periods);
  m= Model(optimizer=ClpSolver() )

  # hourly energy output
  @variable(m, E_out[t=1:t_max] >= 0) # kWh
  # hourly energy input
  @variable(m, E_in[t=1:t_max] >=# optimization problems
   0) # kWh
  # storage level
  @variable(m, Stor_lev[t=1:t_max+1] >= 0) # kWh

  @variable(m,0 <= Stor_init <= 1) # this as a variable ensures

  # maximum battery power
  for t=1:t_max
    @constraint(m, E_out[t] <= P_battery*del_t)
    @constraint(m, E_in[t] <= P_battery*del_t)
  end

  # maximum storage level
  for t=1:t_max+1
    @constraint(m, Stor_lev[t] <= E_battery)
  end

  # battery energy balance
  for t=1:t_max
    @constraint(m,Stor_lev[t+1] == Stor_lev[t] + eff_Storage_in*del_t*E_in[t]-(1/eff_Storage_out)*del_t*E_out[t])
  end

  # initial storage level
  @constraint(m,Stor_lev[1] == Stor_init*E_battery)
  @constraint(m,Stor_lev[t_max+1] >= Stor_lev[1])
  s=:Optimal
  for i =1:num_periods
    #objective
    @objective(m, Max, sum((E_out[t] - E_in[t])*el_price[t,i] for t=1:t_max) )
    status = solve(m)
    if status != :Optimal
      s=:NotSolved
    end
    if weight ==1
      obj[i] = getobjectivevalue(m)
    else
      obj[i] = getobjectivevalue(m) * weight[i]
    end
    E_in_arr[:,i] = getvalue(E_in)
    E_out_arr[:,i] = getvalue(E_out)
    stor[:,i] = getvalue(Stor_lev)
  end
  vars= Dict()
  vars["E_out"] = OptVariable(E_out_arr,"operation")
  vars["E_in"] = OptVariable(E_in_arr,"operation")
  vars["Stor_level"] = OptVariable(stor,"operation")
  res = OptResult(s,sum(obj),vars,Dict())
  return res
end # run_battery_opt()

 ###

"""
    run_gas_opt(data::ClustData)

operational gas turbine optimization problem
runs every day seperately and adds results in the end
"""
function run_gas_opt(data::ClustData)

  prnt=false
  num_periods = data.K # number of periods, 1day, one week, etc.
  num_hours = data.T # hours per period (24 per day, 48 per 2days)
  el_price = data.data["el_price"]
  weight = data.weights
  # time steps
  del_t = 1; # hour

  # example gas turbine
  P_gt = 100; # MW
  eta_t = 0.6; # 40 % efficiency
  if data.region == "GER"
    gas_price = 24.65  # EUR/MWh    7.6$/GJ = 27.36 $/MWh=24.65EUR/MWh with 2015 conversion rate
  elseif data.region == "CA"
    gas_price  = 14.40   # $/MWh        4$/GJ = 14.4 $/MWh
  end

  # optimization
  # Sets
  # time,
  t_max = num_hours;

  E_out_arr = zeros(num_hours,num_periods)

  obj = zeros(num_periods);
  m= Model(optimizer=ClpSolver() )

  # hourly energy output
  @variable(m, 0 <= E_out[t=1:t_max] <= P_gt) # MWh

  s=:Optimal
  for i =1:num_periods
    #objective
    @objective(m, Max, sum(E_out[t]*el_price[t,i] - 1/eta_t*E_out[t]*gas_price for t=1:t_max) )
    status = solve(m)
    if status != :Optimal
      s=:NotSolved
    end

    if weight ==1
      obj[i] = getobjectivevalue(m)
    else
      obj[i] = getobjectivevalue(m) * weight[i]
    end
    E_out_arr[:,i] = getvalue(E_out)
  end

  op_vars= Dict()
  op_vars["E_out"] = E_out_arr
  res = OptResult(s,sum(obj),Dict(),op_vars,Dict())
  return res
end # run_gas_opt()


text:
export Event

abstract OutgoingEvent
abstract Event

macro slackevent(event_type::Symbol, event_name::AbstractString, event_block::Expr)
    quote
        immutable $event_type <: Event
            $(event_block.args...)
        end

        $(esc(:register_event))($event_name, $event_type)
    end
end

macro slackoutevent(event_type::Symbol, event_name::AbstractString, event_block::Expr)
    quote
        immutable $event_type <: OutgoingEvent
            $(event_block.args...)
        end

        $(esc(:outgoing_event_name))(::Type{$(esc(event_type))}) = $event_name
    end
end

function serialize(e::OutgoingEvent)
    s = Dict{UTF8String, Any}(
        utf8("type") => outgoing_event_name(typeof(e)))
    for f in fieldnames(e)
        s[utf8(string(f))] = getfield(e, f)
    end
    s
end


text:
"""
Transform an object into a format that is accepted by MeshCat.
"""
function render_obj(obj::Object)::MeshFileGeometry
    points = transform(obj.mesh.verts, obj.pose)
    lines = ["v " * join(p, " ") for p in eachcol(points)]
    for (a, b, c, d) in eachcol(obj.mesh.tets)
        for (x, y, z) in [(a, b, c), (b, c, d), (c, d, a), (d, a, b)]
            # Only include faces with potential 1, to get the outer surface.
            if obj.mesh.potentials[[x, y, z]] ≈ [0, 0, 0]
                push!(lines, "f $x $y $z")
            end
        end
    end
    MeshFileGeometry(join(lines, "\n"), "obj")
end

export render_obj

max_stars_repo_name: ekzhang/hydroelastics

max_issues_repo_name: ekzhang/hydroelastics

max_forks_repo_name: ekzhang/hydroelastics


text:
@testset "Testing utils.jl" begin

    bucket = [1,2,3,4,5,6]
    @test AbrkgaObop.create_bucket_order(bucket) == " 1 | 2 | 3 | 4 | 5 | 6 "

    bucket = [1,1,1,1,1,1]
    @test AbrkgaObop.create_bucket_order(bucket) == " 1 , 2 , 3 , 4 , 5 , 6 "

    bucket = [1,1,1,2,2,3]
    @test AbrkgaObop.create_bucket_order(bucket) == " 1 , 2 , 3 | 4 , 5 | 6 "

    bucket = [1,2,3,2,1,1]
    @test AbrkgaObop.create_bucket_order(bucket) == " 1 , 5 , 6 | 2 , 4 | 3 "

end

max_stars_repo_name: LuizHNLorena/AbrkgaObop

max_issues_repo_name: LuizHNLorena/AbrkgaObop

max_forks_repo_name: LuizHNLorena/AbrkgaObop


text:
function sweepline(V,EV)
	# event generation and ordering
	segments = presortedlines(V,EV)
	evpairs = [[(v1,v2,"start",k), (v2,v1,"end",k)] for (k,(v1,v2)) in enumerate(segments)]
	isless0 = (x,y) -> [x[1][1],x[1][2],x[2][2]] < [y[1][1],y[1][2],y[2][2]]
	events = sort(cat(evpairs),lt=isless0)
	# Initialize event queue ξ = all segment endpoints; Sort ξ by increasing x and y
	pqkeys = [(e[1],e[3],e[4]) for e in events]
	pqvalues = events
	ξ = PriorityQueue(zip(pqkeys,pqvalues))
@show ξ,0 println()
	# Proper ordering object for the `SortedMultiDic` used by sweep line SL
	##o = Lt((x,y) -> (x,y))
	# Initialize sweep line SL to be empty
	SL = SortedMultiDict{Any,Any}()##o)
@show SL, 0 println()
	# Initialized output intersection list Λ to be empty
	Λ = Array{Array{Float64,1},1}[]
@show Λ,0 println()

	while length(ξ) ≠ 0 # (ξ is nonempty)
		E = peek(ξ)[2]   # the next value of event (k => v) from ξ
@show E, 1 println()
		(v1, v2, nodetype, edgeid) = E
		# segE = E's segment
		e, segE = v1, (v1, v2, nodetype, edgeid) # key, value
		if E[3] == "start" # (E is a left endpoint)
			# Add segE to SL
			i = insert!(SL, e, segE)  # SortedMultiDict, key, value -> semitoken
@show SL, 1 println()
			global xe,ye = e
			if first(SL) !== last(SL) # more than one segment in SL
				# compute segA and/or segB
				segA, segB = selectsegmentneighbor(SL,i)
				if segA ≠ []
				# segA = the segment above segE in SL
					a = Int(segA[3][1]); segA = segA[1:2] # edgeid of segA
					I = intersection(segE,segA)
					# (if Intersect( segE with segA) exists)
					if typeof(I) ≠ Nothing
						# Insert I into ξ
						#(keyE,newsegE),(keyA,newsegA) = swapsegments(SL,I,segE,segA)
						key = (I,"int",edgeid); val = (I,I,"int",a)
						enqueue!(ξ, key,val)
						#enqueue!(ξ, keyA,newsegA)
@show ξ,6 println()
					end
				end
				if segB ≠ []
				# segB = the segment below segE in SL
					b = Int(segB[3][1]); segB = segB[1:2] # edgeid of segB
					I = intersection(segE,segB)
					# (if Intersect( segE with segB) exists)
					if typeof(I) ≠ Nothing
						# Insert I into ξ
						#(keyE,newsegE),(keyB,newsegB) = swapsegments(SL,I,segE,segB)
						key = (I,"int",edgeid); val = (I,I,"int",b)
						enqueue!(ξ, key,val)
						#enqueue!(ξ, keyB,newsegB)
@show ξ,5 println()
					end
				end
			end
		elseif E[3] == "end" # (E is a right endpoint)
			# segE = E's segment
			# e = key(segE); st1 = corresponding semitoken in SL
			st1,st2 = searchequalrange(SL::SortedMultiDict, e)
			# compute segA and/or segB
			segA, segB = selectsegmentneighbor(SL,st1)
			# segA = the segment above segE in SL
			if segA ≠ []
				a = Int(segA[3][1]) # edgeid of segA
				segA = segA[1:2]
			end
			# segB = the segment below segE in SL
			if segB ≠ []
				b = Int(segB[3][1])
				segB = segB[1:2] # edgeid of segB
			end
			# Remove segE from SL
@show st1,st2 println()
			##delete!((SL,st1))
			# (I = Intersect( segA with segB) exists)
			if segA ≠ [] && segB ≠ []
				I = intersection(segA, segB)
				if typeof(I) ≠ Nothing
					# (I is not in ξ already)
					# no problem anyway (in case I is overwritten)
					# Insert I into ξ
					#(keyA,newsegA),(keyB,newsegB) = swapsegments(SL,I,segA,segB)
					key = (I,"int",a); val = (I,I,"int",b)
					enqueue!(ξ, key,val)
					#enqueue!(ξ, keyB,newsegB)
@show ξ,4 println()
				end
			end
		else # E is an intersection event
		@assert E[3] == "int"
			# Add E to the output list Λ
			push!(Λ, [E[1],E[2]])
@show Λ,1 println()
			# segE1 above segE2 be E's intersecting segments in SL
			i = insert!(SL, e, segE) ### ????
			segE1, segE2 = selectsegmentneighbor(SL,i)
			if segE1 ≠ []
				h = Int(segE1[3][1]) # index of segment segE1
				v1 = V[:,EV[h][1]] # key of segE1 in SL
				hst1,hst2 = searchequalrange(SL::SortedMultiDict, v1)
				# extreme semitokens of segE1 in SL
				# segB = the segment below segE1 in SL
				segE1, segB = selectsegmentneighbor(SL,hst1)
				if segE1 ≠ [] && segB ≠ []
					I = intersection(segE1,segB)
					# (I = Intersect(segE1 with segB) exists)
					if typeof(I) ≠ Nothing
						# (if I is not in ξ already) no problem in case
						# Insert I into ξ
						(keyE1,newsegE1),(keyB,newsegB) = swapsegments(SL,I,segE1,segB)
						key = (I,"int",edgeid); val = (I,I,"int",b) # edgeid ??
						enqueue!(ξ, key,val)
						#enqueue!(ξ, keyB,newsegB)
@show ξ,5 println()
					end
				end
			end
			if segE2 ≠ []
				k = Int(segE2[3][1])  # index of segment segE2
				v1 = V[:,EV[k][1]]  # key of segE2 in SL
				kst1,kst2 = searchequalrange(SL::SortedMultiDict, v1)
				# extreme semitokens of segE2 in SL
				# Swap their positions so that segE2 is now above segE1
				# ??? HOW TO DO IT ???
				# segA = the segment above segE2 in SL
@show SL,kst1 println()
				segA, segE2 = selectsegmentneighbor(SL,kst1)
				if segA ≠ [] && segE2 ≠ []
					I = intersection(segA, segE2)
					# (I = Intersect(segE2 with segA) exists)
					if typeof(I) ≠ Nothing
						# (if I is not in ξ already) no problem in case
						# Insert I into ξ
						(keyA,newsegA),(keyE2,newsegE2) = swapsegments(SL,I,segA,segE2)
						##key = (I,"int",edgeid); val = (I,I,"int",a) # edgeid ??
						enqueue!(ξ, keyA,newsegA)
						enqueue!(ξ, keyE2,newsegE2)
@show ξ,4 println()
					end
				end
			end
			delete!((SL, i)) ### ????
		end
		# remove E from ξ
@show ξ,1 println()
		dequeue!(ξ)
@show ξ,2 println()
	end
	return Λ
end


text:
using Oceananigans.Grids: Center, Face

const RG   = RectilinearGrid
const RGX  = XRegRectilinearGrid
const RGY  = YRegRectilinearGrid
const RGZ  = ZRegRectilinearGrid

const CCSG = ConformalCubedSphereFaceGrid

const LLG  = LatitudeLongitudeGrid
const LLGX = XRegLatLonGrid
const LLGY = YRegLatLonGrid
const LLGZ = ZRegLatLonGrid

# On the fly calculations of metrics
const LLGF  = LatitudeLongitudeGrid{<:Any, <:Any, <:Any, <:Any, <:Nothing}
const LLGFX = LatitudeLongitudeGrid{<:Any, <:Any, <:Any, <:Any, <:Nothing, <:Any, <:Number}
const LLGFY = LatitudeLongitudeGrid{<:Any, <:Any, <:Any, <:Any, <:Nothing, <:Any, <:Any, <:Number}

@inline hack_cosd(φ) = cos(π * φ / 180)
@inline hack_sind(φ) = sin(π * φ / 180)

"""
Notes:
This file defines grid lengths, areas, and volumes for staggered structured grids.
Each "reference cell" is associated with an index i, j, k.
The "location" of each reference cell is roughly the geometric centroid of the reference cell.
On the staggered grid, there are 7 cells additional to the "reference cell"
that are staggered with respect to the reference cell in x, y, and/or z.

The staggering is denoted by the locations "Center" and "Face":
  - "Center" is shared with the reference cell;
  - "Face" lies between reference cell centers, roughly at the interface between
    reference cells.

The three-dimensional location of an object is defined by a 3-tuple of locations, and
denoted by a triplet of superscripts. For example, an object `φ` whose cell is located at
(Center, Center, Face) is denoted `φᶜᶜᶠ`. `ᶜᶜᶠ` is Centered in `x`, `Centered` in `y`, and on
reference cell interfaces in `z` (this is where the vertical velocity is located, for example).

The operators in this file fall into three categories:

1. Operators needed for an algorithm valid on rectilinear grids with
   at most a stretched vertical dimension and regular horizontal dimensions.
2. Operators needed for an algorithm on a grid that is curvilinear in the horizontal
   at rectilinear (possibly stretched) in the vertical.

"""

#####
#####
##### Spacings!!
#####
#####

@inline Δxᶠᵃᵃ(i, j, k, grid) =  nothing
@inline Δxᶜᵃᵃ(i, j, k, grid) =  nothing
@inline Δyᵃᶠᵃ(i, j, k, grid) =  nothing
@inline Δyᵃᶜᵃ(i, j, k, grid) =  nothing

const ZRG = Union{LLGZ, RGZ}

@inline Δzᵃᵃᶠ(i, j, k, grid) = @inbounds grid.Δzᵃᵃᶠ[k]
@inline Δzᵃᵃᶜ(i, j, k, grid) = @inbounds grid.Δzᵃᵃᶜ[k]

@inline Δzᵃᵃᶠ(i, j, k, grid::ZRG) = @inbounds grid.Δzᵃᵃᶠ
@inline Δzᵃᵃᶜ(i, j, k, grid::ZRG) = @inbounds grid.Δzᵃᵃᶜ

@inline Δzᵃᵃᶜ(i, j, k, grid::CCSG) = grid.Δz
@inline Δzᵃᵃᶠ(i, j, k, grid::CCSG) = grid.Δz

# Convenience Functions for all grids
for LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ)

    x_spacing_1D = Symbol(:Δx, LX, :ᵃ, :ᵃ)
    x_spacing_2D = Symbol(:Δx, LX, LY, :ᵃ)

    y_spacing_1D = Symbol(:Δy, :ᵃ, LY, :ᵃ)
    y_spacing_2D = Symbol(:Δy, LX, LY, :ᵃ)

    @eval begin
        @inline $x_spacing_2D(i, j, k, grid) = $x_spacing_1D(i, j, k, grid)
        @inline $y_spacing_2D(i, j, k, grid) = $y_spacing_1D(i, j, k, grid)
    end

    for LZ in (:ᶜ, :ᶠ)
        x_spacing_3D = Symbol(:Δx, LX, LY, LZ)
        y_spacing_3D = Symbol(:Δy, LX, LY, LZ)

        z_spacing_1D = Symbol(:Δz, :ᵃ, :ᵃ, LZ)
        z_spacing_3D = Symbol(:Δz, LX, LY, LZ)

        @eval begin
            @inline $x_spacing_3D(i, j, k, grid) = $x_spacing_2D(i, j, k, grid)
            @inline $y_spacing_3D(i, j, k, grid) = $y_spacing_2D(i, j, k, grid)
            @inline $z_spacing_3D(i, j, k, grid) = $z_spacing_1D(i, j, k, grid)
        end
    end
end

#####
##### Rectilinear Grids (Flat grids already have Δ = 1)
#####

@inline Δxᶠᵃᵃ(i, j, k, grid::RG)  =  @inbounds grid.Δxᶠᵃᵃ[i]
@inline Δxᶜᵃᵃ(i, j, k, grid::RG)  =  @inbounds grid.Δxᶜᵃᵃ[i]
@inline Δyᵃᶠᵃ(i, j, k, grid::RG)  =  @inbounds grid.Δyᵃᶠᵃ[j]
@inline Δyᵃᶜᵃ(i, j, k, grid::RG)  =  @inbounds grid.Δyᵃᶜᵃ[j]

@inline Δxᶠᵃᵃ(i, j, k, grid::RGX) =  @inbounds grid.Δxᶠᵃᵃ
@inline Δxᶜᵃᵃ(i, j, k, grid::RGX) =  @inbounds grid.Δxᶜᵃᵃ
@inline Δyᵃᶠᵃ(i, j, k, grid::RGY) =  @inbounds grid.Δyᵃᶠᵃ
@inline Δyᵃᶜᵃ(i, j, k, grid::RGY) =  @inbounds grid.Δyᵃᶜᵃ

#####
##### LatitudeLongitudeGrid
#####

## Pre computed metrics

@inline Δxᶜᶠᵃ(i, j, k, grid::LLG)  = @inbounds grid.Δxᶜᶠᵃ[i, j]
@inline Δxᶠᶜᵃ(i, j, k, grid::LLG)  = @inbounds grid.Δxᶠᶜᵃ[i, j]
@inline Δxᶠᶠᵃ(i, j, k, grid::LLG)  = @inbounds grid.Δxᶠᶠᵃ[i, j]
@inline Δxᶜᶜᵃ(i, j, k, grid::LLG)  = @inbounds grid.Δxᶜᶜᵃ[i, j]
@inline Δxᶠᶜᵃ(i, j, k, grid::LLGX) = @inbounds grid.Δxᶠᶜᵃ[j]
@inline Δxᶜᶠᵃ(i, j, k, grid::LLGX) = @inbounds grid.Δxᶜᶠᵃ[j]
@inline Δxᶠᶠᵃ(i, j, k, grid::LLGX) = @inbounds grid.Δxᶠᶠᵃ[j]
@inline Δxᶜᶜᵃ(i, j, k, grid::LLGX) = @inbounds grid.Δxᶜᶜᵃ[j]

@inline Δyᶜᶠᵃ(i, j, k, grid::LLG)  = @inbounds grid.Δyᶜᶠᵃ[j]
@inline Δyᶠᶜᵃ(i, j, k, grid::LLG)  = @inbounds grid.Δyᶠᶜᵃ[j]
@inline Δyᶜᶠᵃ(i, j, k, grid::LLGY) = @inbounds grid.Δyᶜᶠᵃ
@inline Δyᶠᶜᵃ(i, j, k, grid::LLGY) = @inbounds grid.Δyᶠᶜᵃ
@inline Δyᶜᶜᵃ(i, j, k, grid::LLG)  = Δyᶠᶜᵃ(i, j, k, grid)
@inline Δyᶠᶠᵃ(i, j, k, grid::LLG)  = Δyᶜᶠᵃ(i, j, k, grid)

## On the fly metrics

@inline Δxᶠᶜᵃ(i, j, k, grid::LLGF)  = grid.radius * deg2rad(grid.Δλᶠᵃᵃ[i]) * hack_cosd(grid.φᵃᶜᵃ[j])
@inline Δxᶜᶠᵃ(i, j, k, grid::LLGF)  = grid.radius * deg2rad(grid.Δλᶜᵃᵃ[i]) * hack_cosd(grid.φᵃᶠᵃ[j])
@inline Δxᶠᶠᵃ(i, j, k, grid::LLGF)  = grid.radius * deg2rad(grid.Δλᶠᵃᵃ[i]) * hack_cosd(grid.φᵃᶠᵃ[j])
@inline Δxᶜᶜᵃ(i, j, k, grid::LLGF)  = grid.radius * deg2rad(grid.Δλᶜᵃᵃ[i]) * hack_cosd(grid.φᵃᶜᵃ[j])
@inline Δxᶠᶜᵃ(i, j, k, grid::LLGFX) = grid.radius * deg2rad(grid.Δλᶠᵃᵃ)    * hack_cosd(grid.φᵃᶜᵃ[j])
@inline Δxᶜᶠᵃ(i, j, k, grid::LLGFX) = grid.radius * deg2rad(grid.Δλᶜᵃᵃ)    * hack_cosd(grid.φᵃᶠᵃ[j])
@inline Δxᶠᶠᵃ(i, j, k, grid::LLGFX) = grid.radius * deg2rad(grid.Δλᶠᵃᵃ)    * hack_cosd(grid.φᵃᶠᵃ[j])
@inline Δxᶜᶜᵃ(i, j, k, grid::LLGFX) = grid.radius * deg2rad(grid.Δλᶜᵃᵃ)    * hack_cosd(grid.φᵃᶜᵃ[j])

@inline Δyᶜᶠᵃ(i, j, k, grid::LLGF)  = grid.radius * deg2rad(grid.Δφᵃᶠᵃ[j])
@inline Δyᶠᶜᵃ(i, j, k, grid::LLGF)  = grid.radius * deg2rad(grid.Δφᵃᶜᵃ[j])
@inline Δyᶜᶠᵃ(i, j, k, grid::LLGFY) = grid.radius * deg2rad(grid.Δφᵃᶠᵃ)
@inline Δyᶠᶜᵃ(i, j, k, grid::LLGFY) = grid.radius * deg2rad(grid.Δφᵃᶜᵃ)

#####
#####  ConformalCubedSphereFaceGrid
#####

@inline Δxᶜᶜᵃ(i, j, k, grid::CCSG) = @inbounds grid.Δxᶜᶜᵃ[i, j]
@inline Δxᶠᶜᵃ(i, j, k, grid::CCSG) = @inbounds grid.Δxᶠᶜᵃ[i, j]
@inline Δxᶜᶠᵃ(i, j, k, grid::CCSG) = @inbounds grid.Δxᶜᶠᵃ[i, j]
@inline Δxᶠᶠᵃ(i, j, k, grid::CCSG) = @inbounds grid.Δxᶠᶠᵃ[i, j]

@inline Δyᶜᶜᵃ(i, j, k, grid::CCSG) = @inbounds grid.Δyᶜᶜᵃ[i, j]
@inline Δyᶠᶜᵃ(i, j, k, grid::CCSG) = @inbounds grid.Δyᶠᶜᵃ[i, j]
@inline Δyᶜᶠᵃ(i, j, k, grid::CCSG) = @inbounds grid.Δyᶜᶠᵃ[i, j]
@inline Δyᶠᶠᵃ(i, j, k, grid::CCSG) = @inbounds grid.Δyᶠᶠᵃ[i, j]

#####
#####
##### Areas!!
#####
#####

for LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ)

    x_spacing_2D = Symbol(:Δx, LX, LY, :ᵃ)
    y_spacing_2D = Symbol(:Δy, LX, LY, :ᵃ)
    z_area_2D    = Symbol(:Az, LX, LY, :ᵃ)

    @eval $z_area_2D(i, j, k, grid) = $x_spacing_2D(i, j, k, grid) * $y_spacing_2D(i, j, k, grid)

    for LZ in (:ᶜ, :ᶠ)
        x_spacing_3D = Symbol(:Δx, LX, LY, LZ)
        y_spacing_3D = Symbol(:Δy, LX, LY, LZ)
        z_spacing_3D = Symbol(:Δz, LX, LY, LZ)

        x_area_3D = Symbol(:Ax, LX, LY, LZ)
        y_area_3D = Symbol(:Ay, LX, LY, LZ)
        z_area_3D = Symbol(:Az, LX, LY, LZ)

        @eval begin
            @inline $x_area_3D(i, j, k, grid) = $y_spacing_3D(i, j, k, grid) * $z_spacing_3D(i, j, k, grid)
            @inline $y_area_3D(i, j, k, grid) = $x_spacing_3D(i, j, k, grid) * $z_spacing_3D(i, j, k, grid)
            @inline $z_area_3D(i, j, k, grid) = $z_area_2D(i, j, k, grid)
        end
    end
end

####
#### Special 2D z Areas for LatitudeLongitudeGrid and ConformalCubedSphereFaceGrid
####

@inline Azᶠᶜᵃ(i, j, k, grid::LLGF)  = grid.radius^2 * deg2rad(grid.Δλᶠᵃᵃ[i]) * (hack_sind(grid.φᵃᶠᵃ[j+1]) - hack_sind(grid.φᵃᶠᵃ[j]))
@inline Azᶜᶠᵃ(i, j, k, grid::LLGF)  = grid.radius^2 * deg2rad(grid.Δλᶜᵃᵃ[i]) * (hack_sind(grid.φᵃᶜᵃ[j])   - hack_sind(grid.φᵃᶜᵃ[j-1]))
@inline Azᶠᶠᵃ(i, j, k, grid::LLGF)  = grid.radius^2 * deg2rad(grid.Δλᶠᵃᵃ[i]) * (hack_sind(grid.φᵃᶜᵃ[j])   - hack_sind(grid.φᵃᶜᵃ[j-1]))
@inline Azᶜᶜᵃ(i, j, k, grid::LLGF)  = grid.radius^2 * deg2rad(grid.Δλᶜᵃᵃ[i]) * (hack_sind(grid.φᵃᶠᵃ[j+1]) - hack_sind(grid.φᵃᶠᵃ[j]))
@inline Azᶠᶜᵃ(i, j, k, grid::LLGFX) = grid.radius^2 * deg2rad(grid.Δλᶠᵃᵃ)    * (hack_sind(grid.φᵃᶠᵃ[j+1]) - hack_sind(grid.φᵃᶠᵃ[j]))
@inline Azᶜᶠᵃ(i, j, k, grid::LLGFX) = grid.radius^2 * deg2rad(grid.Δλᶜᵃᵃ)    * (hack_sind(grid.φᵃᶜᵃ[j])   - hack_sind(grid.φᵃᶜᵃ[j-1]))
@inline Azᶠᶠᵃ(i, j, k, grid::LLGFX) = grid.radius^2 * deg2rad(grid.Δλᶠᵃᵃ)    * (hack_sind(grid.φᵃᶜᵃ[j])   - hack_sind(grid.φᵃᶜᵃ[j-1]))
@inline Azᶜᶜᵃ(i, j, k, grid::LLGFX) = grid.radius^2 * deg2rad(grid.Δλᶜᵃᵃ)    * (hack_sind(grid.φᵃᶠᵃ[j+1]) - hack_sind(grid.φᵃᶠᵃ[j]))

for LX in (:ᶠ, :ᶜ), LY in (:ᶠ, :ᶜ)

    z_area_2D = Symbol(:Az, LX, LY, :ᵃ)

    @eval begin
        @inline $z_area_2D(i, j, k, grid::CCSG) = @inbounds grid.$z_area_2D[i, j]
        @inline $z_area_2D(i, j, k, grid::LLG)  = @inbounds grid.$z_area_2D[i, j]
        @inline $z_area_2D(i, j, k, grid::LLGX) = @inbounds grid.$z_area_2D[j]
    end
end

#####
#####
##### Volumes!!
#####
#####

@inline Vᶜᶜᶜ(i, j, k, grid) = Azᶜᶜᶜ(i, j, k, grid) * Δzᶜᶜᶜ(i, j, k, grid)
@inline Vᶠᶜᶜ(i, j, k, grid) = Azᶠᶜᶜ(i, j, k, grid) * Δzᶠᶜᶜ(i, j, k, grid)
@inline Vᶜᶠᶜ(i, j, k, grid) = Azᶜᶠᶜ(i, j, k, grid) * Δzᶜᶠᶜ(i, j, k, grid)
@inline Vᶜᶜᶠ(i, j, k, grid) = Azᶜᶜᶠ(i, j, k, grid) * Δzᶜᶜᶠ(i, j, k, grid)
@inline Vᶠᶠᶜ(i, j, k, grid) = Azᶠᶠᶜ(i, j, k, grid) * Δzᶠᶠᶜ(i, j, k, grid)
@inline Vᶠᶜᶠ(i, j, k, grid) = Azᶠᶜᶠ(i, j, k, grid) * Δzᶠᶜᶠ(i, j, k, grid)
@inline Vᶜᶠᶠ(i, j, k, grid) = Azᶜᶠᶠ(i, j, k, grid) * Δzᶜᶠᶠ(i, j, k, grid)
@inline Vᶠᶠᶠ(i, j, k, grid) = Azᶠᶠᶠ(i, j, k, grid) * Δzᶠᶠᶠ(i, j, k, grid)

#####
##### Generic functions for specified locations
#####
##### For example, Δx(i, j, k, Face, Center, Face) is equivalent to = Δxᶠᶜᶠ(i, j, k, grid).
#####
##### We also use the function "volume" rather than `V`.
#####

location_code(LX, LY, LZ) = Symbol(interpolation_code(LX), interpolation_code(LY), interpolation_code(LZ))

for LX in (:Center, :Face)
    for LY in (:Center, :Face)
        for LZ in (:Center, :Face)
            LXe = @eval $LX
            LYe = @eval $LY
            LZe = @eval $LZ

            volume_function = Symbol(:V, location_code(LXe, LYe, LZe))
            @eval begin
                volume(i, j, k, grid, ::$LX, ::$LY, ::$LZ) = $volume_function(i, j, k, grid)
            end

            for op in (:Δ, :A), dir in (:x, :y, :z)
                func   = Symbol(op, dir)
                metric = Symbol(op, dir, location_code(LXe, LYe, LZe))

                @eval begin
                    $func(i, j, k, grid, ::$LX, ::$LY, ::$LZ) = $metric(i, j, k, grid)
                end
            end
        end
    end
end


text:
# top-level functions to run neutral.jl without parallel map
# Example run:  julia -L neutral.jl nrun.jl examples/sn_example1
# Example run:  julia -L neutral.jl nrun.jl examples/sn_example1 <seed>

function run_trials( simname::AbstractString )
  global num_trials
  global use_population
  global save_populations
  global log_error    # use log normal error:  work in log space
  global wright_fisher_copy
  #println("log_error: ",log_error)
  if !(simtype == 3 || simtype == 4)
    error("simtype must be 3 or 4 for simple neutral evolution!")
  end
  # set defaults for parameters that may not be included in the config file.
  if !isdefined(:use_population)
    use_population = true
  end
  if !isdefined(:save_populations)
    save_populations = false
  end
  if !isdefined(:log_error)
    println("log_error not defined")
    log_error = false    # default is to not use log error
  end
  if !isdefined(:wright_fisher_copy)
    println("wright_fisher_copy not defined")
    wright_fisher_copy = true     # default is to use wright-fisher copy
  end
  sn = simple_neutral_init( simtype, N, mutstddev, ngens, initial_value, num_trials, record_interval, use_population,
      save_populations, log_error, wright_fisher_copy )
  csn = cummulative_neutral_init( sn )
  print_simple_neutral_params( sn )
  for t = 1:num_trials
    if sn.N > 1 || use_population
      sn = simple_neutral_simulation( sn )
    elseif sn.N == 1
      #sn = ces_log( sn )
      sn = ces( sn )
    else
      error(" sn.N must be positive ")
    end
    #println("trial: ",t,"  avg mean: ",sn.average_attr_mean,"  avg coef var: ",sn.average_attr_coef_var )
    #println("mean history: ",sn.attr_mean_history)
    #println("coef_var history: ",sn.attr_coef_var_history)
    accumulate_results( sn, csn )
  end
  #print_cummulative_neutral( csn )
  #print_results( csn )
  #=
  println("saved_populations: ")
  for i = 1:csn.num_trials
    println(csn.saved_populations[:,i])
  end
  =#
  if !csn.save_populations
    writeheader( STDOUT, csn )
    writerows( STDOUT, csn )
    open("$(simname).csv","w") do stream
      writeheader( stream, csn )
      writerows( stream, csn )
    end
  else
    writeheader_populations( STDOUT, csn )
    writerows_populations( STDOUT, csn )
    open("$(simname)_pops.csv","w") do stream
      writeheader_populations( stream, csn )
      writerows_populations( stream, csn )
    end
  end
end

function coef_var( lst )
  return std(lst)/mean(lst)
end

if length(ARGS) == 0
  simname = "examples/sn_example1"
else
  simname = ARGS[1]
  if length(ARGS) >= 2   # second command-line argument is random number seed
    seed = parse(Int,ARGS[2])
    println("seed: ",seed)
    srand(seed)
  end
end
println("simname: ",simname)
include("$(simname).jl")
println("simtype: ",simtype)
run_trials( simname )


text:
########################################################################
# Conditional gradient method
########################################################################
"""
Structure to hold the iterates of the accelerated dual conditional gradient method.

AccDualCGIterable(M, b, A; optTol=√eps)

Construct a dual conditional gradient iterator for the problem

    minimize  ½‖Mx-b‖²  subject to gauge(A|x) ≤ 1.

The absolute tolerance `optTol` is used to determine the stopping condition:

    gap ≤ optTol(1+‖b‖₂)

"""
mutable struct AccDualCGIterable{T1<:Float64, T2<:Vector{T1}, T3<:AbstractArray{T1}, T4<:AbstractLinearOp, T5<:AbstractAtomicSet}
    M::T4           # linear operator
    b::T2           # observation
    A::T5           # atomic set
    r::T2           # residual
    z::T3           # negative gradient
    Ma::T3          # image of the atom
    Mx::T3          # image of primal variable
    My::T3          # image of primal variable 2
    τ::T1           # current scale of atomic set
    optTol::T1      # tolerance
    function AccDualCGIterable(M::AbstractLinearOp, b::Vector{Float64}, A::AbstractAtomicSet, τ::Float64, optTol = sqrt(eps(1.0)))
        r = copy(b)
        z = M'*r
        Ma = M*expose(A, z)
        Mx = copy(Ma); fill!(Mx, 0.0)
        My = copy(Mx)
        new{Float64, Vector{Float64}, AbstractArray{Float64}, AbstractLinearOp, AbstractAtomicSet}(M, b, A, r, z, Ma, Mx, My, τ, optTol)
    end
end

function oracleExit(p::AccDualCGIterable, gap::Float64)
    return gap < p.optTol
end

"""
Defines one complete iteration of the dual conditional gradient method.
"""
function iterate(p::AccDualCGIterable, k::Int=0)
    # set δ
    δ = 2 / (k+3)
    # compute My
    p.My .= (1-δ)*p.Mx + δ*p.Ma
    # compute z
    p.z .*= (1-δ)
    r = p.M'*(p.b - p.My)
    p.z .+= δ*r
    # Compute  M*a  where  a = expose(τA, z)
    a = expose(p.A, p.z)
    p.Ma .= p.M*a
    rmul!(p.Ma, p.τ)
    # check exit
    gap = dot(p.Ma - p.Mx, p.r)
    oracleExit(p, gap) && return (gap, p.r), k+1
    # compute Mx
    p.Mx .*= (1-δ)
    p.Mx .+= δ*p.Ma
    # compute r
    p.r .= p.b - p.Mx

    return (gap, p.r), k+1
end

"""
Return an exit flag for the dual conditional gradient algorithm.
"""
function checkExitCg(gap::Float64, u::Float64, k::Int64, α::Float64, feaTol::Float64, maxIts::Int64)
    u < α + feaTol && return :feasible
    gap ≤ feaTol/2 && return :optimal                   # current lasso problem is fully optimal
    k ≥ maxIts && return :iterations                    # out of iterations
    return :noerror                                     # otherwise
end

"""
Dual conditional gradient method for the problem

    min 1/2||Mx-b||^2 subj to gauge(x | A) ≤ τ
"""
function conditional_graident(M::AbstractLinearOp, b::Vector{Float64}, A::AbstractAtomicSet, τ::Float64;
                   α::Float64 = 0.0,
                   tol::Float64 = 1e-12,
                   gapTol::Float64 = 1+norm(b),
                   maxIts::Int = size(M,2),
                   pr::Bool = true,
                   logger::Bool = true
                   )

    feaTol = tol*(1+norm(b))
    # Structure to hold the iterates
    dcg = DualCGIterable(M, b, A); setRadius!(dcg, τ)
    # Structure to hold the solution
    sol = Solution(M, b, A)

    # head logging
    m, n = size(M)
    logger && logger_head_cg(m, n, feaTol, α, maxIts)

    # Check whether b is zero
    if norm(b) < tol
        return sol
    end

    # ----------------------------------------------------------------
    # min loop
    # ----------------------------------------------------------------
    k = 0
    exitFlag = :noerror
    for (gap, r) in dcg

        exitFlag == :noerror || break
        k = k + 1
        u = norm(r)^2/2

        exitFlag = checkExitCg(gap, u, k, α, feaTol, maxIts)

        # --------------------------------------------------------------
        # Primal recovery.
        # --------------------------------------------------------------
        if pr
            exitFlag = primalrecover!(dcg, sol, α, feaTol, exitFlag)
        end

        # --------------------------------------------------------------
        # Logging and bookkeeping.
        # --------------------------------------------------------------
        logger && logger_level_cg(α, k, u, gap, exitFlag, sol.feas)
    end

    primalrecover!(dcg, sol, α, feaTol, exitFlag)

    # foot logging
    logger&& logger_foot_cg(b, sol.feas, k)

    return sol
end

function logger_head_cg(m, n, feaTol, α, maxIts)
    println()
    println("  -------------------------------------------------------------------------")
    println("  Dual Conditional Gradient Method")
    println("  -------------------------------------------------------------------------")
    @printf("  %-22s %7d %7s %-20s %7d\n","number of variables",n,"","number of constraints",m)
    @printf("  %-22s %7.2e %7s %-20s %7.2e\n","feasibility tolerance",feaTol,"","α",α)
    @printf("  %-22s %7d \n","max iterations",maxIts)
    println("  -------------------------------------------------------------------------")
    @printf("  %8s   %8s   %8s   %8s   %10s\n",
            "iteration","u-α","gap","infeas-α","exitFlag")
end

function logger_foot_cg(b, feas, totIts)
    scale = 1+norm(b)
    normr = sqrt(2*feas)
    println("  -------------------------------------------------------------------------")
    @printf("  %-22s %8.1e \n","residual (abs)",normr)
    @printf("  %-22s %8.1e \n","residual (rel)",normr/scale)
    @printf("  %-22s %8d \n","total iterations",totIts)
    println("  -------------------------------------------------------------------------")
end

function logger_level_cg(α, k, u, gap, exitFlag, feas)
    @printf("  %8d   %8.2e   %8.2e   %8.2e   %s\n",
            k, u-α, gap, feas-α, exitFlag)
end


text:
# Autogenerated wrapper script for libplist_jll for x86_64-apple-darwin
export libplist, plistutil

JLLWrappers.@generate_wrapper_header("libplist")
JLLWrappers.@declare_library_product(libplist, "@rpath/libplist-2.0.3.dylib")
JLLWrappers.@declare_executable_product(plistutil)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libplist,
        "lib/libplist-2.0.3.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        plistutil,
        "bin/plistutil",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module SamIO

abstract SequenceReader

include("ReferenceContigs.jl")
include("BamReader.jl")
include("BedReader.jl")
include("BinningMap.jl")

end # module


text:
using MathOptInterface
const MOI = MathOptInterface

include("permuted_exponential_cone.jl")

MOI.Utilities.@product_of_sets(Zeros, MOI.Zeros)

MOI.Utilities.@product_of_sets(
    Cones,
    MOI.Nonnegatives,
    MOI.SecondOrderCone,
    PermutedExponentialCone
)

MOI.Utilities.@struct_of_constraints_by_set_types(
    ZerosOrNot,
    MOI.Zeros,
    Union{MOI.Nonnegatives,MOI.SecondOrderCone,PermutedExponentialCone},
)

# !!! danger
#     Make sure to use ECOS's definition of pfloat for double and idxint for
#     int. These are _not_ the same across platforms, differ from Julia's
#     Cdouble and Clong, and even differ across patch releases of ECOS :(.

const OptimizerCache = MOI.Utilities.GenericModel{
    pfloat,
    MOI.Utilities.ObjectiveContainer{pfloat},
    MOI.Utilities.VariablesContainer{pfloat},
    ZerosOrNot{pfloat}{
        MOI.Utilities.MatrixOfConstraints{
            pfloat,
            MOI.Utilities.MutableSparseMatrixCSC{
                pfloat,
                idxint,
                MOI.Utilities.ZeroBasedIndexing,
            },
            Vector{pfloat},
            Zeros{pfloat},
        },
        MOI.Utilities.MatrixOfConstraints{
            pfloat,
            MOI.Utilities.MutableSparseMatrixCSC{
                pfloat,
                idxint,
                MOI.Utilities.ZeroBasedIndexing,
            },
            Vector{pfloat},
            Cones{pfloat},
        },
    },
}

mutable struct _Solution
    ret_val::Union{Nothing,idxint}
    primal::Vector{pfloat}
    dual_eq::Vector{pfloat}
    dual_ineq::Vector{pfloat}
    slack::Vector{pfloat}
    slack_eq::Vector{pfloat}
    objective_value::pfloat
    dual_objective_value::pfloat
    solve_time::pfloat
    iter::idxint
end

# The values used by ECOS are from -7 to 10 so -10 should be safe
function _Solution()
    return _Solution(
        nothing,
        pfloat[],
        pfloat[],
        pfloat[],
        pfloat[],
        pfloat[],
        NaN,
        NaN,
        NaN,
        0,
    )
end

mutable struct Optimizer <: MOI.AbstractOptimizer
    zeros::Union{Nothing,Zeros{pfloat}}
    cones::Union{Nothing,Cones{pfloat}}
    sol::_Solution
    silent::Bool
    options::Dict{Symbol,Any}
    function Optimizer(; kwargs...)
        return new(nothing, nothing, _Solution(), false, Dict{Symbol,Any}())
    end
end

function MOI.is_empty(optimizer::Optimizer)
    return optimizer.zeros === nothing && optimizer.cones === nothing
end

function MOI.empty!(optimizer::Optimizer)
    optimizer.zeros = nothing
    optimizer.cones = nothing
    optimizer.sol = _Solution()
    return
end

function MOI.get(::Optimizer, ::MOI.Bridges.ListOfNonstandardBridges)
    return [PermutedExponentialBridge{pfloat}]
end

function _rows(
    optimizer::Optimizer,
    ci::MOI.ConstraintIndex{MOI.VectorAffineFunction{pfloat},MOI.Zeros},
)
    return MOI.Utilities.rows(optimizer.zeros, ci)
end

function _rows(
    optimizer::Optimizer,
    ci::MOI.ConstraintIndex{MOI.VectorAffineFunction{pfloat}},
)
    return MOI.Utilities.rows(optimizer.cones, ci)
end

MOI.get(::Optimizer, ::MOI.SolverName) = "ECOS"

MOI.get(::Optimizer, ::MOI.SolverVersion) = unsafe_string(ECOS_ver())

# MOI.RawOptimizerAttribute

function MOI.supports(::Optimizer, param::MOI.RawOptimizerAttribute)
    return hasfield(settings, Symbol(param.name))
end

function MOI.set(optimizer::Optimizer, param::MOI.RawOptimizerAttribute, value)
    optimizer.options[Symbol(param.name)] = value
    return
end

function MOI.get(optimizer::Optimizer, param::MOI.RawOptimizerAttribute)
    # TODO(odow): This gives a poor error message if the name of the parameter
    # is invalid.
    return optimizer.options[Symbol(param.name)]
end

# MOI.Silent

MOI.supports(::Optimizer, ::MOI.Silent) = true

function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool)
    optimizer.silent = value
    return
end

MOI.get(optimizer::Optimizer, ::MOI.Silent) = optimizer.silent

# MOI.supports

function MOI.supports(
    ::Optimizer,
    ::Union{
        MOI.ObjectiveSense,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{pfloat}},
    },
)
    return true
end

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VectorAffineFunction{pfloat}},
    ::Type{
        <:Union{
            MOI.Zeros,
            MOI.Nonnegatives,
            MOI.SecondOrderCone,
            PermutedExponentialCone,
        },
    },
)
    return true
end

function _optimize!(dest::Optimizer, src::OptimizerCache)
    MOI.empty!(dest)
    Ab = MOI.Utilities.constraints(
        src.constraints,
        MOI.VectorAffineFunction{pfloat},
        MOI.Zeros,
    )
    A = Ab.coefficients
    Gh = MOI.Utilities.constraints(
        src.constraints,
        MOI.VectorAffineFunction{pfloat},
        MOI.Nonnegatives,
    )
    G = Gh.coefficients
    soc_indices = MOI.get(
        Gh,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{pfloat},
            MOI.SecondOrderCone,
        }(),
    )
    q = MOI.dimension.(MOI.get.(Gh, MOI.ConstraintSet(), soc_indices))
    @assert A.n == G.n
    max_sense = MOI.get(src, MOI.ObjectiveSense()) == MOI.MAX_SENSE
    obj =
        MOI.get(src, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{pfloat}}())
    objective_constant = MOI.constant(obj)
    c0 = zeros(A.n)
    for term in obj.terms
        c0[term.variable.value] += term.coefficient
    end
    dest.zeros = deepcopy(Ab.sets) # TODO copy(Ab.sets)
    dest.cones = deepcopy(Gh.sets) # TODO copy(Gh.sets)
    options = copy(dest.options)
    if dest.silent
        options[:verbose] = false
    end
    num_exponential = MOI.get(
        Gh,
        MOI.NumberOfConstraints{
            MOI.VectorAffineFunction{pfloat},
            PermutedExponentialCone,
        }(),
    )
    # This is before ECOS_setup because ECOS modifies the constant vectors
    # in-place!
    slack_eq = copy(Ab.constants)
    inner = ECOS_setup(
        A.n,
        G.m,
        A.m,
        Gh.sets.num_rows[1],
        length(q),
        q,
        num_exponential,
        # Negated because ECOS uses `h - Gx ∈ K`.
        -G.nzval,
        G.colptr,
        G.rowval,
        # Negated because ECOS uses `b -  Ax = 0`.
        -A.nzval,
        A.colptr,
        A.rowval,
        max_sense ? -c0 : c0,
        Gh.constants,
        Ab.constants,
    )
    if inner == C_NULL
        error("ECOS failed to construct problem.")
    end
    unsafe_add_settings(inner, options)
    ret_val = ECOS_solve(inner)
    ecos_prob = unsafe_load(inner)::pwork
    stat = unsafe_load(ecos_prob.info)::stats
    x = copy(unsafe_wrap(Array, ecos_prob.x, ecos_prob.n))
    for col in 1:A.n
        for i in A.colptr[col]:(A.colptr[col+1]-1)
            slack_eq[A.rowval[i+1]+1] += A.nzval[i+1] * x[col]
        end
    end
    dest.sol = _Solution(
        ret_val,
        x,
        copy(unsafe_wrap(Array, ecos_prob.y, ecos_prob.p)),
        copy(unsafe_wrap(Array, ecos_prob.z, ecos_prob.m)),
        copy(unsafe_wrap(Array, ecos_prob.s, ecos_prob.m)),
        slack_eq,
        max_sense ? -stat.pcost : stat.pcost,
        max_sense ? -stat.dcost : stat.dcost,
        stat.tsetup + stat.tsolve,
        stat.iter,
    )
    ECOS_cleanup(inner, 0)
    if !MOI.Utilities.is_ray(MOI.get(dest, MOI.PrimalStatus()))
        dest.sol.objective_value += objective_constant
    end
    if MOI.Utilities.is_ray(MOI.get(dest, MOI.DualStatus()))
        # ECOS can return very large rays here! Without this rescaling it fails
        # many of the MOI tests due to rounding error in floating-point.
        len = sqrt(sum(dest.sol.dual_ineq .^ 2) + sum(dest.sol.dual_eq .^ 2))
        dest.sol.dual_ineq ./= len
        dest.sol.dual_eq ./= len
        dest.sol.dual_objective_value /= len
    else
        dest.sol.dual_objective_value += objective_constant
    end
    return
end

function MOI.optimize!(dest::Optimizer, src::OptimizerCache)
    _optimize!(dest, src)
    return MOI.Utilities.identity_index_map(src), false
end

function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
    cache = OptimizerCache()
    index_map = MOI.copy_to(cache, src)
    _optimize!(dest, cache)
    return index_map, false
end

MOI.get(optimizer::Optimizer, ::MOI.SolveTimeSec) = optimizer.sol.solve_time

function MOI.get(optimizer::Optimizer, ::MOI.BarrierIterations)
    return Int64(optimizer.sol.iter)
end

function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString)
    # Strings from https://github.com/ifa-ethz/ecos/blob/master/include/ecos.h
    flag = optimizer.sol.ret_val
    if flag === nothing
        return "Optimize not called"
    elseif flag == ECOS_OPTIMAL
        return "Problem solved to optimality"
    elseif flag == ECOS_OPTIMAL + ECOS_INACC_OFFSET
        return "Problem solved to inaccurate optimality"
    elseif flag == ECOS_PINF
        return "Found certificate of primal infeasibility"
    elseif flag == ECOS_PINF + ECOS_INACC_OFFSET
        return "Found inaccurate certificate of primal infeasibility"
    elseif flag == ECOS_DINF
        return "Found certificate of dual infeasibility"
    elseif flag == ECOS_DINF + ECOS_INACC_OFFSET
        return "Found inaccurate certificate of dual infeasibility"
    elseif flag == ECOS_MAXIT
        return "Maximum number of iterations reached"
    elseif flag == ECOS_NUMERICS
        return "Search direction unreliable"
    elseif flag == ECOS_OUTCONE
        return "s or z got outside the cone, numerics?"
    elseif flag == ECOS_SIGINT
        return "solver interrupted by a signal/ctrl-c"
    else
        @assert flag == ECOS_FATAL
        return "Unknown problem in solver"
    end
end

# Implements getter for result value and statuses
function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus)
    flag = optimizer.sol.ret_val
    if flag === nothing
        return MOI.OPTIMIZE_NOT_CALLED
    elseif flag == ECOS_OPTIMAL
        return MOI.OPTIMAL
    elseif flag == ECOS_PINF
        return MOI.INFEASIBLE
    elseif flag == ECOS_DINF
        return MOI.DUAL_INFEASIBLE
    elseif flag == ECOS_MAXIT
        return MOI.ITERATION_LIMIT
    elseif flag == ECOS_NUMERICS || flag == ECOS_OUTCONE
        return MOI.NUMERICAL_ERROR
    elseif flag == ECOS_SIGINT
        return MOI.INTERRUPTED
    elseif flag == ECOS_FATAL
        return MOI.OTHER_ERROR
    elseif flag == ECOS_OPTIMAL + ECOS_INACC_OFFSET
        return MOI.ALMOST_OPTIMAL
    elseif flag == ECOS_PINF + ECOS_INACC_OFFSET
        return MOI.ALMOST_INFEASIBLE
    else
        @assert flag == ECOS_DINF + ECOS_INACC_OFFSET
        return MOI.ALMOST_DUAL_INFEASIBLE
    end
end

function MOI.get(optimizer::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.sol.objective_value
end

function MOI.get(optimizer::Optimizer, attr::MOI.DualObjectiveValue)
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.sol.dual_objective_value
end

function MOI.get(optimizer::Optimizer, attr::MOI.PrimalStatus)
    if attr.result_index > MOI.get(optimizer, MOI.ResultCount())
        return MOI.NO_SOLUTION
    end
    flag = optimizer.sol.ret_val
    if flag == ECOS_OPTIMAL
        return MOI.FEASIBLE_POINT
    elseif flag == ECOS_PINF
        return MOI.INFEASIBLE_POINT
    elseif flag == ECOS_DINF
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif flag == ECOS_MAXIT
        return MOI.UNKNOWN_RESULT_STATUS
    elseif flag == ECOS_OPTIMAL + ECOS_INACC_OFFSET
        return MOI.NEARLY_FEASIBLE_POINT
    elseif flag == ECOS_PINF + ECOS_INACC_OFFSET
        return MOI.INFEASIBLE_POINT
    elseif flag == ECOS_DINF + ECOS_INACC_OFFSET
        return MOI.NEARLY_INFEASIBILITY_CERTIFICATE
    else
        return MOI.OTHER_RESULT_STATUS
    end
end

function MOI.get(
    optimizer::Optimizer,
    attr::MOI.VariablePrimal,
    vi::MOI.VariableIndex,
)
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.sol.primal[vi.value]
end

function MOI.get(
    optimizer::Optimizer,
    attr::MOI.ConstraintPrimal,
    ci::MOI.ConstraintIndex{MOI.VectorAffineFunction{pfloat},MOI.Zeros},
)
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.sol.slack_eq[_rows(optimizer, ci)]
end

function MOI.get(
    optimizer::Optimizer,
    attr::MOI.ConstraintPrimal,
    ci::MOI.ConstraintIndex,
)
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.sol.slack[_rows(optimizer, ci)]
end

function MOI.get(optimizer::Optimizer, attr::MOI.DualStatus)
    if attr.result_index > MOI.get(optimizer, MOI.ResultCount())
        return MOI.NO_SOLUTION
    end
    flag = optimizer.sol.ret_val
    if flag == ECOS_OPTIMAL
        return MOI.FEASIBLE_POINT
    elseif flag == ECOS_PINF
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif flag == ECOS_DINF
        return MOI.INFEASIBLE_POINT
    elseif flag == ECOS_MAXIT
        return MOI.UNKNOWN_RESULT_STATUS
    elseif flag == ECOS_OPTIMAL + ECOS_INACC_OFFSET
        return MOI.NEARLY_FEASIBLE_POINT
    elseif flag == ECOS_PINF + ECOS_INACC_OFFSET
        return MOI.NEARLY_INFEASIBILITY_CERTIFICATE
    elseif flag == ECOS_DINF + ECOS_INACC_OFFSET
        return MOI.INFEASIBLE_POINT
    end
    return MOI.OTHER_RESULT_STATUS
end

function MOI.get(
    optimizer::Optimizer,
    attr::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{MOI.VectorAffineFunction{pfloat},MOI.Zeros},
)
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.sol.dual_eq[_rows(optimizer, ci)]
end

function MOI.get(
    optimizer::Optimizer,
    attr::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{MOI.VectorAffineFunction{pfloat}},
)
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.sol.dual_ineq[_rows(optimizer, ci)]
end

MOI.get(optimizer::Optimizer, ::MOI.ResultCount) = 1