text:
"""
`heatmap(z; nargs...)` → `Plot`

Description
============

Draws a heatmap for the given points.
It uses the parameter `z` (which should be a matrix)
as the z value of the heatmap, with the column and row indices of the matrix
as x and y coordinates respectively.

Usage
======

    heatmap(z::AbstractMatrix; title = "", out_stream::Union{Nothing,IO} = nothing, width = 0, height = 0,
            xlabel = "", ylabel = "", zlabel = "", labels = true,
            border = :solid, colormap = :viridis, colorbar_border = :solid, colorbar = :true,
            xfact = 0, yfact = 0, xlim = (0, 0), ylim = (0, 0), zlim = (0, 0),
            xoffset = 0.0, yoffset = 0.0, margin = 3, padding = 1)

Arguments
==========

- **`z`** : The z value for each point.

- **`title`** : Text to display on the top of the plot.

- **`xlabel`** : x-axis label

- **`ylabel`** : y-axis label

- **`zlabel`** : z-axis (colorbar) label

- **`width`** : Number of characters per row that should be used for plotting.

- **`height`** : Number of rows that should be used for plotting. Not applicable to `barplot`.

- **`border`** : The style of the bounding box of the plot. Supports `:solid`, `:bold`, `:dashed`, `:dotted`, `:ascii`, and `:none`.

- **`colormap`** : The colormap to use for the heatmap. Supports `:viridis`, `:plasma`, `:magma`, and `:inferno`.
                   Alternatively, supply a function `f: z, zmin, zmax -> terminal color (Int)`,
                   or a vector of RGB vectors in the format shown [here](https://github.com/BIDS/colormap/blob/master/colormaps.py).

- **`colorbar_border`** : The style of the bounding box of the color bar. Supports `:solid`, `:bold`, `:dashed`, `:dotted`, `:ascii`, and `:none`.

- **`xfact`** : Scale for the x coordinate. Defaults to 0 - i.e. each column in `z` corresponds to one unit, x origin starting at 1. If set to anything else, the x origin will start at 0.

- **`yfact`** : Scale for the y coordinate labels. Defaults to 0 - i.e. each row in `z` corresponds to one unit, y origin starting at 1. If set to anything else, the y origin will start at 0.

- **`xlim`** : Plotting range for the x coordinate (after scaling). `(0, 0)` stands for automatic.

- **`ylim`** : Plotting range for the y coordinate (after scaling). `(0, 0)` stands for automatic.

- **`zlim`** : Data range that the colormap is scaled to. `(0, 0)` stands for automatic.

- **`xoffset`** : Plotting offset for the x coordinate (after scaling).

- **`yoffset`** : Plotting offset for the y coordinate (after scaling).

- **`margin`** : Number of empty characters to the left of the whole plot.

- **`padding`** : Space of the left and right of the plot between the labels and the canvas.

- **`labels`** : Can be used to hide the labels by setting `labels=false`.

- **`fix_ar`** : Fix terminal aspect ratio (experimental)

Returns
========

A plot object of type `Plot{HeatmapCanvas}`

Author(s)
==========

- Rowan Katekar (https://github.com/rjkat)

See also
=========

`Plot`, `scatterplot`, `HeatmapCanvas`
"""
function heatmap(
    z::AbstractMatrix; xlim = (0, 0), ylim = (0, 0), zlim = (0, 0), xoffset = 0., yoffset = 0.,
    out_stream::Union{Nothing,IO} = nothing, width::Int = 0, height::Int = 0, margin::Int = 3,
    padding::Int = 1, colormap=:viridis, xfact=0, yfact=0, labels = true,
    fix_ar::Bool = false, kw...
)
    nrows = size(z, 1)
    ncols = size(z, 2)

    # if scale is auto, use the matrix indices as axis labels
    # otherwise, start axis labels at zero
    X = xfact == 0 ? collect(1:ncols) : collect(0:(ncols-1)) .* xfact
    X .+= xoffset
    xfact = xfact == 0 ? 1 : xfact
    Y = yfact == 0 ? collect(1:nrows) : collect(0:(nrows-1)) .* yfact
    Y .+= yoffset
    yfact = yfact == 0 ? 1 : yfact

    # set the axis limits automatically
    if xlim == (0, 0) && length(X) > 0
        xlim = extrema(X)
    end
    if ylim == (0, 0) && length(Y) > 0
        ylim = extrema(Y)
    end

    # select a subset of z based on the supplied limits
    firstx = findfirst(x -> x >= xlim[1], X)
    lastx = findlast(x -> x <= xlim[2], X)
    xrange = (firstx == nothing || lastx == nothing) ? (1:0) : (firstx:lastx)
    firsty = findfirst(y -> y >= ylim[1], Y)
    lasty = findlast(y -> y <= ylim[2], Y)
    yrange = (firsty == nothing || lasty == nothing) ? (1:0) : (firsty:lasty)
    z = z[yrange, xrange]
    X = X[xrange]
    Y = Y[yrange]

    # allow z to be an array over which min and max is not defined,
    # e.g. an array of RGB color values
    minz, maxz = (0, 0)
    noextrema = true
    try
        minz, maxz = extrema(z)
        noextrema = false
    catch
    end
    if zlim != (0, 0)
        noextrema && throw(ArgumentError("zlim cannot be set when the element type is $(eltype(z))"))
        minz, maxz = zlim
    end

    # if z is an rgb image, translate the colors directly to the terminal
    callback = if length(z) > 0 && all(x -> x in fieldnames(eltype(z)), [:r, :g, :b])
        (z, minz, maxz) -> rgbimgcolor(z)
    else
        colormap_callback(colormap)
    end

    nrows = ceil(Int, (ylim[2] - ylim[1]) / yfact) + 1
    ncols = ceil(Int, (xlim[2] - xlim[1]) / xfact) + 1
    data_ar = ncols / nrows  # data aspect ratio

    max_width = width == 0 ? (height == 0 ? 0 : ceil(Int, height * data_ar)) : width
    max_height = height == 0 ? (width == 0 ? 0 : ceil(Int, width / data_ar)) : height

    # 2nrows: compensate nrows(c::HeatmapCanvas) = div(size(grid(c), 2) + 1, 2)
    width, height, max_width, max_height = get_canvas_dimensions_for_matrix(
        HeatmapCanvas, 2nrows, ncols, max_width, max_height,
        width, height, margin, padding, out_stream, fix_ar
    )

    if height < 7
        # for small plots, don't show colorbar by default
        colorbar = !noextrema && get(kw, :colorbar, false)
    else
        # show colorbar by default, unless set to false, or labels == false
        colorbar = !noextrema && get(kw, :colorbar, labels)
    end
    kw = (; kw..., colorbar=colorbar)

    xs = length(X) > 0 ? [X[1], X[end]] : Float64[0, 0]
    ys = length(Y) > 0 ? [Y[1], Y[end]] : Float64[0, 0]
    plot = Plot(
        xs, ys, HeatmapCanvas;
        grid = false, colormap = callback, colorbar = colorbar, colorbar_lim = (minz, maxz),
        ylim = ylim, xlim = xlim, labels = labels, width = width, height = height,
        min_width = 1, min_height = 1, kw...
    )

    for row = 1:length(Y)
        points!(plot,
            X,
            fill(Y[row], length(X)),
            UserColorType[callback(v, minz, maxz) for v in z[row, :]]
        )
    end
    plot
end


text:
export CUFFTError

struct CUFFTError <: Exception
    code::cufftResult
end

Base.convert(::Type{cufftResult}, err::CUFFTError) = err.code

Base.showerror(io::IO, err::CUFFTError) =
    print(io, "CUFFTError: ", description(err), " (code $(reinterpret(Int32, err.code)), $(name(err)))")

name(err::CUFFTError) = string(err.code)

## COV_EXCL_START
function description(err::CUFFTError)
    if err.code == CUFFT_SUCCESS
        "the operation completed successfully"
    elseif err.code == CUFFT_INVALID_PLAN
        "cuFFT was passed an invalid plan handle"
    elseif err.code == CUFFT_ALLOC_FAILED
        "cuFFT failed to allocate GPU or CPU memory"
    elseif err.code == CUFFT_INVALID_TYPE
        "cuFFT invalid type " # No longer used
    elseif err.code == CUFFT_INVALID_VALUE
        "user specified an invalid pointer or parameter"
    elseif err.code == CUFFT_INTERNAL_ERROR
        "driver or internal cuFFT library error"
    elseif err.code == CUFFT_EXEC_FAILED
        "failed to execute an FFT on the GPU"
    elseif err.code == CUFFT_SETUP_FAILED
        "the cuFFT library failed to initialize"
    elseif err.code == CUFFT_INVALID_SIZE
        "user specified an invalid transform size"
    elseif err.code == CUFFT_UNALIGNED_DATA
        "cuFFT unaligned data" # No longer used
    elseif err.code == CUFFT_INCOMPLETE_PARAMETER_LIST
        "missing parameters in call"
    elseif err.code == CUFFT_INVALID_DEVICE
        "execution of a plan was on different GPU than plan creation"
    elseif err.code == CUFFT_PARSE_ERROR
        "internal plan database error"
    elseif err.code == CUFFT_NO_WORKSPACE
        "no workspace has been provided prior to plan execution"
    elseif err.code == CUFFT_NOT_IMPLEMENTED
        "function does not implement functionality for parameters given."
    elseif err.code == CUFFT_LICENSE_ERROR
        "cuFFT license error" # Used in previous versions.
    elseif err.code == CUFFT_NOT_SUPPORTED
        "operation is not supported for parameters given."
    else
        "no description for this error"
    end
end
## COV_EXCL_STOP

## API call wrapper

# outlined functionality to avoid GC frame allocation
@noinline function throw_api_error(res)
    throw(CUFFTError(res))
end

function initialize_api()
    # make sure the calling thread has an active context
    CUDAnative.initialize_context()
end

macro check(ex)
    quote
        res = @retry_reclaim CUFFT_ALLOC_FAILED $(esc(ex))
        if res != CUFFT_SUCCESS
            throw_api_error(res)
        end

        return
    end
end

ext: jl


text:
"""
Build the data for the model, returning intermediate `DataFrame`s containing parsed data.
The function will not parse data where the corresponding model setting equals `missing`.
"""
function parse_data_to_model!(dtr::AbstractDieterModel; dataname::AbstractString="")

    datapath = dtr.settings[:datapath]

    fileDict = dtr.data["files"]

    dtr.data["dataframes"] = Dict{String,DataFrame}()

    dfDict = dtr.data["dataframes"]

    ## Base data

    dfDict["nodes"] = parse_file(fileDict["nodes"]; dataname=dataname)
    parse_nodes!(dtr)

    # e.g. fileDict["tech"] = joinpath(datapath,"base","technologies.csv")
    dfDict["tech"] = parse_file(fileDict["tech"]; dataname=dataname)
    parse_base_technologies!(dtr)

    # e.g. fileDict["storage"] = joinpath(datapath,"base","storages.csv")
    dfDict["storage"] = parse_file(fileDict["storage"]; dataname=dataname)
    parse_storages!(dtr)

    # e.g. fileDict["load"] = joinpath(datapath,"base","load.csv")
    dfDict["load"] = parse_file(fileDict["load"]; dataname=dataname)
    parse_load!(dtr, dfDict["load"])

    # e.g. fileDict["avail"] = joinpath(datapath,"base","availability.csv")
    dfDict["avail"] = parse_file(fileDict["avail"]; dataname=dataname)
    parse_availibility!(dtr)

    # Relations (i.e set-to-set correspondences)
    dfDict["map_node_demand"] = parse_file(fileDict["map_node_demand"]; dataname=dataname)
    dfDict["map_node_tech"] = parse_file(fileDict["map_node_tech"]; dataname=dataname)
    dfDict["map_node_storages"] = parse_file(fileDict["map_node_storages"]; dataname=dataname)
    dfDict["arcs"] = parse_file(fileDict["arcs"]; dataname=dataname)
    initialise_set_relation_data!(dtr)
    parse_set_relations!(dtr)
    STABLE.parse_arcs!(dtr)

    # Calculated base parameters
    calc_base_parameters!(dtr)

    return dfDict

end

function parse_extensions!(dtr::AbstractDieterModel; dataname::AbstractString="")

    datapath = dtr.settings[:datapath]
    fileDict = dtr.data["files"]
    dfDict = dtr.data["dataframes"]

    ## EV
    if !(dtr.settings[:ev] |> ismissing)
        # e.g. fileDict["ev"] = joinpath(datapath,"ev","ev.csv")
        dfDict["ev"] = parse_file(fileDict["ev"]; dataname=dataname)
        parse_ev_technologies!(dtr, dfDict["ev"])

        # e.g. fileDict["ev_demand"] = joinpath(datapath,"ev","ev_demand.csv")
        dfDict["ev_demand"] = parse_file(fileDict["ev_demand"]; dataname=dataname)
        parse_ev_demand!(dtr, dfDict["ev_demand"])

        # e.g. fileDict["ev_power"] = joinpath(datapath,"ev","ev_power.csv")
        dfDict["ev_power"] = parse_file(fileDict["ev_power"]; dataname=dataname)
        parse_ev_power!(dtr, dfDict["ev_power"])

        calc_ev_quantity!(dtr)
    else
        @info "No Electric Vehicle aspects present in model."
        dtr.sets[:ElectricVehicles] = Array([])
        for s in [:AbsoluteEvDemand,:AbsoluteEvPower,:AbsoluteEvCapacity,:EvFuel,:EvType]
            dtr.parameters[s] = Dict()
        end
    end

    ## Heat
    if !(dtr.settings[:heat] |> ismissing)
        # e.g. fileDict["heat"] = joinpath(datapath,"heat","heat.csv")
        dfDict["heat"] = parse_file(fileDict["heat"]; dataname=dataname)
        parse_heat!(dtr, dfDict["heat"])

        # e.g. fileDict["heat_technologies"] = joinpath(datapath,"heat","heat_technologies.csv")
        dfDict["heat_technologies"] = parse_file(fileDict["heat_technologies"]; dataname=dataname)
        parse_heat_technologies!(dtr, dfDict["heat_technologies"])

        # e.g. fileDict["buildings"] = joinpath(datapath,"heat","buildings.csv")
        dfDict["buildings"] = parse_file(fileDict["buildings"]; dataname=dataname)
        parse_buildings!(dtr, dfDict["buildings"])

        # e.g. fileDict["temperature"] = joinpath(datapath,"heat","temperature.csv")
        dfDict["temperature"] = parse_file(fileDict["temperature"]; dataname=dataname)
        parse_temperature!(dtr, dfDict["temperature"])

        # e.g. fileDict["heat_demand"] = joinpath(datapath,"heat","heat_demand.csv")
        dfDict["heat_demand"] = parse_file(fileDict["heat_demand"]; dataname=dataname)
        parse_heat_demand!(dtr, dfDict["heat_demand"])

        # e.g. fileDict["dhw_demand"] = joinpath(datapath,"heat","dhw_demand.csv")
        dfDict["dhw_demand"] = parse_file(fileDict["dhw_demand"]; dataname=dataname)
        parse_dhw_demand!(dtr, dfDict["dhw_demand"])

        calc_hp_cop!(dtr)
        calc_heat_demand!(dtr)
    else
        @info "No Heat aspects present in model."
        dtr.sets[:BuildingType] = Array([])
        dtr.sets[:HeatingType]  = Array([])
        for s in [:HeatConsumption, :HeatShare, :MaxLevel, :HeatMaxPower, :StaticEfficiency, :CoP]
            dtr.parameters[s] = Dict()
        end
    end

    ## Hydrogen (h2)
    if !(dtr.settings[:h2] |> ismissing)
        # e.g. fileDict["h2"] = joinpath(datapath,"h2","h2_technologies.csv")
        dfDict["h2_technologies"] = parse_file(fileDict["h2_technologies"]; dataname=dataname)

        # Overwrite values from "h2_technologies" table by Scenario:
        ScenCostPower = dtr.parameters[:ScenCostPower]

        tech_types = dtr.settings[:tech_types]

        H2ElectrolyserTypes = tech_types[:H2ElectrolyserTypes]
        H2RecipEngTypes = tech_types[:H2RecipEngTypes]

        ## Map OvernightCost data into table
        dfDict["h2_technologies"] = @eachrow dfDict["h2_technologies"] begin
            if :H2Technologies in keys(ScenCostPower)
                :H2OvernightCost = ScenCostPower[:H2Technologies]
                h2_electrolyser = :H2Technologies
                region = :Region
                @info "H2 Technology ($(h2_electrolyser)) cost was overwritten in $(region) with scenario value $(ScenCostPower[h2_electrolyser])."
            end
        end
        # else
        #     @warn "H2 Electrolyser ($(h2_electrolyser)) cost was NOT overwritten by scenario value"
        # end

        parse_h2_technologies!(dtr, dfDict["h2_technologies"])

        # # H2 Reciprocating Engine Fuel prices
        InstanceYear = dtr.settings[:inst_year]
        InstanceYear_Sym = Symbol("FYE$(InstanceYear)")
        node2DemReg = dtr.parameters[:node_demreg_dict]

        fuel_prices_sql_db_path = dtr.settings[:datapaths]["Fuel_Price_SQL_db"]
        H2_Cost_Tech_Map = dtr.settings[:h2_cost_tech_map]

        if dtr.settings[:H2_G2P_Allowed_flag] == true
            fileDict["fuel_prices_h2"] = joinpath(datapath,"base","fuel_prices_h2.sql")
            dfDict["fuel_prices_h2"] = parse_file(fileDict["fuel_prices_h2"]; dataname=fuel_prices_sql_db_path)

            df_fuel_h2 = @subset(dfDict["fuel_prices_h2"],:TechID .== H2_Cost_Tech_Map["N_RecipH2"])
            # Rename the InstanceYear column as `:FuelCost`:
            select!(df_fuel_h2,:Region, (InstanceYear_Sym => :FuelCost))
            h2_fuel_dict = Dict(eachrow(df_fuel_h2))

            for ((n,t),v) in dtr.parameters[:FuelCost]
                    if t == H2_Cost_Tech_Map["N_RecipH2"]
                        dtr.parameters[:FuelCost][(n,t)] = h2_fuel_dict[node2DemReg[n]]
                    end
            end
        end

        calc_inv_gas!(dtr)
    else
        @info "No Power-to-Gas/Hydrogen aspects present in model."
        dtr.sets[:H2Technologies] = Array([])
        dtr.sets[:P2G] = Array([])
        dtr.sets[:G2P] = Array([])
        dtr.sets[:GasStorages] = Array([])
        dtr.sets[:Nodes_H2Tech] = Array([])
        dtr.sets[:Nodes_P2G] = Array([])
        dtr.sets[:Nodes_G2P] = Array([])
        dtr.sets[:Nodes_GasStorages] = Array([])

        dtr.sets[:Nodes_P2G_PEM] = Array([])
        dtr.sets[:Nodes_P2G_AE]  = Array([])

        dtr.parameters[:H2Conversion] = Dict{Any,Any}()
        # for s in []
        #     dtr.parameters[s] = Dict()
        # end
    end

    return nothing
end

# Expected input e.g.  df = dtr.data["dataframes"]["nodes"]
# Columns: | Nodes │ NodeType │ NodePromote │ DemandRegion | IncludeLevel |
# Expected that NodeType has (at least) entries including ["TxZone", "REZone"]
function parse_nodes!(dtr::DieterModel)

    df = dtr.data["dataframes"]["nodes"]

    dtr.sets[:Nodes] = disallowmissing(unique(df[!, :Nodes]))
    dtr.sets[:NodeType] = disallowmissing(unique(df[!, :NodeType]))
    dtr.sets[:DemandRegions] = collect(skipmissing(unique(df[!, :DemandRegion])))
    # dtr.sets[:Regions] = disallowmissing(unique(df[!, :Region]))

    # Obtain the Transmission Zones
    dtr.sets[:TxZones] = (@subset(df,:NodeType .== "TxZone"))[!,:Nodes]
    # Obtain the Renewable Energy Zones
    dtr.sets[:REZones] = (@subset(df,:NodeType .== "REZone"))[!,:Nodes]

    # Equate Demand Zones to Transmission Zones
    dtr.sets[:DemandZones] = dtr.sets[:TxZones]

    # Provide a dict (lookup map) between nodes and the associated demand region
    dtr.parameters[:node_demreg_dict] = Dict(zip(df[!,:Nodes],df[!,:DemandRegion]))

    # Provide a dict (lookup map) between nodes and the associated "IncludeLevel", a relative measure of cost increase:
    dtr.parameters[:include_level_dict] = Dict(zip(df[!,:Nodes],df[!,:IncludeLevel]))

    return nothing
end

"""
Expected Dieter model parameters read into `dtr.parameters`:
CapacityFactor, FuelCost, FixedCost, VariableCost, OvernightCost,
CurtailmentCost, LoadIncreaseCost, LoadDecreaseCost, Lifetime,
MaxCapacity, MaxEnergy,
TechTypeCategory, Renewable, Dispatchable,
FuelSource, Efficiency, CarbonContent, CO2_price
"""
function parse_base_technologies!(dtr::DieterModel)

    df = dtr.data["dataframes"]["tech"]

    dtr.sets[:Technologies] = disallowmissing(unique(df[!, :Technologies]))

    dtr.sets[:Renewables] = disallowmissing(unique([row[:Technologies] for row in eachrow(df)
        if row[:Renewable] == 1]))

    dtr.sets[:Conventional] = disallowmissing(unique([row[:Technologies] for row in eachrow(df)
        if row[:Renewable] == 0]))

    dtr.sets[:Dispatchable] = disallowmissing(unique([row[:Technologies]
        for row in eachrow(df) if row[:Dispatchable] == 1]))

    dtr.sets[:NonDispatchable] = disallowmissing(unique([row[:Technologies]
        for row in eachrow(df) if row[:Dispatchable] == 0]))

    # params = map_idcol(df, skip_cols=[:Region, :Renewable, :Dispatchable])
    params = map_idcol(df, [:Region, :Technologies], skip_cols=Symbol[])
    merge!(dtr.parameters, params)

    # If a column of parameters is entirely missing, convert container to a Union type.
    # Otherwise, we get a parameter Dict of type Dict{Tuple{R,T},Missing} :
    DefaultType = Float64 # default type to contain general parameters
    R = typeof(df[1,:Region])
    T = typeof(df[1,:Technologies])

    for (k,v) in dtr.parameters
        if typeof(v) <: Dict{Tuple{R,T},Missing}
            dtr.parameters[k] = convert(Dict{Tuple{R,T},Union{DefaultType,Missing}},v)
        end
    end

    return nothing
end

function parse_storages!(dtr::DieterModel)

    df = dtr.data["dataframes"]["storage"]

    dtr.sets[:Storages] = disallowmissing(unique(df[!,:Storages]))

    params = map_idcol(df, [:Region, :Storages], skip_cols=Symbol[])
    for (k,v) in params
        try
            update_dict!(dtr.parameters, k, v)
        catch e
            print(e)
            error(
                "Update failed: dtr.parameters[$k] is type
                    $(typeof(dtr.parameters[k])) \n while new parameter values are of type
                    $(typeof(v))."
            )
        end
    end

    return nothing
end

function parse_load!(dtr::DieterModel, df::DataFrame)

    params = map_idcol(df, [:Nodes, :TimeIndex], skip_cols=Symbol[])
    # params = map_idcol(df, [:DemandRegion, :TimeIndex], skip_cols=Symbol[])
    for (k,v) in params update_dict!(dtr.parameters, k, v) end

    return nothing
end

function parse_availibility!(dtr::DieterModel)

    df = dtr.data["dataframes"]["avail"]

    params = map_idcol(df, [:RenewRegionID, :TechTypeID, :TimeIndex], skip_cols=Symbol[])
    for (k,v) in params update_dict!(dtr.parameters, k, v) end

    # Construct the set of node/technology pairs that have availability traces.
    # This is used to validate the data before constraint construction.
    availpairs = unique(df[!,[:RenewRegionID,:TechTypeID]])
    dtr.sets[:Nodes_Avail_Techs] = [Tuple(x) for x in eachrow(availpairs)]

    if !isempty( setdiff(dtr.sets[:Nodes_Avail_Techs],dtr.sets[:Nodes_Techs]) )
          println(setdiff(dtr.sets[:Nodes_Avail_Techs],dtr.sets[:Nodes_Techs]))
          error("Available trace data includes a technology an excluded region.")
    end

    return nothing
end

"Create the mathematical set relations for doing set-correspondences in models."
function initialise_set_relation_data!(dtr)
# See util.jl/create_relation(...) for the function definition
# See struct.jl/parse_data_to_model!(...) for construction of the Dict dtr.data["maps"]

# map_<...> is a a reference to a DataFrame containing side-by-side columns of set correspondences
# rel_<...> is a filter function created from the DataFrame map_<...>
# rel_<...> applied is a Julia array object with Tuples containing the mathematical relation for later iteration.

    # Input for this scope (the inputs need to have been parsed previously using initialise_data_file_dict! & parse_data_to_model! )
    dfDict = dtr.data["dataframes"]

    # Output stored in model's data field:
    dtr.data["relations"] = Dict{String,Any}()

    # Relation between a node and associated technologies
    rel_node_tech = create_relation(dfDict["map_node_tech"],:Nodes,:Technologies,:IncludeFlag)
    # rel_node_tech = filter(x->f_node_tech(x[1],x[2]), [(n,t) for n in N for t in T] )

    # Relation between a node and associated storage technologies
    rel_node_storages = create_relation(dfDict["map_node_storages"],:Nodes,:Technologies,:IncludeFlag)
    # rel_node_storages = filter(x->f_node_storages(x[1],x[2]),[(n,s) for n in N for s in S])

    # Relation between a node and associated type
    rel_node_type = create_relation(dfDict["nodes"],:Nodes,:NodeType,:IncludeLevel)

    # Relation between a node and associated type above (its promotion)
    rel_node_promote = create_relation(dfDict["nodes"],:Nodes,:NodePromote,:IncludeLevel)

    # Relation between a node and demand
    rel_node_demand = create_relation(dfDict["map_node_demand"],:Nodes,:DemandZone,:IncludeFlag)
    # rel_node_demand = create_relation(dfDict["map_node_demand"],:Nodes,:DemandRegion,:IncludeFlag)

    # Relation between a node and incident/connected nodes (inter-connectors)
    rel_node_incidence = create_relation(dfDict["arcs"],:FromNode,:ToNode,:TransferCapacity)

    # Store the results:
    dtr.data["relations"]["rel_node_tech"] = rel_node_tech
    dtr.data["relations"]["rel_node_storages"] = rel_node_storages
    dtr.data["relations"]["rel_node_type"] = rel_node_type
    dtr.data["relations"]["rel_node_promote"] = rel_node_promote
    dtr.data["relations"]["rel_node_demand"] = rel_node_demand
    dtr.data["relations"]["rel_node_incidence"] = rel_node_incidence

    return nothing

end

function parse_set_relations!(dtr)

    Nodes = dtr.sets[:Nodes]
    NodeType = dtr.sets[:NodeType]
    DemandRegions = dtr.sets[:DemandRegions]
    DemandZones = dtr.sets[:DemandZones]

    Technologies = dtr.sets[:Technologies]
    Storages = dtr.sets[:Storages]

    Renewables = dtr.sets[:Renewables]
    Conventional = dtr.sets[:Conventional]
    Dispatchable = dtr.sets[:Dispatchable]
    NonDispatchable = dtr.sets[:NonDispatchable]

    # Relation between a node and associated technologies
    rel_node_tech = dtr.data["relations"]["rel_node_tech"]
    dtr.sets[:Nodes_Techs] = tuple2_filter(rel_node_tech, Nodes, Technologies)
    dtr.sets[:Nodes_Renew] = tuple2_filter(rel_node_tech, Nodes, Renewables)
    dtr.sets[:Nodes_Conven] = tuple2_filter(rel_node_tech, Nodes, Conventional)
    dtr.sets[:Nodes_Dispatch] = tuple2_filter(rel_node_tech, Nodes, Dispatchable)
    dtr.sets[:Nodes_NonDispatch] = tuple2_filter(rel_node_tech, Nodes, NonDispatchable)

    # Relation between a node and associated storage technologies
    rel_node_storages = dtr.data["relations"]["rel_node_storages"]
    dtr.sets[:Nodes_Storages] = tuple2_filter(rel_node_storages, Nodes, Storages)

    # Relation between a node and associated type
    rel_node_type = dtr.data["relations"]["rel_node_type"]
    dtr.sets[:Nodes_Types] = tuple2_filter(rel_node_type, Nodes, NodeType)

    # Obtain the Transmission Zones
    # dtr.sets[:TxZones] = [x[1] for x in filter(x -> (x[2]=="TxZone"),dtr.sets[:Nodes_Types])]
    # Obtain the Renewable Energy Zones
    # dtr.sets[:REZones] = [x[1] for x in filter(x -> (x[2]=="REZone"),dtr.sets[:Nodes_Types])]

    # Relation between a node and associated node above it (its "promotion") if that is also a Node
    # Typically, creates a relation of REZone Nodes to the connected TxZone.
    rel_node_promote = dtr.data["relations"]["rel_node_promote"]
    dtr.sets[:Nodes_Promotes] = tuple2_filter(rel_node_promote, Nodes, Nodes)

    # Relation between a node and associated demand
    rel_node_demand = dtr.data["relations"]["rel_node_demand"]
    dtr.sets[:Nodes_Demand] = tuple2_filter(rel_node_demand, Nodes, DemandZones)
    # dtr.sets[:Nodes_Demand] = tuple2_filter(rel_node_demand, Nodes, DemandRegions)

    return nothing

end

function parse_arcs!(dtr::DieterModel)

    df = dtr.data["dataframes"]["arcs"]

    Nodes = dtr.sets[:Nodes]

    # Relation between a node and incident/connected nodes (inter-connectors)
    rel_node_incidence = dtr.data["relations"]["rel_node_incidence"]
    dtr.sets[:Arcs] = tuple2_filter(rel_node_incidence, Nodes, Nodes)

    dtr.sets[:Arcs_From] = disallowmissing(unique(df[!,:FromNode]))

    dtr.sets[:Arcs_REZones] = [(rez,txz) for (rez,txz) in dtr.sets[:Nodes_Promotes]
                                if rez in dtr.sets[:REZones]]

    params = map_idcol(df, [:FromNode, :ToNode], skip_cols=Symbol[])
    for (k,v) in params update_dict!(dtr.parameters, k, v) end

    return nothing
end

function calc_inv_trans_exp!(dtr::DieterModel)

    tec = dtr.parameters[:TransExpansionCost]
    lt_Tx = dtr.settings[:lifetime_Tx]
    i = dtr.settings[:interest]

    dict = Dict((from,to) => tec[from,to]*annuity(i, lt_Tx) for (from,to) in keys(tec))
    update_dict!(dtr.parameters, :InvestmentCostTransExp, dict)

    return nothing
end

annuity(i,lifetime) = i*((1+i)^lifetime) / (((1+i)^lifetime)-1)
# Equivalent alternative: annuity(i,lifetime) = i +  i/( ((1+i)^lifetime)-1 )

"Calculate the amortised investment cost of technologies from capital costs."
function calc_inv_tech!(dtr::DieterModel)
    Nodes_Techs = dtr.sets[:Nodes_Techs]
    # Nodes = dtr.sets[:Nodes]
    # Technologies = dtr.sets[:Technologies]

    oc = dtr.parameters[:OvernightCostPower]
    lt = dtr.parameters[:Lifetime]
    i = dtr.settings[:interest]

    dict = Dict((n,t) => oc[n,t]*annuity(i, lt[n,t]) for (n,t) in Nodes_Techs)
    update_dict!(dtr.parameters, :InvestmentCost, dict)

    return nothing
end

"Calculate the amortised investment cost of technologies from capital costs."
function calc_inv_storages!(dtr::DieterModel)
    Nodes_Storages = dtr.sets[:Nodes_Storages]
    # Nodes = dtr.sets[:Nodes]
    # Storages = dtr.sets[:Storages]

    oce = dtr.parameters[:OvernightCostEnergy]
    ocp = dtr.parameters[:OvernightCostPower]
    lt = dtr.parameters[:Lifetime]
    i = dtr.settings[:interest]

    cost_energy = Dict((n,s) => oce[n,s]*annuity(i, lt[n,s]) for (n,s) in Nodes_Storages)
    cost_power  = Dict((n,s) => ocp[n,s]*annuity(i, lt[n,s]) for (n,s) in Nodes_Storages)

    update_dict!(dtr.parameters, :InvestmentCostEnergy, cost_energy)
    update_dict!(dtr.parameters, :InvestmentCostPower, cost_power)

    return nothing
end

"Calculated marginal costs for plants including variable costs"
function calc_mc!(dtr::DieterModel)
    Nodes_Techs = dtr.sets[:Nodes_Techs]
    Nodes_Storages = dtr.sets[:Nodes_Storages]
    # Nodes = dtr.sets[:Nodes]
    # Technologies = dtr.sets[:Technologies]
    # Storages = dtr.sets[:Storages]

    fc = dtr.parameters[:FuelCost]      # FuelCost in currency/MWh-thermal-input
    eff = dtr.parameters[:Efficiency]   # Efficiency is unitless (MWh-output/MWh-thermal-input)
    cc = dtr.parameters[:CarbonContent] # CarbonContent in units of t-CO2/MWh-thermal-input
    co2 = dtr.settings[:co2]            # Price on carbon in currency/t-CO2
    vc = dtr.parameters[:VariableCost]  # VariableCost in currency/MWh-output

    marginalcost = Dict((n,t) => fc[n,t]/eff[n,t] + (cc[n,t]*co2)/eff[n,t] + vc[n,t] for (n,t) in Nodes_Techs)
    update_dict!(dtr.parameters, :MarginalCost, marginalcost)

    marginalcost_sto = Dict((n,s) => fc[n,s]/eff[n,s] + vc[n,s] for (n,s) in Nodes_Storages)
    update_dict!(dtr.parameters, :MarginalCost, marginalcost_sto)

    return nothing
end

function calc_base_parameters!(dtr::DieterModel)
    calc_inv_tech!(dtr)
    calc_inv_storages!(dtr)
    calc_mc!(dtr)

    return nothing
end


text:
#######################
# MathAST Definitions

abstract type MathExpressionFragment end

struct MathSequence <:MathExpressionFragment
    fragments::Vector{MathExpressionFragment}
end
Base.length(ms::MathSequence) = length(ms.fragments)

const MathSymbolInner = Union{Number, Symbol}
struct MathSymbol{T<:MathSymbolInner} <: MathExpressionFragment
    sym::T
end
MathSymbol(s::AbstractString) = MathSymbol(Symbol(s))

struct MathFrac <: MathExpressionFragment
    numer::MathExpressionFragment
    denom::MathExpressionFragment
end

struct MathSuperscript <: MathExpressionFragment
    base::MathExpressionFragment
    expo::MathExpressionFragment
end

struct MathSubscript <: MathExpressionFragment
    base::Symbol
    sub::Symbol
end

struct MathFunctionCall <: MathExpressionFragment
    fname::Symbol
    arg::MathSequence
end

struct MathMatrix <: MathExpressionFragment
    mat::Matrix{MathExpressionFragment}
end


text:
const rand_str_length = 8

make_random_label() = "ARGNode_$(randstring(rand_str_length))"
make_random_label(str::String) = str * "_" * randstring(rand_str_length)
make_random_label(i::Int) = "ARGNode_$(i)"
make_random_label(str, i) = str * "_" * i

function check(n::ARGNode; strict=true)
	#General topology
	## No ancestor ==> tau is missing
	@assert mapreduce(
		(a,τ) -> !isnothing(a) || (isnothing(a) && ismissing(τ)),
		&,
		n.anc, n.tau
	) "No ancestor and non-missing `tau` field (node $(n.label))"

	# Colors
	## color[i] <==> (root[i] or (ancestor[i] and the anc has color `i` too))
	## (i.e. not color --> root is false and ancestor is nothing)
	for (i,c) in enumerate(color(n))
		if c
			@assert isroot(n, i) || !isnothing(ancestor(n, i)) && hascolor(ancestor(n,i), i) "Inconsistent color / root / ancestor for $(label(n))"
		end
	end

	# If it's a leaf, it should have all colors
	if isleaf(n)
		@assert degree(n) == length(color(n))
	end

end

function Base.show(io::IO, n::ARGNode)
    if !get(io, :compact, false)
        nodeinfo(io, n)
    end
end
Base.show(n::ARGNode) = show(stdout, n)

"""
    nodeinfo(io, node)

Print information about `node`.
"""
function nodeinfo(io, node)
    print(io, "$(label(node)) ")
  	if node.hybrid
  		print(io, "(hybrid) ")
  	end
  	println(io, ":")
  	println(io, "Color: $(color(node))")
  	println(io, "Ancestors: $([label(x) for x in ancestors(node)])")
  	println(io, "Branch length: $(tau(node))")
  	println(io, "Children: $([label(x) for x in children(node)])")
end

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/TreeKnit

max_forks_repo_name: PierreBarrat/TreeKnit


text:
@testset "ACEtb.SlaterKoster Basics" begin

   using ACEtb.SlaterKoster, Test
   using ACEtb.SlaterKoster: get_l, get_bidx

   function test_fail(f)
      try
         f()
         @info("We should have thrown an error here")
         return false
      catch
         return true
      end
   end

   @info("Testing `get_l` for orbitals")
   println(@test get_l(Val{:s}()) == 0)
   println(@test get_l(Val{:p}()) == 1)
   println(@test get_l(Val{:d}()) == 2)
   println(@test get_l(Val{:f}()) == 3)
   println(@test get_l(Val{:g}()) == 4)

   @info("Testing `get_l, get_bidx` for bonds")
   println(@test get_l(skb"spσ") == (0, 1))
   println(@test get_bidx(skb"spσ") == 0)
   println(@test get_l(skb"2p2pπ") == (1, 1))
   println(@test get_bidx(skb"2p2pπ") == 1)
   println(@test get_l(skb"1s2pσ") == (0,1))
   println(@test get_bidx(skb"1s2pσ") == 0)
end


text:
module Monadics

export mbind, mreturn, mjoin, ↣, @monadic

function mbind end
function mreturn end
function mjoin end

macro monadic(x)
    quote
        $(esc(x))
        $(esc(_dispatch_monadic(x)))
    end
end

function _dispatch_monadic(x::Expr)
    _dispatch_monadic(Val(x.head), x.args)
end

function _dispatch_monadic(s::Symbol)
    quote
        @inline mreturn(::Type{$s}, x) = $s(x)
        @inline mjoin(m::$s) = mbind(identity, m)
        @inline Base.map(f::Function, m::$s) =
            mbind(m) do x
                mreturn($s, f(x))
            end
        @inline ↣(m::$s, f::Function) = mbind(f, m)
    end
end

function _dispatch_monadic(::Val{:abstract}, s::Symbol)
    quote
        @inline mreturn(::Type{T}, x) where {T<:$s} = T(x)
        @inline mjoin(m::T) where {T<:$s} = mbind(identity, m)
        @inline Base.map(f::Function, m::T) where {T<:$s} =
            mbind(m) do x
                mreturn(T, f(x))
            end
        @inline ↣(m::T, f::Function) where {T<:$s} = mbind(f, m)
    end
end

function _dispatch_monadic(::Val{:struct}, args::Vector{Any})
    _dispatch_monadic(args[2])
end

function _dispatch_monadic(::Val{:mutable}, args::Vector{Any})
    _dispatch_monadic(args[2])
end

function _dispatch_monadic(::Val{:curly}, args::Vector{Any})
    _dispatch_monadic(args[1])
end

function _dispatch_monadic(::Val{:const}, args::Vector{Any})
    _dispatch_monadic(args[1])
end

function _dispatch_monadic(::Val{:abstract}, args::Vector{Any})
    _dispatch_monadic(Val(:abstract), args[1])
end

function _dispatch_monadic(::Val{:(=)}, args::Vector{Any})
    _dispatch_monadic(args[1])
end

@monadic const Optional{T} = Union{Nothing,T}
@monadic const Maybe{T} = Optional{Some{T}}

end # module


text:
immutable PointWithAbstract
x::AbstractFloat
y::AbstractFloat
end

max_stars_repo_path: Module 2/Chapter03/PointWithAbstract.jl

max_issues_repo_path: Module 2/Chapter03/PointWithAbstract.jl

max_forks_repo_path: Module 2/Chapter03/PointWithAbstract.jl


text:
#=
Generate MarkovChain and DiscreteDP instances randomly.

@author : Daisuke Oyama

=#
import StatsBase: sample
import QuantEcon: MarkovChain, DiscreteDP

# random_markov_chain

"""
Return a randomly sampled MarkovChain instance with n states.

##### Arguments

- `n::Integer` : Number of states.

##### Returns

- `mc::MarkovChain` : MarkovChain instance.

##### Examples

```julia
julia> using QuantEcon

julia> mc = random_markov_chain(3)
Discrete Markov Chain
stochastic matrix:
3x3 Array{Float64,2}:
 0.281188  0.61799   0.100822
 0.144461  0.848179  0.0073594
 0.360115  0.323973  0.315912

```

"""
function random_markov_chain(n::Integer)
    p = random_stochastic_matrix(n)
    mc = MarkovChain(p)
    return mc
end

"""
Return a randomly sampled MarkovChain instance with n states, where each state
has k states with positive transition probability.

##### Arguments

- `n::Integer` : Number of states.

##### Returns

- `mc::MarkovChain` : MarkovChain instance.

##### Examples

```julia
julia> using QuantEcon

julia> mc = random_markov_chain(3, 2)
Discrete Markov Chain
stochastic matrix:
3x3 Array{Float64,2}:
 0.369124  0.0       0.630876
 0.519035  0.480965  0.0
 0.0       0.744614  0.255386

```

"""
function random_markov_chain(n::Integer, k::Integer)
    p = random_stochastic_matrix(n, k)
    mc = MarkovChain(p)
    return mc
end

# random_stochastic_matrix

"""
Return a randomly sampled n x n stochastic matrix with k nonzero entries for
each row.

##### Arguments

- `n::Integer` : Number of states.
- `k::Union{Integer, Void}(nothing)` : Number of nonzero entries in each
column of the matrix. Set to n if note specified.

##### Returns

- `p::Array` : Stochastic matrix.

"""
function random_stochastic_matrix(n::Integer, k::Union{Integer, Void}=nothing)
    if !(n > 0)
        throw(ArgumentError("n must be a positive integer"))
    end
    if k != nothing && !(k > 0 && k <= n)
        throw(ArgumentError("k must be an integer with 0 < k <= n"))
    end

    p = _random_stochastic_matrix(n, n, k=k)

    return transpose(p)
end

"""
Generate a "non-square column stochstic matrix" of shape (n, m), which contains
as columns m probability vectors of length n with k nonzero entries.

##### Arguments

- `n::Integer` : Number of states.
- `m::Integer` : Number of probability vectors.
- `;k::Union{Integer, Void}(nothing)` : Number of nonzero entries in each
column of the matrix. Set to n if note specified.

##### Returns

- `p::Array` : Array of shape (n, m) containing m probability vectors of length
n as columns.

"""
function _random_stochastic_matrix(n::Integer, m::Integer;
                                   k::Union{Integer, Void}=nothing)
    if k == nothing
        k = n
    end
    probvecs = random_probvec(k, m)

    k == n && return probvecs

    # if k < n
    # Randomly sample row indices for each column for nonzero values
    row_indices = Vector{Int}(k*m)
    for j in 1:m
        row_indices[(j-1)*k+1:j*k] = sample(1:n, k, replace=false)
    end

    p = zeros(n, m)
    for j in 1:m
        for i in 1:k
            p[row_indices[(j-1)*k+i], j] = probvecs[i, j]
        end
    end

    return p
end

# random_discrete_dp

"""
Generate a DiscreteDP randomly. The reward values are drawn from the normal
distribution with mean 0 and standard deviation `scale`.

##### Arguments

- `num_states::Integer` : Number of states.
- `num_actions::Integer` : Number of actions.
- `beta::Union{Float64, Void}(nothing)` : Discount factor. Randomly chosen from
[0, 1) if not specified.
- `;k::Union{Integer, Void}(nothing)` : Number of possible next states for each
state-action pair. Equal to `num_states` if not specified.

- `scale::Real(1)` : Standard deviation of the normal distribution for the
reward values.

##### Returns

- `ddp::DiscreteDP` : An instance of DiscreteDP.

"""
function random_discrete_dp(num_states::Integer,
                            num_actions::Integer,
                            beta::Union{Real, Void}=nothing;
                            k::Union{Integer, Void}=nothing,
                            scale::Real=1)
    L = num_states * num_actions
    R = scale * randn(L)
    Q = _random_stochastic_matrix(num_states, L; k=k)
    if beta == nothing
        beta = rand()
    end

    R = reshape(R, num_states, num_actions)
    Q = reshape(transpose(Q), num_states, num_actions, num_states)

    ddp = DiscreteDP(R, Q, beta)
    return ddp
end

# random_probvec

"""
Return m randomly sampled probability vectors of size k.

##### Arguments

- `k::Integer` : Size of each probability vector.
- `m::Integer` : Number of probability vectors.

##### Returns

- `a::Array` : Array of shape (k, m) containing probability vectors as colums.

"""
function random_probvec(k::Integer, m::Integer)
    k == 1 && return ones((k, m))

    # if k >= 2
    x = Array(Float64, (k, m))

    r = rand(k-1, m)
    x[1:end-1, :] = sort(r, 1)

    for j in 1:m
        x[end, j] = 1 - x[end-1, j]
        for i in k-1:-1:2
            x[i, j] -= x[i-1, j]
        end
    end

    return x
end


text:
using InventoryManagementWebApp

build_app(parse(Int,ARGS[1]))


text:
using DataFrames

"""
    encode_onehot(column::T[, categories::T[, prefix::String]]) where T <: AbstractArray

Converts a categorical column into a `DataFrame` of one-hot-encoded columns --
with one binary-encoded column per unique value in `column`.

# Examples

The basic version of this function makes a column for each unique
value in `column`.

```julia-repl
julia> data = [3,1,2,4];
julia> encode_onehot(data)
4×4 DataFrame
 Row │ 1      2      3      4
     │ Bool   Bool   Bool   Bool
─────┼────────────────────────────
   1 │ false  false   true  false
   2 │  true  false  false  false
   3 │ false   true  false  false
   4 │ false  false  false   true
```

You can also specify a prefix for each column.

```julia-repl
julia> data = [3,1,2,4];
julia> encode_onehot(data,"col_")
4×4 DataFrame
 Row │ col_1  col_2  col_3  col_4
     │ Bool   Bool   Bool   Bool
─────┼────────────────────────────
   1 │ false  false   true  false
   2 │  true  false  false  false
   3 │ false   true  false  false
   4 │ false  false  false   true
```

Additionally, you can specify the categories to convert to columns,
regardless of whether it exists in `column`.

```julia-repl
julia> data = [3,1,2,4];
julia> encode_onehot(data,[1:6;],"c")
4×6 DataFrame
 Row │ c1     c2     c3     c4     c5     c6
     │ Bool   Bool   Bool   Bool   Bool   Bool
─────┼──────────────────────────────────────────
   1 │ false  false   true  false  false  false
   2 │  true  false  false  false  false  false
   3 │ false   true  false  false  false  false
   4 │ false  false  false   true  false  false
```

See also: [`encode_dummy`](@ref), [`encode_hash`](@ref)
"""
function encode_onehot(column::T, categories::T, prefix::String) where T <: AbstractArray
    DataFrame(Dict([
        "$(prefix)$(l)" => column .== l
        for l = categories
    ]))
end

function encode_onehot(column::T, categories::T) where T <: AbstractArray
    encode_onehot(column,categories,"")
end

function encode_onehot(column::T, prefix::String) where T <: AbstractArray
    encode_onehot(column,levels(column),prefix)
end

function encode_onehot(column::T) where T <: AbstractArray
    encode_onehot(column,levels(column),"")
end

"""
    encode_dummy(column::T[, categories::T[, prefix::String]]) where T <: AbstractArray

Same as [`encode_onehot`](@ref) except that it drops the first column (to help prevent
issues caused by [multicollinearity](https://en.wikipedia.org/wiki/Multicollinearity)).

See also: [`encode_onehot`](@ref), [`encode_hash`](@ref)
"""
function encode_dummy(column::T, categories::T, prefix::String) where T <: AbstractArray
    df = encode_onehot(column,categories,prefix)
    if length(names(df)) > 1
        df[2:end]
    else df
    end
end

function encode_dummy(column::T, categories::T) where T <: AbstractArray
    df = encode_onehot(column,categories)
    if length(names(df)) > 1
        df[2:end]
    else df
    end
end

function encode_dummy(column::T, prefix::String) where T <: AbstractArray
    df = encode_onehot(column,prefix)
    if length(names(df)) > 1
        df[2:end]
    else df
    end
end

function encode_dummy(column::T) where T <: AbstractArray
    df = encode_onehot(column)
    if length(names(df)) > 1
        df[2:end]
    else df
    end
end

"""
    encode_hash(column::T, n_cols::Int = 8, prefix::String = "c") where T <: AbstractArray

Deterministically encode categorical features with high cardinality as a
`DataFrame` with `n_cols` columns.

```julia-repl
julia> data = [1:100;1000;];
julia> encode_hash(data)
10×8 DataFrame
 Row │ c1     c2     c3     c4     c5     c6     c7     c8
     │ Int64  Int64  Int64  Int64  Int64  Int64  Int64  Int64
─────┼────────────────────────────────────────────────────────
   1 │     1      1      1      1      1      1      1      0
   2 │     1      1      0      1      1      0      1      0
   3 │     1      0      1      0      1      1      0      1
   4 │     1      1      1      1      0      1      1      1
   5 │     1      0      0      0      0      0      1      0
   6 │     1      0      1      0      0      1      0      1
   7 │     0      1      1      0      0      1      1      1
   8 │     1      1      0      1      0      1      0      0
   9 │     0      1      0      1      0      1      1      0
  10 │     0      1      0      0      0      1      1      0
```

See also: [`encode_onehot`](@ref), [`encode_dummy`](@ref)
"""
function encode_hash(column::T, n_cols::Int = 8, prefix::String = "c") where T <: AbstractArray
    DataFrame(
        hcat(digits.(hash.(column).&(2^n_cols-1),base=2,pad=n_cols)...)',
        ["$(prefix)$(i)" for i = 1:n_cols]
    )
end


text:
__precompile__()

module GPU
	using OpenCL

	struct GPUOpts
		scr_width::Int32
		scr_height::Int32
		ignore_vertical_min::Int32
		ignore_vertical_max::Int32
		ignore_horizontal_min::Int32
		ignore_horizontal_max::Int32
		sector_width_top::Int32
		sector_width_right::Int32
		sector_width_bottom::Int32
		sector_width_left::Int32
		sector_offset_top::Int32
		sector_offset_right::Int32
		sector_offset_bottom::Int32
		sector_offset_left::Int32
		sector_count_bottom::Int32
		sector_count_left::Int32
		sector_total::Int32
	end

	struct GPUContext
		device::cl.Device
		ctx::cl.Context
		queue::cl.CmdQueue
		kernel::cl.Kernel
	end

	export GPUOpts, GPUContext

	function init_kernel()
		bin = Array{UInt8, 1}(undef, filesize("./bin/kernel.bin"))

		open("./bin/kernel.bin", "r") do f
			read!(f, bin)
		end

		device, ctx, queue = cl.create_compute_context()
		binaries = Dict(device => bin)
		p = cl.Program(ctx, binaries=binaries) |> cl.build!
		k = cl.Kernel(p, "avg")

		return GPUContext(device, ctx, queue, k)
	end

	function process_data(context::GPUContext, opts::GPUOpts, data::Array{UInt8, 1})
		zerosArr = zeros(Int32, opts.sector_total * 4)

		opts_buff = cl.Buffer(GPUOpts, context.ctx, (:r, :copy), hostbuf=[opts])
		raw_buff = cl.Buffer(UInt8, context.ctx, (:r, :copy), hostbuf=data)
		out_buff = cl.Buffer(Int32, context.ctx, (:rw, :copy), hostbuf=zerosArr)

		context.queue(context.kernel, opts.scr_height * opts.scr_width, nothing, opts_buff, raw_buff, out_buff)
		return cl.read(context.queue, out_buff)
	end

	export init_kernel, process_data
end


text:
g2 = evolving_graph(String, Int)
add_edge!(g2, "a", "b", 1)
add_edge!(g2, "b", "c", 2)
add_edge!(g2, "c", "d", 2)
add_edge!(g2, "a", "b", 2)
add_edge!(g2, "b", "a", 3)

k =  global_temporal_efficiency(g2, 1, 2)


text:
@testset "Literate examples" begin
    SKIPFILE = []

    EXAMPLEDIR = joinpath(@__DIR__, "../docs/src/examples")

    myfilter(str) = occursin(r"\.jl$", str) && !(str in SKIPFILE)
    for file in filter!(myfilter, readdir(EXAMPLEDIR))
        @testset "$file" begin
            mktempdir() do dir
                cd(dir) do
                    include(joinpath(EXAMPLEDIR, file))
                end
            end
        end
    end
end


text:
# Byte Set
# ========

struct ByteSet <: Base.AbstractSet{UInt8}
    a::UInt64  # 0x00:0x3F
    b::UInt64  # 0x40:0x7F
    c::UInt64  # 0x80:0xBF
    d::UInt64  # 0xC0:0xFF
end

function ByteSet()
    z = UInt64(0)
    return ByteSet(z, z, z, z)
end

function ByteSet(bytes::Union{UInt8,AbstractVector{UInt8},Set{UInt8}})
    a = b = c = d = UInt64(0)
    for byte in bytes
        if byte < 0x40
            a |= UInt64(1) << (byte - 0x00)
        elseif byte < 0x80
            b |= UInt64(1) << (byte - 0x40)
        elseif byte < 0xc0
            c |= UInt64(1) << (byte - 0x80)
        else
            d |= UInt64(1) << (byte - 0xc0)
        end
    end
    return ByteSet(a, b, c, d)
end

function Base.:(==)(s1::ByteSet, s2::ByteSet)
    return s1.a == s2.a && s1.b == s2.b && s1.c == s2.c && s1.d == s2.d
end

function Base.in(byte::UInt8, set::ByteSet)
    if byte < 0x40
        return set.a & (UInt64(1) << (byte - 0x00)) != 0
    elseif byte < 0x80
        return set.b & (UInt64(1) << (byte - 0x40)) != 0
    elseif byte < 0xc0
        return set.c & (UInt64(1) << (byte - 0x80)) != 0
    else
        return set.d & (UInt64(1) << (byte - 0xc0)) != 0
    end
end

Base.eltype(::Type{ByteSet}) = UInt8

function Base.show(io::IO, set::ByteSet)
    print(io, summary(set), "([", join(repr.(collect(set)), ','), "])")
end

function Base.length(set::ByteSet)
    return count_ones(set.a) + count_ones(set.b) + count_ones(set.c) + count_ones(set.d)
end

function Base.iterate(set::ByteSet, abcd=(set.a, set.b, set.c, set.d))
    a, b, c, d = abcd
    if a != 0
        byte = UInt8(trailing_zeros(a))
        a = xor(a, UInt64(1) << byte)
        byte += 0x00
    elseif b != 0
        byte = UInt8(trailing_zeros(b))
        b = xor(b, UInt64(1) << byte)
        byte += 0x40
    elseif c != 0
        byte = UInt8(trailing_zeros(c))
        c = xor(c, UInt64(1) << byte)
        byte += 0x80
    elseif d != 0
        byte = UInt8(trailing_zeros(d))
        d = xor(d, UInt64(1) << byte)
        byte += 0xc0
    else
        return nothing
    end
    return byte, (a, b, c, d)
end

function Base.union(s1::ByteSet, s2::ByteSet)
    return ByteSet(s1.a | s2.a, s1.b | s2.b, s1.c | s2.c, s1.d | s2.d)
end

function Base.intersect(s1::ByteSet, s2::ByteSet)
    return ByteSet(s1.a & s2.a, s1.b & s2.b, s1.c & s2.c, s1.d & s2.d)
end

function Base.setdiff(s1::ByteSet, s2::ByteSet)
    return ByteSet(s1.a & ~s2.a, s1.b & ~s2.b, s1.c & ~s2.c, s1.d & ~s2.d)
end

function Base.minimum(set::ByteSet)
    if set.a != 0x00
        return UInt8(trailing_zeros(set.a))
    elseif set.b != 0x00
        return UInt8(trailing_zeros(set.b)) + 0x40
    elseif set.c != 0x00
        return UInt8(trailing_zeros(set.c)) + 0x80
    elseif set.d != 0x00
        return UInt8(trailing_zeros(set.d)) + 0xc0
    else
        throw(ArgumentError("empty set"))
    end
end

function Base.maximum(set::ByteSet)
    if set.d != 0x00
        return UInt8(63 - leading_zeros(set.d)) + 0xc0
    elseif set.c != 0x00
        return UInt8(63 - leading_zeros(set.c)) + 0x80
    elseif set.b != 0x00
        return UInt8(63 - leading_zeros(set.b)) + 0x40
    elseif set.a != 0x00
        return UInt8(63 - leading_zeros(set.a))
    else
        throw(ArgumentError("empty set"))
    end
end

function isdisjoint(s1::ByteSet, s2::ByteSet)
    return isempty(intersect(s1, s2))
end

# Encode a byte set into a non-empty sequence of ranges.
function range_encode(set::ByteSet)
    labels = collect(set)
    labels′ = UnitRange{UInt8}[]
    while !isempty(labels)
        lo = popfirst!(labels)
        hi = lo
        while !isempty(labels) && first(labels) == hi + 1
            hi = popfirst!(labels)
        end
        push!(labels′, lo:hi)
    end
    return labels′
end


text:
# # Collinear spin and magnetic systems
#
# In this example we consider iron in the BCC phase.
# To show that this material is ferromagnetic we will model it once
# allowing collinear spin polarization and once without
# and compare the resulting SCF energies. In particular
# the ground state can only be found if collinear spins are allowed.
#
# First we setup BCC iron without spin polarization
# using a single iron atom inside the unit cell.
using DFTK

a = 5.42352  # Bohr
lattice = a / 2 * [[-1  1  1];
                   [ 1 -1  1];
                   [ 1  1 -1]]
Fe = ElementPsp(:Fe, psp=load_psp("hgh/lda/Fe-q8.hgh"))
atoms = [Fe => [zeros(3)]];

# To get the ground-state energy we use an LDA model and rather moderate
# discretisation parameters.

kgrid = [3, 3, 3]  # k-point grid (Regular Monkhorst-Pack grid)
Ecut = 15          # kinetic energy cutoff in Hartree
model_nospin = model_LDA(lattice, atoms, temperature=0.01)
basis_nospin = PlaneWaveBasis(model_nospin; kgrid, Ecut)

scfres_nospin = self_consistent_field(basis_nospin, tol=1e-6, mixing=KerkerMixing());
#-
scfres_nospin.energies

# Since we did not specify any initial magnetic moment on the iron atom,
# DFTK will automatically assume that a calculation with only spin-paired
# electrons should be performed. As a result the obtained ground state
# features no spin-polarization.

# Now we repeat the calculation, but give the iron atom an initial magnetic moment.
# For specifying the magnetic moment pass the desired excess of spin-up over spin-down
# electrons at each centre to the `Model` and the guess density functions.
# In this case we seek the state with as many spin-parallel
# ``d``-electrons as possible. In our pseudopotential model the 8 valence
# electrons are 2 pair of ``s``-electrons, 1 pair of ``d``-electrons
# and 4 unpaired ``d``-electrons giving a desired magnetic moment of `4` at the iron centre.
# The structure (i.e. pair mapping and order) of the `magnetic_moments` array needs to agree
# with the `atoms` array and `0` magnetic moments need to be specified as well.

magnetic_moments = [Fe => [4, ]];

# !!! tip "Units of the magnetisation and magnetic moments in DFTK"
#     Unlike all other quantities magnetisation and magnetic moments in DFTK
#     are given in units of the Bohr magneton ``μ_B``, which in atomic units has the
#     value ``\frac{1}{2}``. Since ``μ_B`` is (roughly) the magnetic moment of
#     a single electron the advantage is that one can directly think of these
#     quantities as the excess of spin-up electrons or spin-up electron density.
#
# We repeat the calculation using the same model as before. DFTK now detects
# the non-zero moment and switches to a collinear calculation.

model = model_LDA(lattice, atoms, magnetic_moments=magnetic_moments, temperature=0.01)
basis = PlaneWaveBasis(model; Ecut, kgrid)
ρ0 = guess_density(basis, magnetic_moments)
scfres = self_consistent_field(basis, tol=1e-6; ρ=ρ0, mixing=KerkerMixing());
#-
scfres.energies

# !!! note "Model and magnetic moments"
#     DFTK does not store the `magnetic_moments` inside the `Model`, but only uses them
#     to determine the lattice symmetries. This step was taken to keep `Model`
#     (which contains the physical model) independent of the details of the numerical details
#     such as the initial guess for the spin density.
#
# In direct comparison we notice the first, spin-paired calculation to be
# a little higher in energy
println("No magnetization: ", scfres_nospin.energies.total)
println("Magnetic case:    ", scfres.energies.total)
println("Difference:       ", scfres.energies.total - scfres_nospin.energies.total);
# Notice that with the small cutoffs we use to generate the online
# documentation the calculation is far from converged.
# With more realistic parameters a larger energy difference of about
# 0.1 Hartree is obtained.

# The spin polarization in the magnetic case is visible if we
# consider the occupation of the spin-up and spin-down Kohn-Sham orbitals.
# Especially for the ``d``-orbitals these differ rather drastically.
# For example for the first ``k``-point:

iup   = 1
idown = iup + length(scfres.basis.kpoints) ÷ 2
@show scfres.occupation[iup][1:7]
@show scfres.occupation[idown][1:7];

# Similarly the eigenvalues differ
@show scfres.eigenvalues[iup][1:7]
@show scfres.eigenvalues[idown][1:7];

# !!! note "k-points in collinear calculations"
#     For collinear calculations the `kpoints` field of the `PlaneWaveBasis` object contains
#     each ``k``-point coordinate twice, once associated with spin-up and once with down-down.
#     The list first contains all spin-up ``k``-points and then all spin-down ``k``-points,
#     such that `iup` and `idown` index the same ``k``-point, but differing spins.

# We can observe the spin-polarization by looking at the density of states (DOS)
# around the Fermi level, where the spin-up and spin-down DOS differ.

using Plots
plot_dos(scfres)

# Similarly the band structure shows clear differences between both spin components.
using Unitful
using UnitfulAtomic
plot_bandstructure(scfres; kline_density=6)


text:
"""
Convenience code for plotting spectra, associated data types and results of analysis

Author: Eric Ford
Created: August 2020
"""

""" plot_spectrum_chunks( chunklist, chunks )
Returns a plot based of flux versus wavelength for selected chunks from the spectra in chunklist.
"""
function plot_spectrum_chunks(clt::CLT, chunks::Union{Integer,AbstractUnitRange};
    time_idx::Union{Integer,AbstractUnitRange}=1, color = nothing,
    plt::PT = Plots.plot(legend=:none)) where {CLT<:AbstractChunkListTimeseries, PT<:AbstractPlot}
    xmin = Inf # minimum(order_list_df.lambda_lo[chunk_idx])
    xmax = 0   # maximum(order_list_df.lambda_hi[chunk_idx])
    for c in chunks
        t = 1
        if color == nothing
            plot!(plt,clt.chunk_list[t].data[c].λ ,clt.chunk_list[t].data[c].flux)
        else
            colorval = typeof(color)<:Union{Integer,Symbol} ? color : color[c]
            plot!(plt,clt.chunk_list[t].data[c].λ ,clt.chunk_list[t].data[c].flux, linecolor=colorval)
        end
        xmin = min(xmin,minimum(clt.chunk_list[t].data[c].λ))
        xmax = max(xmax,maximum(clt.chunk_list[t].data[c].λ))
    end
    xlabel!(plt,"λ (Å)")
    ylabel!(plt,"Flux)")
    #xlims!(xmin,xmax)
    plt
end

"""" get_λs( grids::Vector{AbstractRange}, r::AbstractUnitRange )
Return range of wavelengths corresponding to indices in r from a vector of ranges.
Useful for extracting wavelengths from a chunklist.
"""
function get_λs(grids::V1, r::AR) where { T1<:AbstractRange, V1<:AbstractVector{T1}, AR<:AbstractUnitRange }
  local idx_offset = 1
  #println("first(r) = ",first(r), "   last(r) = ",last(r) )
  for chunk_range in grids
    #println("c: ",c," chunk_range=", chunk_range, " idx_offset=",idx_offset," - ", idx_offset+length(chunk_range)-1)
    if (first(r)>=idx_offset) && (last(r)<=idx_offset+length(chunk_range)-1)
       idx = (first(r)-idx_offset+1):(last(r)-idx_offset+1)
       return chunk_range[idx]
    else
      idx_offset += length(chunk_range)
    end
  end
  @error "get_λs: Failed to find range (",r,")."
  return
end

""" `add_time_gap_lines(plt, times )`
Inputs:
- plt:  AbstractPlot to add to
- times: observation times
Optional Inputs:
- Δt_boundaries: list of interval boundaries for different line types
- lweights: weights of lines dividing different intervals
- lstyles: styles of lines  dividing different intervals
"""
function add_time_gap_lines(plt::AbstractPlot, times::AbstractVector{T} ; Δt_boundaries::AbstractVector{T2} = [0.5, 1.5, 7, 14, 42],
                                    lweights::AbstractVector{T3} = [ 1, 1, 1, 1, 2, 2, 3, 3  ],
                                    lstyles::AbstractVector{Symbol} = [ :dot, :dashdot, :dash, :solid, :dash, :solid, :dash, :solid ]
                                          ) where { T<:Real, T2<:Real, T3<:Real  }
   @assert 2 <= length(Δt_boundaries) <= length(lstyles)-1
   Δt = times[2:end] .- times[1:end-1]
   xvals = [xlims(plt)[1], xlims(plt)[2]]
   for i in 1:(length(Δt_boundaries)-1)
      idx = findall(x -> Δt_boundaries[i]<=x<Δt_boundaries[i+1], Δt)
      if length(idx) >= 1
         map( y->plot!(plt,xvals, (y+0.5).*[1, 1],label=:none, color=:black, linewidth=lweights[i], linestyle=lstyles[i]), idx )
      end
   end
   idx = findall(x -> x > last(Δt_boundaries), Δt)
   if length(idx) >= 1
      map(y->plot!(plt,xvals, (y+0.5).*[1, 1],label=:none, color=:black, linewidth=lweights[length(idx)], linestyle=lstyles[length(idx)]), idx )
   end
   return plt
end


text:
"""
Create a random partition for a matrix `A`. Each node is assigned to a module at
random. If the smallest dimension of the matrix is *m*, there will be between 1
and *m* modules.

This method tends to work well for small to medium networks (and is faster than
`partition_lp`).
"""
function partition_random(A::Array{Int64, 2})
  # Smaller margin of A determines the number of modules
  min_part = minimum(size(A))
  # Generates a matrix with int(0) for the module partition
  S = zeros(Int64, sum(size(A)), min_part)
  # Pick modules at random and put them in the partition
  init_part = rand(1:1:min_part,sum(size(A)))
  for i = 1:length(init_part)
    S[i,init_part[i]] = 1
  end
  # Only non-empty modules are allowed
  # And it returns a Modular object
  return no_empty_modules!(Modular(A, S))
end


text:
# eval hamiltonian
function hamiltoniank!(mat::AbstractMatrix, sys::BlochHamiltonian, k::AbstractVector)
    @assert length(k) == D "length of k vector much be the same as the dimensionality of the system"

    sys.hk!(mat, SVector{D}(k))

end

function hamiltoniank(sys::BlochHamiltonian, k::AbstractVector)
    @assert length(k) == D "length of k vector much be the same as the dimensionality of the system"

    dim = norbs()
    mat = zeros(Complex{Float64}, )

    sys.hk!(mat, SVector{D}(k))

end


text:
export AbstractBlock, PrimitiveBlock, CompositeBlock, AbstractContainer, TagBlock

"""
    AbstractBlock

Abstract type for quantum circuit blocks.
"""
abstract type AbstractBlock{N} end

"""
    PrimitiveBlock{N} <: AbstractBlock{N}

Abstract type that all primitive block will subtype from. A primitive block
is a concrete block who can not be decomposed into other blocks. All composite
block can be decomposed into several primitive blocks.

!!! note

    subtype for primitive block with parameter should implement `hash` and `==`
    method to enable key value cache.
"""
abstract type PrimitiveBlock{N} <: AbstractBlock{N} end

"""
    CompositeBlock{N} <: AbstractBlock{N}

Abstract supertype which composite blocks will inherit from. Composite blocks
are blocks composited from other [`AbstractBlock`](@ref)s, thus it is a `AbstractBlock`
as well.
"""
abstract type CompositeBlock{N} <: AbstractBlock{N} end

"""
    AbstractContainer{BT, N} <: CompositeBlock{N}

Abstract type for container block. Container blocks are blocks contain a single
block. Container block should have a
"""
abstract type AbstractContainer{BT<:AbstractBlock,N} <: CompositeBlock{N} end

"""
    TagBlock{BT, N} <: AbstractContainer{BT, N}

`TagBlock` is a special kind of Container block, it forwards most of the methods
but tag the block with some extra information.
"""
abstract type TagBlock{BT,N} <: AbstractContainer{BT,N} end

"""
    content(x)

Returns the content of `x`.
"""
@interface content

"""
    chcontent(x, blk)

Create a similar block of `x` and change its content to blk.
"""
@interface chcontent

"""
    apply!(register, block)

Apply a block (of quantum circuit) to a quantum register.
"""

@interface apply!

"""
    occupied_locs(x)

Return a tuple of occupied locations of `x`.
"""
@interface occupied_locs

"""
    subblocks(x)

Returns an iterator of the sub-blocks of a composite block. Default is empty.
"""
@interface subblocks

"""
    chsubblocks(composite_block, itr)

Change the sub-blocks of a [`CompositeBlock`](@ref) with given iterator `itr`.
"""
@interface chsubblocks

@interface print_block

"""
    mat([T=ComplexF64], blk)

Returns the matrix form of given block.
"""
@interface mat

# parameters
"""
    getiparams(block)

Returns the intrinsic parameters of node `block`, default is an empty tuple.
"""
@interface getiparams

"""
    setiparams!(block, itr)
    setiparams!(block, params...)

Set the parameters of `block`.
"""
@interface setiparams!

"""
    parameters(block)

Returns all the parameters contained in block tree with given root `block`.
"""
@interface parameters

"""
    nparameters(block) -> Int

Return number of parameters in `block`. See also [`niparams`](@ref).
"""
@interface nparameters

"""
    niparam(block) -> Int

Return number of intrinsic parameters in `block`. See also [`nparameters`](@ref).
"""
@interface niparams

"""
    iparams_eltype(block)

Return the element type of [`getiparams`](@ref).
"""
@interface iparams_eltype

"""
    parameters_eltype(x)

Return the element type of [`parameters`](@ref).
"""
@interface parameters_eltype

"""
    dispatch!(x::AbstractBlock, collection)

Dispatch parameters in collection to block tree `x`.

!!! note

    it will try to dispatch the parameters in collection first.
"""
@interface dispatch!

"""
    render_params(r::AbstractBlock, params)

`params` can be a number or a symbol like `:zero` and `:random`.
This function renders the input parameter to a consumable type to `r`.
"""
@interface render_params

"""
    expect(op::AbstractBlock, reg) -> Vector
    expect(op::AbstractBlock, reg => circuit) -> Vector
    expect(op::AbstractBlock, density_matrix) -> Vector

Get the expectation value of an operator, the second parameter can be a register `reg` or a pair of input register and circuit `reg => circuit`.

expect'(op::AbstractBlock, reg=>circuit) -> Pair
expect'(op::AbstractBlock, reg) -> AbstracRegister

Obtain the gradient with respect to registers and circuit parameters.
For pair input, the second return value is a pair of `gψ=>gparams`,
with `gψ` the gradient of input state and `gparams` the gradients of circuit parameters.
For register input, the return value is a register.

!!! note

    For batched register, `expect(op, reg=>circuit)` returns a vector of size number of batch as output. However, one can not differentiate over a vector loss, so `expect'(op, reg=>circuit)` accumulates the gradient over batch, rather than returning a batched gradient of parameters.
"""
@interface expect

"""
    operator_fidelity(b1::AbstractBlock, b2::AbstractBlock) -> Number

Operator fidelity defined as

```math
F^2 = \\frac{1}{d^2}\\left[{\\rm Tr}(b1^\\dagger b2)\\right]
```

Here, `d` is the size of the Hilbert space. Note this quantity is independant to global phase.
See arXiv: 0803.2940v2, Equation (2) for reference.
"""
@interface operator_fidelity

"""
    apply_back!((ψ, ∂L/∂ψ*), circuit::AbstractBlock, collector) -> AbstractRegister

back propagate and calculate the gradient ∂L/∂θ = 2*Re(∂L/∂ψ*⋅∂ψ*/∂θ), given ∂L/∂ψ*.
`ψ` is the output register, ∂L/∂ψ* should also be register type.

Note: gradients are stored in `Diff` blocks, it can be access by either `diffblock.grad` or `gradient(circuit)`.
Note2: now `apply_back!` returns the inversed gradient!
"""
@interface apply_back!

"""
    mat_back!(T, rb::AbstractBlock, adjy, collector)

Back propagate the matrix gradients.
"""
@interface mat_back!


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

module Resolve

include("resolve/versionweight.jl")
include("resolve/interface.jl")
include("resolve/maxsum.jl")

using ..Types, ..Query, .Pkg25ToMaxSumInterface, .MaxSum
import ...Pkg25.Pkg25Error

export resolve, sanity_check

# Use the max-sum algorithm to resolve packages dependencies
function resolve(reqs::Requires, deps::Dict{String,Dict{VersionNumber,Available}})
    # init interface structures
    interface = Interface(reqs, deps)

    # attempt trivial solution first
    ok, sol = greedysolver(interface)
    if !ok
        # trivial solution failed, use maxsum solver
        graph = Graph(interface)
        msgs = Messages(interface, graph)

        try
            sol = maxsum(graph, msgs)
        catch err
            if isa(err, UnsatError)
                p = interface.pkgs[err.info]
                msg = "unsatisfiable package requirements detected: " *
                      "no feasible version could be found for package: $p"
                if msgs.num_nondecimated != graph.np
                    msg *= "\n  (you may try increasing the value of the" *
                           "\n   JULIA_PKGRESOLVE_ACCURACY environment variable)"
                end
                throw(Pkg25Error(msg))
            end
            rethrow(err)
        end

        # verify solution (debug code) and enforce its optimality
        @assert verify_solution(sol, interface)
        enforce_optimality!(sol, interface)
    end

    # return the solution as a Dict mapping package_name => sha1
    return compute_output_dict(sol, interface)
end

# Scan dependencies for (explicit or implicit) contradictions
function sanity_check(deps::Dict{String,Dict{VersionNumber,Available}},
                      pkgs::Set{String} = Set{String}())
    isempty(pkgs) || (deps = Query.undirected_dependencies_subset(deps, pkgs))

    deps, eq_classes = Query.prune_versions(deps)

    ndeps = Dict{String,Dict{VersionNumber,Int}}()

    for (p,depsp) in deps
        ndeps[p] = ndepsp = Dict{VersionNumber,Int}()
        for (vn,a) in depsp
            ndepsp[vn] = length(a.requires)
        end
    end

    vers = Array{Tuple{String,VersionNumber,VersionNumber}}(0)
    for (p,d) in deps, vn in keys(d)
        lvns = VersionNumber[filter(vn2->(vn2>vn), keys(d))...]
        nvn = isempty(lvns) ? typemax(VersionNumber) : minimum(lvns)
        push!(vers, (p,vn,nvn))
    end
    sort!(vers, by=pvn->(-ndeps[pvn[1]][pvn[2]]))

    nv = length(vers)

    svdict = Dict{Tuple{String,VersionNumber},Int}(vers[i][1:2]=>i for i = 1:nv)

    checked = falses(nv)

    problematic = Array{Tuple{String,VersionNumber,String}}(0)
    i = 1
    psl = 0
    for (p,vn,nvn) in vers
        if ndeps[p][vn] == 0
            break
        end
        if checked[i]
            i += 1
            continue
        end

        sub_reqs = Dict{String,VersionSet}(p=>VersionSet([vn, nvn]))
        sub_deps = Query.prune_dependencies(sub_reqs, deps)
        interface = Interface(sub_reqs, sub_deps)

        red_pkgs = interface.pkgs
        red_np = interface.np
        red_spp = interface.spp
        red_pvers = interface.pvers

        ok, sol = greedysolver(interface)

        if !ok
            try
                graph = Graph(interface)
                msgs = Messages(interface, graph)
                sol = maxsum(graph, msgs)
                ok = verify_solution(sol, interface)
                @assert ok
            catch err
                isa(err, UnsatError) || rethrow(err)
                pp = red_pkgs[err.info]
                for vneq in eq_classes[p][vn]
                    push!(problematic, (p, vneq, pp))
                end
            end
        end
        if ok
            let
                p0 = interface.pdict[p]
                svn = red_pvers[p0][sol[p0]]
                @assert svn == vn
            end

            for p0 = 1:red_np
                s0 = sol[p0]
                if s0 != red_spp[p0]
                    j = svdict[(red_pkgs[p0], red_pvers[p0][s0])]
                    checked[j] = true
                end
            end
            checked[i] = true
        end
        i += 1
    end

    return sort!(problematic)
end

end # module


text:
using Decimals
using Base.Test

# Number/string to decimal conversions
@test decimal("0.01") == decimal(0.01) == Decimal(0, 1, -2)
@test decimal(".001") == decimal(.001) == Decimal(0, 1, -3)
@test decimal("15.23") == decimal(15.23) == Decimal(0, 1523, -2)
@test decimal("543") == decimal(543) == Decimal(0, 543, 0)
@test decimal("-345") == decimal(-345) == Decimal(1, 345, 0)
@test decimal("000123") == decimal(000123) == Decimal(0, 123, 0)
@test decimal("-00032") == decimal(-00032) == Decimal(1, 32, 0)
@test decimal("200100") == decimal(200100) == Decimal(0, 2001, 2)
@test decimal("-.123") == decimal(-.123) == Decimal(1, 123, -3)
@test decimal("1.23000") == decimal(1.23000) == Decimal(0, 123, -2)
@test decimal("4734.612") == decimal(4734.612) == Decimal(0, 4734612, -3)
@test decimal("541724.2") == decimal(541724.2) == Decimal(0,5417242,-1)
@test decimal("2.5e6") == decimal(2.5e6) == Decimal(0, 25, 5)
@test decimal("2.385350e8") == decimal(2.385350e8) == Decimal(0, 238535, 3)
@test decimal("12.3e-4") == decimal(12.3e-4) == Decimal(0, 123, -5)

@test decimal("-12.3e4") == decimal(-12.3e4) == Decimal(1, 123, 3)

@test decimal("-12.3e-4") == decimal(-12.3e-4) == Decimal(1, 123, -5)

@test decimal("0.1234567891") == decimal(0.1234567891) == Decimal(0,1234567891, -10)
@test decimal("0.12345678912") == decimal(0.12345678912) == Decimal(0,12345678912, -11)

# Array-to-decimal conversion
@test decimal([0.1 0.2 0.3]) == [decimal(0.1) decimal(0.2) decimal(0.3)]

# Decimal-to-string conversions
@test string(Decimal(0, 1, -2)) == "0.01"
@test string(Decimal(0, 1, -3)) == "0.001"
@test string(Decimal(0, 1523, -2)) == "15.23"
@test string(Decimal(0, 543, 0)) == "543"
@test string(Decimal(1, 345, 0)) == "-345"
@test string(Decimal(0, 123, 0)) == "123"
@test string(Decimal(1, 32, 0)) == "-32"
@test string(Decimal(0, 2001, 2)) == "200100"
@test string(Decimal(1, 123, -3)) == "-0.123"
@test string(Decimal(0, 123, -2)) == "1.23"

# Decimal-to-number conversions
@test float(Decimal(0, 1, -2)) == 0.01
@test float(Decimal(0, 1, -3)) == 0.001
@test float(Decimal(0, 1523, -2)) == 15.23
@test float(Decimal(0, 543, 0)) == 543
@test float(Decimal(1, 345, 0)) == -345
@test float(Decimal(0, 123, 0)) == 123
@test float(Decimal(1, 32, 0)) == -32
@test float(Decimal(0, 2001, 2)) == 200100
@test float(Decimal(1, 123, -3)) == -0.123
@test float(Decimal(0, 123, -2)) == 1.23
@test number(Decimal(0, 1, -2)) == 0.01
@test number(Decimal(0, 1, -3)) == 0.001
@test number(Decimal(0, 1523, -2)) == 15.23
@test number(Decimal(0, 543, 0)) == 543
@test number(Decimal(1, 345, 0)) == -345
@test number(Decimal(0, 123, 0)) == 123
@test number(Decimal(1, 32, 0)) == -32
@test number(Decimal(0, 2001, 2)) == 200100
@test number(Decimal(1, 123, -3)) == -0.123
@test number(Decimal(0, 123, -2)) == 1.23
@test string(float(Decimal(0, 543, 0))) == "543.0"
@test string(number(Decimal(0, 543, 0))) == "543"
@test string(number(Decimal(0, 543, -1))) == "54.3"


text:
module WebDataset

using TarIterators
using Images
using FileIO, ImageIO, ImageMagick
using Transducers
using Transducers: Transducer, R_, next, inner, xform
using HTTP
using JSON
using ResumableFunctions
using Base.Iterators
using Test
using Base.Threads
using Parameters

export Continue
export DatasetDescriptor
export Rename
export allsplitext
export braceexpand
export clsdecode
export dataloader
export debugloader
export default_decoders
export generic_open
export imdecode
export jsondecode
export loadproc
export map_by_rules
export samplebatching
export sampleiterator
export sampleshuffle
export sampletransforms
export tariterator
export transpose_batch

global dataloader_tids = nothing
global last_exception = nothing
global num_exceptions = Threads.Atomic{Int}(0)
global max_exceptions = 999999999

# path names

function allsplitext(s)
    if s == ""; return ("", ""); end
    slash = findlast("/", s)
    slash = slash == nothing ? (0,) : slash
    fname = s[slash[1]+1:length(s)]
    dot = findfirst(".", fname)
    dot = dot == nothing ? (0,) : dot
    split = slash[1]+dot[1]
    return s[1:split-1], s[split:length(s)]
end

# generic I/O

function generic_open(shard; verbose=false)
    if startswith(shard, "pipe:")
        shard = Cmd(Base.shell_split(shard[6:length(shard)]))
    end
    verbose && @show shard
    return open(shard)
end

# brace expansion for filenames

function braceexpand(s::String)
    expansion = r"[{][0-9]+[.][.][0-9]+[}]"
    m = match(expansion, s)
    if m == nothing
        return [s]
    end
    prefix = s[1:m.offset-1]
    rest = braceexpand(s[m.offset+length(m.match):length(s)])
    range = s[m.offset+1:m.offset+length(m.match)-2]
    lohi = split(range, "..")
    (lo, hi) = [parse(Int, x) for x in lohi]
    result = []
    for i in lo:hi
        for r in rest
            expanded = prefix * string(i, pad=length(lohi[1])) * r
            push!(result, expanded)
        end
    end
    return result
end

function braceexpand(l::Vector{String})
    expanded = braceexpand(l[1])
    if length(l) < 2
        return expanded
    end

end

# .tar file reader

function substr(header, from, size)
    lo = from+1
    hi = from+size
    s = String(header[lo:hi])
    if length(s) == 0; return s; end
    return rstrip(s, ['\0'])
end

function decode_header(header)
    htype = substr(header, 156, 1)
    fname = substr(header, 0, 100)
    prefix = substr(header, 345, 155)
    raw_size = substr(header, 124, 12)
    size = 0
    @test !occursin("\0", raw_size)
    if raw_size != ""; size = parse(Int64, raw_size, base=8); end
    blocks = Int(floor((size + 511) / 512))
    return (htype, prefix*fname, size, blocks)
end

function next_file(stream)
    header = Array{UInt8}(undef, 512)
    read!(stream, header)
    htype, fname, n, blocks = decode_header(header)
    buffer = Array{UInt8}(undef, blocks*512)
    read!(stream, buffer)
    return (htype, fname, buffer[1:n])
end

@resumable function tariterator(stream)
    count = 0
    while(!eof(stream))
        htype, fname, buffer = next_file(stream)
        if htype != "0"; continue; end
        @yield (fname, buffer)
        count = count + 1
    end
end

# grouping tar entries into samples

valid_sample(d) = (d != nothing && length(d) > 1)

@resumable function sampleiterator(stream)
    current_prefix = nothing
    current_sample = nothing
    for (fname, value) in tariterator(stream)
        (prefix, suffix) = allsplitext(fname)
        if prefix == ""; continue; end
        if current_sample === nothing || prefix != current_prefix
            if valid_sample(current_sample)
                @yield current_sample
            end
            current_sample = Dict{String,Any}("__key__"=>prefix)
            current_prefix = prefix
        end
        current_sample[suffix] = value
    end
    if valid_sample(current_sample)
        @yield current_sample
    end
end

# shuffling

@resumable function sampleshuffle(source, bufsize=1000)
    buffer = []
    iter = iterate(source)
    while iter != nothing
        if length(buffer) < bufsize
            push!(buffer, iter[1])
            iter = iterate(source, iter[2])
        end
        iter == nothing && break
        index = rand(UInt64) % length(buffer) + 1
        result = buffer[index]
        buffer[index] = iter[1]
        @yield result
        iter = iterate(source, iter[2])
    end
    while length(buffer) > 0
        @yield pop!(buffer)
    end
end

# decoding and augmentation

struct EOF
    key
end

abstract type Renamer end

struct Rename <: Renamer
    key
    new_key
end

struct Continue <: Renamer
    key
    new_key
end

function kv_apply(f, key, value)
    try
        return f(value)
    catch e
        if e isa MethodError
            return f(key, value)
        end
        rethrow()
    end
end

function maprule(key::String, value, rule::String, f)
    !endswith(lowercase(key), rule) && return nothing
    return (key, f(value))
end

function maprule(key::String, value, rule::Regex, f)
    match(rule, lowercase(key)) == nothing && return nothing
    return (key, f(value))
end

function maprule(key::String, value, rule::Renamer, f)
    out = maprule(key, value, rule.key, f)
    out == nothing && return nothing
    return (rule.new_key, out[2])
end

function maprule(key::String, value, rule::Type, f)
    !isa(value, rule) && return nothing
    return (key, f(value))
end

function map_by_rules(sample, transformers)
    result = Dict{String,Any}()
    for (key, value) in sample
        if key[1] == '_'
            result[key] = value
            continue
        end
        for (rule, f) in transformers
            output = maprule(key, value, rule, f)
            output == nothing && continue
            new_key, new_value = output
            result[new_key] = new_value
            !isa(rule, Continue) && break
        end
    end
    return result
end

@resumable function sampletransforms(source, transformers; handle="ignore")
    for sample in source
        result = missing
        try
            result = map_by_rules(sample, transformers)
        catch exn
            @show exn
            global last_exception
            global num_exceptions
            last_exception = exn
            Threads.atomic_add!(num_exceptions, 1)
            if num_exceptions[] > max_exceptions
                error("too many exceptions")
            end
            continue
        end
        @yield result
    end
end

function imdecode(data)
    return ImageMagick.load_(data)
end

function clsdecode(data; lo=-999999, hi=999999)
    value = parse(Int64, String(data))
    if value < lo || value > hi
        error("cls value out of range")
    end
    return value
end

function jsondecode(data; query="")
    query == "" || error("unimplemented")
    return JSON.parse(String(data))
end

default_decoders = [
    ".cls" => clsdecode,
    Rename(r".(jpg|jpg|png|p?m)$", ".img") => imdecode,
    ".json" => jsondecode,
]

# batching

function transpose_batch(sample)
    keys = Set{String}(key for d in sample for (key, _) in d)
    result = Dict{String,Any}()
    for key in keys
        result[key] = [d[key] for d in sample]
    end
    if length(Set(length(l) for (_, l) in result)) != 1
        throw(error("not all keys are present in all samples"))
    end
    return result
end

@resumable function samplebatching(source, batchsize; minsize=1)
    current_batch = []
    for sample in source
        if batchsize == 0
            @yield sample
        elseif length(current_batch) >= batchsize
            @yield transpose_batch(current_batch)
            current_batch = []
        else
            push!(current_batch, sample)
        end
    end
    if length(current_batch) >= minsize
        @yield transpose_batch(current_batch)
    end
end

# multithreaded loader

@with_kw mutable struct DatasetDescriptor
    # list of shards
    sources::Array{Union{String,Cmd}} = []
    # size of the inline shuffle buffer
    shuffle::Int = 1000
    # batchsize (0 for no batching)
    batchsize::Int = 16
    # decoding functions
    decoding::Array{Pair{Union{String,Regex,Renamer},Function}} = [""=>x->x]
    # augmentation functions
    augmenting::Array{Pair{Union{String,Regex,Renamer},Function}} = [""=>x->x]
    # a general mapping function running after augmentation
    mapping::Function = x->x
    # collation functions running after batching
    collating::Array{Pair{Union{String,Regex,Renamer},Function}} = [""=>x->x]
    # channel size for output channel
    csize::Int = 100
    # number of concurrent tasks/shards for loading
    ntasks::Int = 4
    # more verbose output
    verbose::Bool = false
    # even more verbose output
    debug::Bool = false
end

@resumable function ldebug(it)
    for sample in it
        @show sample["__key__"], string(sample)[1:100]
        @yield sample
    end
end

function loadproc(desc::DatasetDescriptor, inch::Channel, outch::Channel, eof; maxcount=1e30)
    count = 0
    for source in inch
        desc.debug && @show source
        stream = generic_open(source)
        raw = sampleiterator(stream)
        shuffled = sampleshuffle(raw, desc.shuffle)
        decoded = sampletransforms(shuffled, desc.decoding)
        augmented = sampletransforms(decoded, desc.augmenting)
        mapped = (desc.mapping(x) for x in augmented)
        batched = samplebatching(mapped, desc.batchsize)
        collated = sampletransforms(batched, desc.collating)
        for batch in collated
            desc.debug && @show batch["__key__"]
            put!(outch, batch)
            sleep(0.0001)
            count += 1
            if count > maxcount; return; end
        end
        close(stream)
    end
    put!(outch, eof)
end

function checkfailed(tids)
    for t in tids
        if istaskfailed(t)
            wait(t)
        end
    end
end

@resumable function dataloader(desc::DatasetDescriptor, inch::Channel)
    global num_exceptions
    global dataloader_tids
    Threads.atomic_xchg!(num_exceptions, 0)
    outch = Channel(desc.csize)
    tids = [Threads.@spawn loadproc(desc, inch, outch, EOF(i)) for i in 1:desc.ntasks]
    dataloader_tids = tids
    running = Set(1:desc.ntasks)
    while true
        checkfailed(tids)
        tids = [t for t in tids if !istaskdone(t)]
        sample = take!(outch)
        if isa(sample, EOF)
            delete!(running, sample.key)
            length(running) == 0 && break
            continue
        end
        @yield sample
    end
end

function dataloader(desc::DatasetDescriptor)
    fnames = desc.sources
    sources = Channel(length(fnames))
    for fname in fnames
        put!(sources, fname)
    end
    close(sources)
    return dataloader(desc, sources)
end

function debugloader(desc::DatasetDescriptor, sources; maxcount=10)
    inch = Channel(length(sources))
    for s in sources; put!(inch, s); end
    close(inch)
    outch = Channel(min(maxcount+1, 100))
    loadproc(desc, inch, outch, EOF(0); maxcount=maxcount)
    return outch
end

end


text:
using MittagLefflerFunctions
using LinearAlgebra: cond, svd, lu, norm
using PyPlot
using Printf

α = 0.2
β = 1.0
r = 8
m = r+1
n = r
C = MittagLefflerFunctions.Pade.pade_matrix(α, β, m, n)
Ctilde = [ C[:,1:r] C[:,r+2:2r+1] ]
@printf("Condition number of Ctilde = %8.2e\n", cond(Ctilde))
@printf("Condition number of C      = %8.2e\n", cond(C))
F_lu = lu(C)
@printf("Condition number of U      = %8.2e\n", cond(F_lu.U))

R_fix = MLPade(α, β, m, n, :fix)
R_svd = MLPade(α, β, m, n, :svd)
R_lu  = MLPade(α, β, m, n, :lu)

@printf("\t||p_fix-p_svd||_oo = %8.2e\n", norm(R_fix.p-R_svd.p, Inf))
@printf("\t||q_fix-q_svd||_oo = %8.2e\n", norm(R_fix.q-R_svd.q, Inf))
@printf("\t||p_lu-p_svd||_oo  = %8.2e\n", norm(R_lu.p-R_svd.p, Inf))
@printf("\t||q_lu-q_svd||_oo  = %8.2e\n", norm(R_lu.q-R_svd.q, Inf))

figure(1)
x = range(0, 5, length=201)
plot(x, R_fix.(x) - R_svd.(x), x, R_lu.(x) - R_svd.(x))
legend(("fix - svd", "lu - svd"))
grid(true)


text:
# ExoplanetsSysSim/src/summary_plots.jl
## (c) 2015 Darin Ragozzine

try
if !(Pkg.installed("PyPlot") >= v"2.1.0")  # version chosen arbitrarily
  error("You need to upgrade the Optim package for this example.")
end
catch
  error("You need to install the Optim package for this example.")
end

verbose = true

export plot_all_summary_stats

if length(ARGS)<2
   warn("usage:  syssim_plots.jl param_file comp_file")
   println("ASSUMING!!: param_file = demo_param.in AND comp_file = demo_ss.out")

   ## This hack allows this command to be run from inside the REPL or at the
   ## command line, though Julia understandably doesn't like it.
   ARGS=["demo_param.in","demo_ss.out"]
   # exit(255)
end

param_file = convert(String,ARGS[1])
if !isreadable(param_file) || !isfile(param_file)
   warn(string("Can't read ",param_file))
   exit(255)
end

comp_file = convert(String,ARGS[2])
if !isreadable(comp_file) || !isfile(comp_file)
   warn(string("Can't read ",comp_file))
   exit(255)
end

if verbose println("# Loading ExoplanetsSysSim module."); end
using ExoplanetsSysSim
if verbose println("# Loading PyPlot Module; that's how we'll do plots."); end
using PyPlot
pygui(false)  # Tell PyPlot not to use the GUI, but to return objects

if verbose println("# Reading parameter file ", param_file,"."); end
include(param_file)

if verbose println("# Calling setup_sim_param( ", ARGS[3:end], " )."); end
sim_param = setup_sim_param( convert(Array{String,1},ARGS[3:end]) )
if verbose println("# Active parameter values: ", make_vector_of_sim_param(sim_param) ); end

if verbose println("# Loading summary statistics from ", comp_file,"."); end
ss_comp = load_summary_stats(comp_file)

# Maybe?
type PlotSummaryStatistics
  plotnames::Dict{String,Any}  	# For storing names of plots
  plotfilenames::Dict{String,Any}  # For storing filenames of plots
end
function PlotSummaryStatistics()
  PlotSummaryStatistics( Dict{String,Any}(), Dict{String,Any}() )
end

function plot_all_summary_stats(ss::CatalogSummaryStatistics,
      sim_param::SimParam, verbose::Bool = false)
## TODO: Figure out the best filename to output to (user defined?) and implement.

# Print P vs. depth log-log scatter plot
if haskey(ss.stat,"P list") && haskey(ss.stat,"depth list")
	figure()
	## Some searching and outputting directly to a file without popping
	## up a figure window is not obvious, at least using PyPlot. Could use Gadfly.
	## TODO: Allow to output directly to a file without popping up a window
	loglog(ss.stat["P list"],ss.stat["depth list"],".")
    savefig("P_vs_depth.png")
else
	if verbose println("Cant plot P vs. depth because these are not found."); end
end #plot P vs. depth

# Plot multiplicity histogram in semilogy histogram-y scatter plot
if haskey(ss.stat,"num_sys_tranets")
	figure()
	semilogy(ss.stat["num_sys_tranets"],drawstyle="steps")
    savefig("num_sys_tranets.png")
else
	if verbose println("Cant plot multiplicity histogram because not found."); end
end # plot multiplicity histogram

end # plot_all_summary_stats

plot_all_summary_stats(ss_comp, sim_param, verbose)


text:
import Base: ==, *, /, +, -, Broadcast
import SparseArrays: sparse

const SparseOpPureType{BL<:Basis,BR<:Basis} = Operator{BL,BR,<:SparseMatrixCSC}
const SparseOpAdjType{BL<:Basis,BR<:Basis} = Operator{BL,BR,<:Adjoint{<:Number,<:SparseMatrixCSC}}
const SparseOpType{BL<:Basis,BR<:Basis} = Union{SparseOpPureType{BL,BR},SparseOpAdjType{BL,BR}}

"""
    SparseOperator(b1[, b2, data])

Sparse array implementation of Operator.

The matrix is stored as the julia built-in type `SparseMatrixCSC`
in the `data` field.
"""
SparseOperator(b1::Basis, b2::Basis, data) = Operator(b1, b2, SparseMatrixCSC{ComplexF64,Int}(data))
SparseOperator(b1::Basis, b2::Basis, data::SparseMatrixCSC{ComplexF64,Int}) = Operator(b1, b2, data)
SparseOperator(b::Basis, data) = SparseOperator(b, b, data)
SparseOperator(op::DataOperator) = SparseOperator(op.basis_l, op.basis_r, op.data)

SparseOperator(b1::Basis, b2::Basis) = SparseOperator(b1, b2, spzeros(ComplexF64, length(b1), length(b2)))
SparseOperator(b::Basis) = SparseOperator(b, b)

"""
    sparse(op::AbstractOperator)

Convert an arbitrary operator into a [`SparseOperator`](@ref).
"""
sparse(a::AbstractOperator) = throw(ArgumentError("Direct conversion from $(typeof(a)) not implemented. Use sparse(full(op)) instead."))
sparse(a::DataOperator) = Operator(a.basis_l, a.basis_r, sparse(a.data))

function ptrace(op::SparseOpPureType, indices::Vector{Int})
    check_ptrace_arguments(op, indices)
    shape = [op.basis_l.shape; op.basis_r.shape]
    data = ptrace(op.data, shape, indices)
    b_l = ptrace(op.basis_l, indices)
    b_r = ptrace(op.basis_r, indices)
    Operator(b_l, b_r, data)
end

function expect(op::SparseOpPureType{B1,B2}, state::Operator{B2,B2}) where {B1<:Basis,B2<:Basis}
    check_samebases(op, state)
    result = zero(promote_type(eltype(op),eltype(state)))
    @inbounds for colindex = 1:op.data.n
        for i=op.data.colptr[colindex]:op.data.colptr[colindex+1]-1
            result += op.data.nzval[i]*state.data[colindex, op.data.rowval[i]]
        end
    end
    result
end

function permutesystems(rho::SparseOpPureType{B1,B2}, perm::Vector{Int}) where {B1<:CompositeBasis,B2<:CompositeBasis}
    @assert length(rho.basis_l.bases) == length(rho.basis_r.bases) == length(perm)
    @assert isperm(perm)
    shape = [rho.basis_l.shape; rho.basis_r.shape]
    data = permutedims(rho.data, shape, [perm; perm .+ length(perm)])
    SparseOperator(permutesystems(rho.basis_l, perm), permutesystems(rho.basis_r, perm), data)
end

identityoperator(::Type{T}, b1::Basis, b2::Basis) where {T<:DataOperator} = SparseOperator(b1, b2, sparse(ComplexF64(1)*I, length(b1), length(b2)))
identityoperator(b1::Basis, b2::Basis) = identityoperator(DataOperator, b1, b2)
identityoperator(b::Basis) = identityoperator(b, b)

"""
    diagonaloperator(b::Basis)

Create a diagonal operator of type [`SparseOperator`](@ref).
"""
function diagonaloperator(b::Basis, diag::Vector{T}) where T <: Number
  @assert 1 <= length(diag) <= length(b)
  SparseOperator(b, sparse(Diagonal(Vector{ComplexF64}(diag))))
end

# Fast in-place multiplication implementations
mul!(result::DenseOpType{B1,B3},M::SparseOpType{B1,B2},b::DenseOpType{B2,B3},alpha,beta) where {B1<:Basis,B2<:Basis,B3<:Basis} = (gemm!(alpha,M.data,b.data,beta,result.data); result)
mul!(result::DenseOpType{B1,B3},a::DenseOpType{B1,B2},M::SparseOpType{B2,B3},alpha,beta) where {B1<:Basis,B2<:Basis,B3<:Basis} = (gemm!(alpha,a.data,M.data,beta,result.data); result)
mul!(result::Ket{B1},M::SparseOpPureType{B1,B2},b::Ket{B2},alpha,beta) where {B1<:Basis,B2<:Basis} = (gemv!(alpha,M.data,b.data,beta,result.data); result)
mul!(result::Bra{B2},b::Bra{B1},M::SparseOpPureType{B1,B2},alpha,beta) where {B1<:Basis,B2<:Basis} = (gemv!(alpha,b.data,M.data,beta,result.data); result)


text:
using Random

"""
    poisson(baserate, theta, deltav, v; dt::Real, rng::AbstractRNG)
    poisson(baserate::AbstractArray{<:Real}, theta::AbstractArray{<:Real}, deltav::AbstractArray{<:Real}, v::AbstractArray{<:Real}; dt::Real, rng::AbstractRNG)
    poisson(baserate::CuVector{<:Real}, theta::CuVector{<:Real}, deltav::CuVector{<:Real}, v::CuVector{<:Real}; dt::Real, rng::AbstractRNG)

Evaluate inhomogeneous Poisson process threshold functions.
Modeled as

``X < \\mathrm{d}t \\rho_0 \\exp\\left(\\frac{v - \\Theta}{\\Delta_u}\\right)``

where ``X \\sim \\mathrm{Unif}([0, 1])``.

Use `CuVector` instead of `Vector` to evaluate on GPU.

# Fields
- `baserate`: base line firing rate
- `theta`: threshold potential
- `deltav`: potential resolution
- `v`: current membrane potential
- `dt`: simulation timestep
- `rng`: random number generation
"""
function poisson(baserate, theta, deltav, v; dt::Real, rng::AbstractRNG = Random.GLOBAL_RNG)
    rho = baserate * exp((v - theta) / deltav)

    return rand(rng) < rho * dt
end
function poisson(baserate::AbstractArray{<:Real}, theta::AbstractArray{<:Real}, deltav::AbstractArray{<:Real}, v::AbstractArray{<:Real}; dt::Real, rng::AbstractRNG = Random.GLOBAL_RNG)
    rho = baserate .* exp.((v .- theta) ./ deltav)

    return rand(rng, length(rho)) .< rho .* dt
end
function poisson(baserate::CuVector{<:Real}, theta::CuVector{<:Real}, deltav::CuVector{<:Real}, v::CuVector{<:Real}; dt::Real, rng::AbstractRNG = Random.GLOBAL_RNG)
    rho = baserate .* exp.((v .- theta) ./ deltav)

    return CuArrays.rand(rng, length(rho)) .< rho .* dt
end


text:
using Pipe
using Test

_macroexpand(q) = macroexpand(Main, q)

rml! = Base.remove_linenums!
# performs linenum removal and temp variable replacing to avoid different names of temp variables in different julia versions
stringify_expr(e::Expr) = replace(string(rml!(e)), r"##\d{3}"=>"##000")
pipe_equals(e1::Expr, e2::Expr) = stringify_expr(_macroexpand(e1)) == stringify_expr(e2)

#No change to nonpipes functionality
@test _macroexpand( :(@pipe a) ) == :a #doesn't change single inputs
@test _macroexpand( :(@pipe b(a)) ) == :(b(a)) #doesn't change inputs that a function applications

#Compatable with Julia 1.3 piping functionality
@test _macroexpand( :(@pipe a|>b) ) == :(b(a)) #basic
@test _macroexpand( :(@pipe a|>b|>c) ) == :(c(b(a)))  #Keeps chaining 3
@test _macroexpand( :(@pipe a|>b|>c|>d) ) == :(d(c(b(a)))) #Keeps chaining 4

@test _macroexpand( :(@pipe a|>b(x)) ) == :((b(x))(a))  #applying to function calls returning functions
@test _macroexpand( :(@pipe a(x)|>b ) ) == :(b(a(x)))   #feeding functioncall results on wards

@test _macroexpand(:(@pipe 1|>a)) ==:(a(1)) #Works with literals (int)
@test _macroexpand(:(@pipe "foo"|>a)) == :(a("foo")) #Works with literal (string)
@test _macroexpand( :(@pipe a|>bb[2])) == :((bb[2])(a)) #Should work with RHS that is a array reference

#Marked locations
@test _macroexpand( :(@pipe a |> _)) == :(a) #Identity works
@test _macroexpand( :(@pipe a |> _[b])) == :(a[b]) #Indexing works

@test _macroexpand( :(@pipe a|>b(_) ) ) == :(b(a)) #Marked location only
@test _macroexpand( :(@pipe a|>b(x,_) ) ) == :(b(x,a)) # marked 2nd (and last)
@test _macroexpand( :(@pipe a|>b(_,x) ) ) == :(b(a,x)) # marked first
@test _macroexpand( :(@pipe a|>b(_,_) ) ) == :(b(a,a)) # marked double (Not certain if this is a good idea)
@test _macroexpand( :(@pipe a|>bb[2](x,_))) == :((bb[2])(x,a)) #Should work with RHS that is a array reference

#Macros and blocks
macro testmacro(arg, n)
    esc(:($arg + $n))
end
@test _macroexpand( :(@pipe a |> @testmacro _ 3 ) ) == :(a + 3) # Can pipe into macros
@test _macroexpand( :(@pipe a |> begin b = _; c + b + _ end )) == :(
                                 begin b = a; c + b + a end)

#marked Unpacking
@test _macroexpand( :(@pipe a|>b(_...) ) ) == :(b(a...)) # Unpacking
@test _macroexpand( :(@pipe a|>bb[2](_...))) == :((bb[2])(a...)) #Should work with RHS of arry ref and do unpacking

#Mixing modes
@test _macroexpand( :(@pipe a|>b|>c(_) ) ) == :(c(b(a)))
@test _macroexpand( :(@pipe a|>b(x,_)|>c|>d(_,y) ) ) == :(d(c(b(x,a)),y))
@test _macroexpand( :(@pipe a|>b(xb,_)|>c|>d(_,xd)|>e(xe) |>f(xf,_,yf)|>_[i] ) ) == :(f(xf,(e(xe))(d(c(b(xb,a)),xd)),yf)[i]) #Very Complex

# broadcasting
vars = 1:10 .|> y->gensym() # Julia < 1.3 changes how Symbols are stringified so we compute the representation here
@test pipe_equals(:(@pipe 1:10 .|> _*2 ), :(1:10 .|> $(vars[1])->$(vars[1]) * 2))
@test pipe_equals(:(@pipe 1:10 .|> fn ), :(1:10 .|> $(vars[2])->fn($(vars[2]))))
@test pipe_equals(:(@pipe a .|> fn .|> _*2 ), :(a .|> ($(vars[3])->fn($(vars[3]))) .|> ($(vars[4])->$(vars[4])*2)))
@test pipe_equals(:(@pipe a .|> fn |> _*2 ), :((a .|> $(vars[5])->fn($(vars[5]))) * 2))
@test pipe_equals(:(@pipe [1,2,2] |> atan.([10,20,30], _) ), :(atan.([10,20,30], [1,2,2])))
@test pipe_equals(:(@pipe [1,2,2] .|> atan.([10,20,30], _) ), :([1,2,2] .|> $(vars[6])->atan.([10,20,30], $(vars[6]))))
@test pipe_equals(:(@pipe fn |> _.(1:2) ), :(fn.(1:2)))
@test pipe_equals(:(@pipe fn .|> _.(1:2) ), :(fn .|> $(vars[7])->$(vars[7]).(1:2)))

@test pipe_equals(:(@pipe [true,false] .|> ! ), :([true, false] .|> $(vars[8])->!$(vars[8])))
@test pipe_equals(:(@pipe [1, 2] |> .+(_, x) ), :([1, 2] .+ x))
@test pipe_equals(:(@pipe [1, 2] |>  _ .+ x ), :([1, 2] .+ x))
@test pipe_equals(:(@pipe [1, 2] .|> .+(_, x) ), :([1, 2] .|> $(vars[9])->$(vars[9]).+x))
@test pipe_equals(:(@pipe [1, 2] .|>  _ .+ x ), :([1, 2] .|> $(vars[10])->$(vars[10]).+x))