text:
# math functionality corresponding to SpecialFunctions.jl

## error

@device_override SpecialFunctions.erf(x::Float64) = ccall("extern __nv_erf", llvmcall, Cdouble, (Cdouble,), x)
@device_override SpecialFunctions.erf(x::Float32) = ccall("extern __nv_erff", llvmcall, Cfloat, (Cfloat,), x)

@device_override SpecialFunctions.erfinv(x::Float64) = ccall("extern __nv_erfinv", llvmcall, Cdouble, (Cdouble,), x)
@device_override SpecialFunctions.erfinv(x::Float32) = ccall("extern __nv_erfinvf", llvmcall, Cfloat, (Cfloat,), x)

@device_override SpecialFunctions.erfc(x::Float64) = ccall("extern __nv_erfc", llvmcall, Cdouble, (Cdouble,), x)
@device_override SpecialFunctions.erfc(x::Float32) = ccall("extern __nv_erfcf", llvmcall, Cfloat, (Cfloat,), x)

@device_override SpecialFunctions.erfcinv(x::Float64) = ccall("extern __nv_erfcinv", llvmcall, Cdouble, (Cdouble,), x)
@device_override SpecialFunctions.erfcinv(x::Float32) = ccall("extern __nv_erfcinvf", llvmcall, Cfloat, (Cfloat,), x)

@device_override SpecialFunctions.erfcx(x::Float64) = ccall("extern __nv_erfcx", llvmcall, Cdouble, (Cdouble,), x)
@device_override SpecialFunctions.erfcx(x::Float32) = ccall("extern __nv_erfcxf", llvmcall, Cfloat, (Cfloat,), x)

## gamma function

@device_override SpecialFunctions.lgamma(x::Float64) = ccall("extern __nv_lgamma", llvmcall, Cdouble, (Cdouble,), x)
@device_override SpecialFunctions.lgamma(x::Float32) = ccall("extern __nv_lgammaf", llvmcall, Cfloat, (Cfloat,), x)

@device_override SpecialFunctions.gamma(x::Float64) = ccall("extern __nv_tgamma", llvmcall, Cdouble, (Cdouble,), x)
@device_override SpecialFunctions.gamma(x::Float32) = ccall("extern __nv_tgammaf", llvmcall, Cfloat, (Cfloat,), x)

## Bessel

@device_override SpecialFunctions.besselj0(x::Float64) = ccall("extern __nv_j0", llvmcall, Cdouble, (Cdouble,), x)
@device_override SpecialFunctions.besselj0(x::Float32) = ccall("extern __nv_j0f", llvmcall, Cfloat, (Cfloat,), x)

@device_override SpecialFunctions.besselj1(x::Float64) = ccall("extern __nv_j1", llvmcall, Cdouble, (Cdouble,), x)
@device_override SpecialFunctions.besselj1(x::Float32) = ccall("extern __nv_j1f", llvmcall, Cfloat, (Cfloat,), x)

@device_override SpecialFunctions.besselj(n::Int32, x::Float64) = ccall("extern __nv_jn", llvmcall, Cdouble, (Int32, Cdouble), n, x)
@device_override SpecialFunctions.besselj(n::Int32, x::Float32) = ccall("extern __nv_jnf", llvmcall, Cfloat, (Int32, Cfloat), n, x)

@device_override SpecialFunctions.bessely0(x::Float64) = ccall("extern __nv_y0", llvmcall, Cdouble, (Cdouble,), x)
@device_override SpecialFunctions.bessely0(x::Float32) = ccall("extern __nv_y0f", llvmcall, Cfloat, (Cfloat,), x)

@device_override SpecialFunctions.bessely1(x::Float64) = ccall("extern __nv_y1", llvmcall, Cdouble, (Cdouble,), x)
@device_override SpecialFunctions.bessely1(x::Float32) = ccall("extern __nv_y1f", llvmcall, Cfloat, (Cfloat,), x)

@device_override SpecialFunctions.bessely(n::Int32, x::Float64) = ccall("extern __nv_yn", llvmcall, Cdouble, (Int32, Cdouble), n, x)
@device_override SpecialFunctions.bessely(n::Int32, x::Float32) = ccall("extern __nv_ynf", llvmcall, Cfloat, (Int32, Cfloat), n, x)


text:
module SQLiteFs

using SQLite
using DataFrames
mutable struct FStatus
    db::SQLite.DB
    dir::String
    ino::Int
    exception::Union{Nothing,Exception}
    FStatus(db) = new(db, DIRROOT, ROOT_INO, nothing)
end

include("layout.jl")
include("cstruct.jl")
include("fuseapi.jl")
include("fusebridge.jl")
include("initdb.jl")
include("filesystem.jl")

include("examplefs.jl")

function _umask()
    umask = mktempdir() do tmp
        s = mkdir(joinpath(tmp, "stest"))
        st = stat(s)
        ~st.mode & X_UGO
    end
    umask
end

function __init__()
    eval(:(const UMASK = _umask()))
end
end

ext: jl


text:
# test 3 hardening tests
using Distributed;
addprocs(20);
@everywhere include("loadMod.jl");
@everywhere const GUROBI_ENV = Gurobi.Env();

@everywhere function obtainStat(fData,NN,T,solSDDP)
    # obtain the load mismatch and generation
    spSDDP = Dict();
    sqSDDP = Dict();
    LSDDP = Dict();
    spList = [];
    sqList = [];
    LList = [];
    for i in 1:NN
        spSDDP[i] = Dict();
        sqSDDP[i] = Dict();
        LSDDP[i] = Dict();
        for g in fData.genIDList
            spSDDP[i][g] = Dict();
            sqSDDP[i][g] = Dict();
        end
        for ib in fData.IDList
            LSDDP[i][ib] = Dict();
        end

        dSDDPList = [solSDDP[i][j][2] for j in 1:length(solSDDP[i])];
        push!(dSDDPList,T+1);
        for j in 1:length(solSDDP[i])
            # obtain the inventory & actual injection to network until the next solution
            for t in dSDDPList[j]:(dSDDPList[j+1] - 1)
                for g in fData.genIDList
                    spSDDP[i][g][t] = solSDDP[i][j][1].sp[g,t];
                    sqSDDP[i][g][t] = solSDDP[i][j][1].sq[g,t];
                end
                for ib in fData.IDList
                    LSDDP[i][ib][t] = solSDDP[i][j][1].lp[ib,t] + solSDDP[i][j][1].lq[ib,t];
                end
            end
        end
        # obtain the statistics to show total generation and power mismatch
        push!(spList,sum(sum(spSDDP[i][g][t] for t in 1:T) for g in fData.genIDList));
        push!(sqList,sum(sum(sqSDDP[i][g][t] for t in 1:T) for g in fData.genIDList));
        push!(LList,sum(sum(LSDDP[i][ib][t] for t in 1:T) for ib in fData.IDList));
    end
    return spList,sqList,LList;
end

N = 5;

caseList = [13,33,123];
Δt = 0.25;
iterMax = 20;
T = 24;
NN = 1000;
pathTrain = load("pathHist_600.jld");
pathListDRaw = load("pathHist_1000.jld");

for ci in 1:length(caseList)
    pathDict = pathTrain["pathDict"][ci][T];
    pathDictTest = pathListDRaw["pathDict"][ci][T];

    dataList = Dict();
    τ = Int64(1/6*T);
    for j in procs()
        remotecall_fetch(readInData,j,ci,caseList,T,τ,1e4);
        #remotecall_fetch(readInData_old,j,T,[(2,9),(8,12),(10,13)],10000,0);
    end
    cutDictPG = preGen(T, Δt, N, iterMax, false, Dict(), []);
    # cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(τ, T, Δt, N, true, false,
    #     20, 20, cutDictPG, false, 0, [], 0, pathDict);
    cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(T, Δt, N, false, false,
        max(Int64(round(500/N)),20), max(Int64(round(500/N)),20), cutDictPG, false, 0, []);
    solSDDP, LBSDDP, costSDDP = exeForward(T, Δt, NN, false, pathDictTest);
    spList,sqList,LList = obtainStat(fData,NN,T,solSDDP);

    dataList["NoD"] = [LBHist,UBHist,UBuHist,UBlHist,timeHist,LBSDDP,costSDDP,spList,sqList,LList];

    for ω in keys(pDistr.ωDistrn)
    #for ω in [(2,9),(8,12),(10,13)]
        cutDictPG = preGen(T, Δt, N, iterMax, false, Dict(), pDistr.ωDict[ω]);
        # cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(τ, T, Δt, N, true, false,
        #     20, 20, cutDictPG, false, 0, [ω], 0, pathDict);
        cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(T, Δt, N, false, false,
            max(Int64(round(500/N)),20), max(Int64(round(500/N)),20), cutDictPG, false, 0, pDistr.ωDict[ω]);
        solSDDP, LBSDDP, costSDDP = exeForward(T, Δt, NN, false, pathDictTest, pDistr.ωDict[ω]);
        spList,sqList,LList = obtainStat(fData,NN,T,solSDDP);

        dataList[ω] = [LBHist,UBHist,UBuHist,UBlHist,timeHist,LBSDDP,costSDDP,spList,sqList,LList];
        save("hardResults_$(ci).jld","data",dataList);
    end
end


text:
"""
    LabeledValue{T}

Value of type `T` that is associated with a dictionary of value labels.
If a value `v` is not a key in the dictionary,
then `string(v)` is taken as the label.
See also [`LabeledArray`](@ref).

The value underlying a `LabeledValue` can be accessed via [`unwrap`](@ref).
The label can be obtained by converting `LabeledValue` to `String`
or calling [`labels`](@ref) (notice the `s` at the end).

Comparison operators `==`, `isequal`, `<`, `isless` and `isapprox`
compare the underlying value of type `T`.
An exception is that when a `LabeledValue` and a string are compared with `==`,
the comparison is based on the label.

# Examples
```jldoctest
julia> lbls = Dict{Union{Int,Missing},String}(0=>"a", 1=>"a");

julia> v0 = LabeledValue(0, lbls)
0 => a

julia> v1 = LabeledValue(1, lbls)
1 => a

julia> vm = LabeledValue(missing, lbls)
missing => missing

julia> v0 == v1
false

julia> v1 == 1
true

julia> v1 == "a"
true

julia> isequal(vm, missing)
true
```
"""
struct LabeledValue{T}
    value::T
    labels::Dict{T, String}
end

# A value is not guaranteed to have a label defined in labels
function _getlabel(x::LabeledValue)
    lbl = get(x.labels, x.value, nothing)
    lbl === nothing && return string(x.value)
    return lbl
end

# Comparison between LabeledValues ignores labels
Base.:(==)(x::LabeledValue, y::LabeledValue) = x.value == y.value
Base.isequal(x::LabeledValue, y::LabeledValue) = isequal(x.value, y.value)
Base.:(<)(x::LabeledValue, y::LabeledValue) = x.value < y.value
Base.isless(x::LabeledValue, y::LabeledValue) = isless(x.value, y.value)
Base.isapprox(x::LabeledValue, y::LabeledValue; kwargs...) =
    isapprox(x.value, y.value; kwargs...)

# Comparison with String is interpreted as comparing labels
# Comparison with any other type is interpreted as comparing values
Base.:(==)(x::LabeledValue, y) = x.value == y
Base.:(==)(x, y::LabeledValue) = x == y.value
Base.:(==)(x::LabeledValue, y::AbstractString) = _getlabel(x) == y
Base.:(==)(x::AbstractString, y::LabeledValue) = x == _getlabel(y)
Base.:(==)(::LabeledValue, ::Missing) = missing
Base.:(==)(::Missing, ::LabeledValue) = missing

Base.isequal(x::LabeledValue, y) = isequal(x.value, y)
Base.isequal(x, y::LabeledValue) = isequal(x, y.value)
Base.isequal(x::LabeledValue, y::Missing) = isequal(x.value, y)
Base.isequal(x::Missing, y::LabeledValue) = isequal(x, y.value)
Base.:(<)(x::LabeledValue, y) = x.value < y
Base.:(<)(x, y::LabeledValue) = x < y.value
Base.:(<)(x::LabeledValue, y::Missing) = x.value < y
Base.:(<)(x::Missing, y::LabeledValue) = x < y.value
Base.isless(x::LabeledValue, y) = isless(x.value, y)
Base.isless(x, y::LabeledValue) = isless(x, y.value)
Base.isless(x::LabeledValue, y::Missing) = isless(x.value, y)
Base.isless(x::Missing, y::LabeledValue) = isless(x, y.value)
Base.isapprox(x::LabeledValue, y; kwargs...) = isapprox(x.value, y; kwargs...)
Base.isapprox(x, y::LabeledValue; kwargs...) = isapprox(x, y.value; kwargs...)

Base.hash(x::LabeledValue, h::UInt=zero(UInt)) = hash(x.value, h)

"""
    unwrap(x::LabeledValue)

Get the value wrapped by `x`.
"""
unwrap(x::LabeledValue) = x.value

"""
    labels(x::LabeledValue)

Return the label associated with `x`.
"""
labels(x::LabeledValue) = _getlabel(x)

Base.show(io::IO, x::LabeledValue) = print(io, _getlabel(x))
Base.show(io::IO, ::MIME"text/plain", x::LabeledValue) =
    print(io, x.value, " => ", _getlabel(x))

Base.convert(::Type{String}, x::LabeledValue) = _getlabel(x)

"""
    LabeledArray{V, N, T<:LabeledValue} <: AbstractArray{T, N}

`N`-dimensional dense array with elements associated with labels.

`LabeledArray` provides functionality that is similar to
what value labels achieve in statistical software such as Stata.
When printed to REPL, a `LabeledArray` just looks like an array of labels.
Yet, only the underlying values of type `V` are stored in an `Array`.
The associated labels are looked up
from a dictionary of type `Dict{V, String}`.
If a value `v` is not a key in the dictionary,
then `string(v)` is taken as the label.
The elements of type [`LabeledValue`](@ref)
are only constructed lazily when retrieved.

The array of values underlying a `LabeledArray`
can be accessed with [`refarray`](@ref);
while an iterator over the labels for each element
is returned by [`labels`](@ref).

Equality comparison as defined by `==` involving a `LabeledArray`
only compares the underlying values
unless the element type of the other array is a subtype of `AbstractString`.
The labels are used for comparison in the latter case.

# Examples
```jldoctest
julia> lbls1 = Dict(1=>"a", 2=>"b");

julia> lbls2 = Dict(1.0=>"p", 2.0=>"q");

julia> x = LabeledArray([0, 1, 2], lbls1)
3-element LabeledVector{Int64, LabeledValue{Int64}}:
 0 => 0
 1 => a
 2 => b

julia> y = LabeledArray([0.0, 1.0, 2.0], lbls2)
3-element LabeledVector{Float64, LabeledValue{Float64}}:
 0.0 => 0.0
 1.0 => p
 2.0 => q

julia> x == y
true

julia> x == 0:2
true

julia> x == ["0", "a", "b"]
true
```
"""
struct LabeledArray{V, N, T<:LabeledValue} <: AbstractArray{T, N}
    values::Array{V, N}
    labels::Dict{V, String}
    function LabeledArray{V,N}(values::Array{V,N}, labels::Dict{V,String}) where {V,N}
        V <: AbstractString && throw(ArgumentError("values of type $V are not accepted"))
        return new{V,N,LabeledValue{V}}(values, labels)
    end
end

"""
    LabeledVector{V, T}

Alias for [`LabeledArray{V, 1, T}`](@ref).
"""
const LabeledVector{V, T} = LabeledArray{V, 1, T}

LabeledArray(values::Array{V,N}, labels::Dict{V,String}) where {V,N} =
    LabeledArray{V,N}(values, labels)

Base.size(x::LabeledArray) = size(x.values)
Base.IndexStyle(::Type{<:LabeledArray}) = IndexLinear()

Base.@propagate_inbounds function Base.getindex(x::LabeledArray, i::Int)
    val = x.values[i]
    return LabeledValue(val, x.labels)
end

"""
    refarry(x::LabeledArray)
    refarray(x::SubArray{<:Any, <:Any, <:LabeledArray})
    refarray(x::Base.ReshapedArray{<:Any, <:Any, <:LabeledArray})
    refarray(x::AbstractArray{<:LabeledValue})

Return the array of values underlying a [`LabeledArray`](@ref).
"""
refarray(x::LabeledArray) = x.values
refarray(x::SubArray{<:Any, <:Any, <:LabeledArray}) =
    view(parent(x).values, x.indices...)
refarray(x::Base.ReshapedArray{<:Any, <:Any, <:LabeledArray}) =
    reshape(parent(x).values, size(x))
refarray(x::AbstractArray{<:LabeledValue}) = collect(v.value for v in x)

struct LabelIterator{A, N} <: AbstractArray{String, N}
    a::A
    LabelIterator(a::AbstractArray{<:LabeledValue, N}) where N = new{typeof(a), N}(a)
end

"""
    labels(x::AbstractArray{<:LabeledValue})

Return an iterator over the labels for each element in `x`.
The returned object is a subtype of `AbstractArray` with the same size of `x`.

The iterator can be used to collect the labels while discarding the underlying values.

# Examples
```jldoctest
julia> x = LabeledArray([1, 2, 3], Dict(1=>"a", 2=>"b"))
3-element LabeledVector{Int64, LabeledValue{Int64}}:
 1 => a
 2 => b
 3 => 3

julia> lbls = labels(x)
3-element ReadStatTables.LabelIterator{LabeledVector{Int64, LabeledValue{Int64}}, 1}:
 "a"
 "b"
 "3"

julia> collect(lbls)
3-element Vector{String}:
 "a"
 "b"
 "3"

julia> CategoricalArray(lbls)
3-element CategoricalArray{String,1,UInt32}:
 "a"
 "b"
 "3"
```
"""
labels(x::AbstractArray{<:LabeledValue}) = LabelIterator(x)

Base.size(lbls::LabelIterator) = size(lbls.a)
Base.IndexStyle(::Type{<:LabelIterator{A}}) where A = IndexStyle(A)

Base.@propagate_inbounds function Base.getindex(lbls::LabelIterator{<:LabeledArray}, i::Int)
    val = lbls.a.values[i]
    lbl = get(lbls.a.labels, val, nothing)
    lbl === nothing && return string(val)
    return lbl
end

Base.@propagate_inbounds Base.getindex(lbls::LabelIterator{<:AbstractArray}, i::Int) =
    _getlabel(lbls.a[i])

Base.:(==)(x::LabeledArray, y::LabeledArray) = x.values == y.values
Base.:(==)(x::LabeledArray, y::AbstractArray) = x.values == y
Base.:(==)(x::AbstractArray, y::LabeledArray) = x == y.values
Base.:(==)(x::LabeledArray, y::AbstractArray{<:AbstractString}) = labels(x) == y
Base.:(==)(x::AbstractArray{<:AbstractString}, y::LabeledArray) = x == labels(y)

Base.copy(x::LabeledArray{V,N}) where {V,N} =
    LabeledArray{V,N}(copy(x.values), copy(x.labels))

Base.convert(::Type{Array{V1,N}}, x::LabeledArray{V2,N}) where {V1,V2,N} =
    convert(Array{V1}, x.values)

# Assume VERSION >= v"1.3.0"
function compact_type_str(::Type{LabeledValue{V}}) where V
    str = V >: Missing ? string(nonmissingtype(V)) * "?" : string(V)
    return replace("Labeled{" * str * "}", "Union" => "U")
end


text:
module AtomShell

using Sockets
using WebIO
using RelocatableFolders

abstract type Shell end

include("install.jl")
include("process.jl")
include("window.jl")
include("webio.jl")

end


text:
# FASTA File Format
# =================

module FASTA

import Automa
import Automa.RegExp: @re_str
import Bio: Bio, isfilled
import Bio.Exceptions: missingerror
import BufferedStreams

include("record.jl")
include("index.jl")
include("reader.jl")
include("writer.jl")

end


text:
@testset "Graph coloring" begin

    function normal_49_states(;
        reverse_constraint = false,
        tests = true,
        time_limit = Inf,
        simplify = true,
    )
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "keep_logs" => true,
            "logging" => [],
            "time_limit" => time_limit,
            "simplify" => simplify,
            "seed" => 1,
        ))
        num_colors = 20

        @variable(m, 1 <= states[1:50] <= num_colors, Int)
        @variable(m, 1 <= max_color <= num_colors, Int)

        washington,
        montana,
        maine,
        north_dakota,
        south_dakota,
        wyoming,
        wisconsin,
        idaho,
        vermont,
        minnesota,
        oregon,
        new_hampshire,
        iowa,
        massachusetts,
        nebraska,
        new_york,
        pennsylvania,
        connecticut,
        rhode_island,
        new_jersey,
        indiana,
        nevada,
        utah,
        california,
        ohio,
        illinois,
        washington_dc,
        delaware,
        west_virginia,
        maryland,
        colorado,
        kentucky,
        kansas,
        virginia,
        missouri,
        arizona,
        oklahoma,
        north_carolina,
        tennessee,
        texas,
        new_mexico,
        alabama,
        mississippi,
        georgia,
        south_carolina,
        arkansas,
        louisiana,
        florida,
        michigan = states

        @constraint(m, washington != oregon)
        @constraint(m, washington != idaho)
        @constraint(m, oregon != idaho)
        @constraint(m, oregon != nevada)
        @constraint(m, oregon != california)
        @constraint(m, california != nevada)
        @constraint(m, california != arizona)
        @constraint(m, nevada != idaho)
        @constraint(m, nevada != utah)
        @constraint(m, nevada != arizona)
        @constraint(m, idaho != montana)
        @constraint(m, idaho != wyoming)
        @constraint(m, idaho != utah)
        @constraint(m, utah != wyoming)
        @constraint(m, utah != colorado)
        @constraint(m, utah != new_mexico)
        @constraint(m, utah != arizona)
        @constraint(m, arizona != colorado)
        @constraint(m, arizona != new_mexico)
        @constraint(m, montana != north_dakota)
        @constraint(m, montana != south_dakota)
        @constraint(m, montana != wyoming)
        @constraint(m, wyoming != south_dakota)
        @constraint(m, wyoming != nebraska)
        @constraint(m, wyoming != colorado)
        @constraint(m, colorado != nebraska)
        @constraint(m, colorado != kansas)
        @constraint(m, colorado != oklahoma)
        @constraint(m, colorado != new_mexico)
        @constraint(m, new_mexico != oklahoma)
        @constraint(m, new_mexico != texas)
        @constraint(m, north_dakota != minnesota)
        @constraint(m, north_dakota != south_dakota)
        @constraint(m, south_dakota != minnesota)
        @constraint(m, south_dakota != iowa)
        @constraint(m, south_dakota != nebraska)
        @constraint(m, nebraska != iowa)
        @constraint(m, nebraska != missouri)
        @constraint(m, nebraska != kansas)
        @constraint(m, kansas != missouri)
        @constraint(m, kansas != oklahoma)
        @constraint(m, oklahoma != arkansas)
        @constraint(m, oklahoma != texas)
        @constraint(m, texas != arkansas)
        @constraint(m, texas != louisiana)
        @constraint(m, minnesota != wisconsin)
        @constraint(m, minnesota != iowa)
        @constraint(m, iowa != wisconsin)
        @constraint(m, iowa != illinois)
        @constraint(m, iowa != missouri)
        @constraint(m, missouri != illinois)
        @constraint(m, missouri != kentucky)
        @constraint(m, missouri != tennessee)
        @constraint(m, missouri != arkansas)
        @constraint(m, arkansas != tennessee)
        @constraint(m, arkansas != mississippi)
        @constraint(m, arkansas != louisiana)
        @constraint(m, louisiana != mississippi)
        @constraint(m, wisconsin != illinois)
        @constraint(m, wisconsin != michigan)
        @constraint(m, illinois != indiana)
        @constraint(m, illinois != kentucky)
        @constraint(m, kentucky != indiana)
        @constraint(m, kentucky != ohio)
        @constraint(m, kentucky != west_virginia)
        @constraint(m, kentucky != virginia)
        @constraint(m, kentucky != tennessee)
        @constraint(m, tennessee != virginia)
        @constraint(m, tennessee != north_carolina)
        @constraint(m, tennessee != georgia)
        @constraint(m, tennessee != alabama)
        @constraint(m, tennessee != mississippi)
        @constraint(m, mississippi != alabama)
        @constraint(m, michigan != indiana)
        @constraint(m, michigan != ohio)
        @constraint(m, indiana != ohio)
        @constraint(m, alabama != georgia)
        @constraint(m, alabama != florida)
        @constraint(m, ohio != pennsylvania)
        @constraint(m, ohio != west_virginia)
        @constraint(m, maine != new_hampshire)
        @constraint(m, new_hampshire != vermont)
        @constraint(m, new_hampshire != massachusetts)
        @constraint(m, vermont != massachusetts)
        @constraint(m, vermont != new_york)
        @constraint(m, massachusetts != new_york)
        @constraint(m, massachusetts != rhode_island)
        @constraint(m, massachusetts != connecticut)
        @constraint(m, connecticut != rhode_island)
        @constraint(m, new_york != connecticut)
        @constraint(m, new_york != pennsylvania)
        @constraint(m, new_york != new_jersey)
        @constraint(m, pennsylvania != new_jersey)
        @constraint(m, pennsylvania != delaware)
        @constraint(m, pennsylvania != maryland)
        @constraint(m, pennsylvania != west_virginia)
        @constraint(m, new_jersey != delaware)
        @constraint(m, maryland != washington_dc)
        @constraint(m, maryland != west_virginia)
        @constraint(m, maryland != virginia)
        @constraint(m, washington_dc != virginia)
        @constraint(m, west_virginia != virginia)
        @constraint(m, virginia != north_carolina)
        @constraint(m, north_carolina != south_carolina)
        @constraint(m, south_carolina != georgia)
        @constraint(m, georgia != florida)

        if reverse_constraint
            @constraint(m, max_color .>= states)
        else
            @constraint(m, states .<= max_color)
        end

        @objective(m, Min, max_color)

        optimize!(m)

        status = JuMP.termination_status(m)

        com = CS.get_inner_model(m)

        CS.save_logs(
            com,
            "graph_color_optimize.json",
            :states => states,
            :max_color => max_color,
        )
        rm("graph_color_optimize.json")

        if tests
            @test status == MOI.OPTIMAL
            @test com.best_sol == 4 == JuMP.objective_value(m) == JuMP.objective_bound(m)
            @test maximum([JuMP.value(var) for var in states]) == 4 == JuMP.value(max_color)
            @test 0 <= MOI.get(m, MOI.SolveTime()) < 5
        end
        return com, m
    end

    @testset "49 US states + DC time limit" begin
        com1, m1 = normal_49_states(; time_limit = 0.001, tests = false, simplify = false)
        com2, m2 = normal_49_states(; simplify = false)
        @test test_string([constraint.indices for constraint in com1.constraints]) == test_string([constraint.indices for constraint in com2.constraints])
        info_1 = com1.info
        info_2 = com2.info
        @test info_1.pre_backtrack_calls == info_2.pre_backtrack_calls
        @test JuMP.termination_status(m1) == MOI.TIME_LIMIT
        # can take longer as preprocessing takes time
        @test 0 <= MOI.get(m1, MOI.SolveTime()) < 0.1
    end

    @testset "49 US states + DC" begin
        com1, _ = normal_49_states(; reverse_constraint = true)
        com2, _ = normal_49_states()
        info_1 = com1.info
        info_2 = com2.info
        @test info_1.pre_backtrack_calls == info_2.pre_backtrack_calls
        @test info_1.backtrack_fixes == info_2.backtrack_fixes
        @test info_1.in_backtrack_calls == info_2.in_backtrack_calls
        @test info_1.backtrack_reverses == info_2.backtrack_reverses
        logs_1 = CS.get_logs(com1)
        logs_2 = CS.get_logs(com2)
        @test CS.sanity_check_log(logs_1[:tree])
        @test CS.same_logs(logs_1[:tree], logs_2[:tree])
    end

    @testset "49 US states + DC without sorting + some states same color" begin
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "keep_logs" => true,
            "backtrack_sorting" => false,
            "logging" => [],
        ))
        num_colors = 8
        @variable(m, 1 <= max_color <= num_colors, Int)

        @variable(m, 1 <= states[1:50] <= num_colors, Int)
        washington,
        montana,
        maine,
        north_dakota,
        south_dakota,
        wyoming,
        wisconsin,
        idaho,
        vermont,
        minnesota,
        oregon,
        new_hampshire,
        iowa,
        massachusetts,
        nebraska,
        new_york,
        pennsylvania,
        connecticut,
        rhode_island,
        new_jersey,
        indiana,
        nevada,
        utah,
        california,
        ohio,
        illinois,
        washington_dc,
        delaware,
        west_virginia,
        maryland,
        colorado,
        kentucky,
        kansas,
        virginia,
        missouri,
        arizona,
        oklahoma,
        north_carolina,
        tennessee,
        texas,
        new_mexico,
        alabama,
        mississippi,
        georgia,
        south_carolina,
        arkansas,
        louisiana,
        florida,
        michigan = states

        @constraint(m, washington != oregon)
        @constraint(m, washington != idaho)
        @constraint(m, oregon != idaho)
        @constraint(m, oregon != nevada)
        @constraint(m, oregon != california)
        @constraint(m, california != nevada)
        @constraint(m, california != arizona)
        @constraint(m, nevada != idaho)
        @constraint(m, nevada != utah)
        @constraint(m, nevada != arizona)
        @constraint(m, idaho != montana)
        @constraint(m, idaho != wyoming)
        @constraint(m, idaho != utah)
        @constraint(m, utah != wyoming)
        @constraint(m, utah != colorado)
        @constraint(m, utah != new_mexico)
        @constraint(m, utah != arizona)
        @constraint(m, arizona != colorado)
        @constraint(m, arizona != new_mexico)
        @constraint(m, montana != north_dakota)
        @constraint(m, montana != south_dakota)
        @constraint(m, montana != wyoming)
        @constraint(m, wyoming != south_dakota)
        @constraint(m, wyoming != nebraska)
        @constraint(m, wyoming != colorado)
        @constraint(m, colorado != nebraska)
        @constraint(m, colorado != kansas)
        @constraint(m, colorado != oklahoma)
        @constraint(m, colorado != new_mexico)
        @constraint(m, new_mexico != oklahoma)
        @constraint(m, new_mexico != texas)
        @constraint(m, north_dakota != minnesota)
        @constraint(m, north_dakota != south_dakota)
        @constraint(m, south_dakota != minnesota)
        @constraint(m, south_dakota != iowa)
        @constraint(m, south_dakota != nebraska)
        @constraint(m, nebraska != iowa)
        @constraint(m, nebraska != missouri)
        @constraint(m, nebraska != kansas)
        @constraint(m, kansas != missouri)
        @constraint(m, kansas != oklahoma)
        @constraint(m, oklahoma != arkansas)
        @constraint(m, oklahoma != texas)
        @constraint(m, texas != arkansas)
        @constraint(m, texas != louisiana)
        @constraint(m, minnesota != wisconsin)
        @constraint(m, minnesota != iowa)
        @constraint(m, iowa != wisconsin)
        @constraint(m, iowa != illinois)
        @constraint(m, iowa != missouri)
        @constraint(m, missouri != illinois)
        @constraint(m, missouri != kentucky)
        @constraint(m, missouri != tennessee)
        @constraint(m, missouri != arkansas)
        @constraint(m, arkansas != tennessee)
        @constraint(m, arkansas != mississippi)
        @constraint(m, arkansas != louisiana)
        @constraint(m, louisiana != mississippi)
        @constraint(m, wisconsin != illinois)
        @constraint(m, wisconsin != michigan)
        @constraint(m, illinois != indiana)
        @constraint(m, illinois != kentucky)
        @constraint(m, kentucky != indiana)
        @constraint(m, kentucky != ohio)
        @constraint(m, kentucky != west_virginia)
        @constraint(m, kentucky != virginia)
        @constraint(m, kentucky != tennessee)
        @constraint(m, tennessee != virginia)
        @constraint(m, tennessee != north_carolina)
        @constraint(m, tennessee != georgia)
        @constraint(m, tennessee != alabama)
        @constraint(m, tennessee != mississippi)
        @constraint(m, mississippi != alabama)
        @constraint(m, michigan != indiana)
        @constraint(m, michigan != ohio)
        @constraint(m, indiana != ohio)
        @constraint(m, alabama != georgia)
        @constraint(m, alabama != florida)
        @constraint(m, ohio != pennsylvania)
        @constraint(m, ohio != west_virginia)
        @constraint(m, maine != new_hampshire)
        @constraint(m, new_hampshire != vermont)
        @constraint(m, new_hampshire != massachusetts)
        @constraint(m, vermont != massachusetts)
        @constraint(m, vermont != new_york)
        @constraint(m, massachusetts != new_york)
        @constraint(m, massachusetts != rhode_island)
        @constraint(m, massachusetts != connecticut)
        @constraint(m, connecticut != rhode_island)
        @constraint(m, new_york != connecticut)
        @constraint(m, new_york != pennsylvania)
        @constraint(m, new_york != new_jersey)
        @constraint(m, pennsylvania != new_jersey)
        @constraint(m, pennsylvania != delaware)
        @constraint(m, pennsylvania != maryland)
        @constraint(m, pennsylvania != west_virginia)
        @constraint(m, new_jersey != delaware)
        @constraint(m, maryland != washington_dc)
        @constraint(m, maryland != west_virginia)
        @constraint(m, maryland != virginia)
        @constraint(m, washington_dc != virginia)
        @constraint(m, west_virginia != virginia)
        @constraint(m, virginia != north_carolina)
        @constraint(m, north_carolina != south_carolina)
        @constraint(m, south_carolina != georgia)
        @constraint(m, georgia != florida)

        # test for AllEqual constraint
        @constraint(m, [california, new_york, florida] in CS.AllEqual())

        @constraint(m, max_color .>= states)

        # test for constant in objective
        @objective(m, Min, max_color + 1.1)

        optimize!(m)

        status = JuMP.termination_status(m)

        com = CS.get_inner_model(m)
        @test com.info.n_constraint_types.equality == 1
        @test com.info.n_constraint_types.inequality == length(states)

        CS.save_logs(com, "graph_color_optimize.json")
        rm("graph_color_optimize.json")

        @test status == MOI.OPTIMAL

        # all values fixed
        @test com.best_sol ≈ 5.1
        @test maximum([JuMP.value(var) for var in states]) == JuMP.value(max_color) == 4
        @test is_solved(com)
    end

    @testset "49 US states + DC only 3 colors" begin
        m = Model(CSJuMPTestOptimizer())
        num_colors = 3

        @variable(m, 1 <= states[1:50] <= num_colors, Int)
        washington,
        montana,
        maine,
        north_dakota,
        south_dakota,
        wyoming,
        wisconsin,
        idaho,
        vermont,
        minnesota,
        oregon,
        new_hampshire,
        iowa,
        massachusetts,
        nebraska,
        new_york,
        pennsylvania,
        connecticut,
        rhode_island,
        new_jersey,
        indiana,
        nevada,
        utah,
        california,
        ohio,
        illinois,
        washington_dc,
        delaware,
        west_virginia,
        maryland,
        colorado,
        kentucky,
        kansas,
        virginia,
        missouri,
        arizona,
        oklahoma,
        north_carolina,
        tennessee,
        texas,
        new_mexico,
        alabama,
        mississippi,
        georgia,
        south_carolina,
        arkansas,
        louisiana,
        florida,
        michigan = states

        @constraint(m, washington != oregon)
        @constraint(m, washington != idaho)
        @constraint(m, oregon != idaho)
        @constraint(m, oregon != nevada)
        @constraint(m, oregon != california)
        @constraint(m, california != nevada)
        @constraint(m, california != arizona)
        @constraint(m, nevada != idaho)
        @constraint(m, nevada != utah)
        @constraint(m, nevada != arizona)
        @constraint(m, idaho != montana)
        @constraint(m, idaho != wyoming)
        @constraint(m, idaho != utah)
        @constraint(m, utah != wyoming)
        @constraint(m, utah != colorado)
        @constraint(m, utah != new_mexico)
        @constraint(m, utah != arizona)
        @constraint(m, arizona != colorado)
        @constraint(m, arizona != new_mexico)
        @constraint(m, montana != north_dakota)
        @constraint(m, montana != south_dakota)
        @constraint(m, montana != wyoming)
        @constraint(m, wyoming != south_dakota)
        @constraint(m, wyoming != nebraska)
        @constraint(m, wyoming != colorado)
        @constraint(m, colorado != nebraska)
        @constraint(m, colorado != kansas)
        @constraint(m, colorado != oklahoma)
        @constraint(m, colorado != new_mexico)
        @constraint(m, new_mexico != oklahoma)
        @constraint(m, new_mexico != texas)
        @constraint(m, north_dakota != minnesota)
        @constraint(m, north_dakota != south_dakota)
        @constraint(m, south_dakota != minnesota)
        @constraint(m, south_dakota != iowa)
        @constraint(m, south_dakota != nebraska)
        @constraint(m, nebraska != iowa)
        @constraint(m, nebraska != missouri)
        @constraint(m, nebraska != kansas)
        @constraint(m, kansas != missouri)
        @constraint(m, kansas != oklahoma)
        @constraint(m, oklahoma != arkansas)
        @constraint(m, oklahoma != texas)
        @constraint(m, texas != arkansas)
        @constraint(m, texas != louisiana)
        @constraint(m, minnesota != wisconsin)
        @constraint(m, minnesota != iowa)
        @constraint(m, iowa != wisconsin)
        @constraint(m, iowa != illinois)
        @constraint(m, iowa != missouri)
        @constraint(m, missouri != illinois)
        @constraint(m, missouri != kentucky)
        @constraint(m, missouri != tennessee)
        @constraint(m, missouri != arkansas)
        @constraint(m, arkansas != tennessee)
        @constraint(m, arkansas != mississippi)
        @constraint(m, arkansas != louisiana)
        @constraint(m, louisiana != mississippi)
        @constraint(m, wisconsin != illinois)
        @constraint(m, wisconsin != michigan)
        @constraint(m, illinois != indiana)
        @constraint(m, illinois != kentucky)
        @constraint(m, kentucky != indiana)
        @constraint(m, kentucky != ohio)
        @constraint(m, kentucky != west_virginia)
        @constraint(m, kentucky != virginia)
        @constraint(m, kentucky != tennessee)
        @constraint(m, tennessee != virginia)
        @constraint(m, tennessee != north_carolina)
        @constraint(m, tennessee != georgia)
        @constraint(m, tennessee != alabama)
        @constraint(m, tennessee != mississippi)
        @constraint(m, mississippi != alabama)
        @constraint(m, michigan != indiana)
        @constraint(m, michigan != ohio)
        @constraint(m, indiana != ohio)
        @constraint(m, alabama != georgia)
        @constraint(m, alabama != florida)
        @constraint(m, ohio != pennsylvania)
        @constraint(m, ohio != west_virginia)
        @constraint(m, maine != new_hampshire)
        @constraint(m, new_hampshire != vermont)
        @constraint(m, new_hampshire != massachusetts)
        @constraint(m, vermont != massachusetts)
        @constraint(m, vermont != new_york)
        @constraint(m, massachusetts != new_york)
        @constraint(m, massachusetts != rhode_island)
        @constraint(m, massachusetts != connecticut)
        @constraint(m, connecticut != rhode_island)
        @constraint(m, new_york != connecticut)
        @constraint(m, new_york != pennsylvania)
        @constraint(m, new_york != new_jersey)
        @constraint(m, pennsylvania != new_jersey)
        @constraint(m, pennsylvania != delaware)
        @constraint(m, pennsylvania != maryland)
        @constraint(m, pennsylvania != west_virginia)
        @constraint(m, new_jersey != delaware)
        @constraint(m, maryland != washington_dc)
        @constraint(m, maryland != west_virginia)
        @constraint(m, maryland != virginia)
        @constraint(m, washington_dc != virginia)
        @constraint(m, west_virginia != virginia)
        @constraint(m, virginia != north_carolina)
        @constraint(m, north_carolina != south_carolina)
        @constraint(m, south_carolina != georgia)
        @constraint(m, georgia != florida)

        optimize!(m)

        @test JuMP.termination_status(m) == MOI.INFEASIBLE
    end

    @testset "Maximization objective" begin
        m = Model(CSJuMPTestOptimizer())
        num_colors = 20

        @variable(m, 1 <= states[1:50] <= num_colors, Int)
        @variable(m, 1 <= max_color <= num_colors, Int)

        washington,
        montana,
        maine,
        north_dakota,
        south_dakota,
        wyoming,
        wisconsin,
        idaho,
        vermont,
        minnesota,
        oregon,
        new_hampshire,
        iowa,
        massachusetts,
        nebraska,
        new_york,
        pennsylvania,
        connecticut,
        rhode_island,
        new_jersey,
        indiana,
        nevada,
        utah,
        california,
        ohio,
        illinois,
        washington_dc,
        delaware,
        west_virginia,
        maryland,
        colorado,
        kentucky,
        kansas,
        virginia,
        missouri,
        arizona,
        oklahoma,
        north_carolina,
        tennessee,
        texas,
        new_mexico,
        alabama,
        mississippi,
        georgia,
        south_carolina,
        arkansas,
        louisiana,
        florida,
        michigan = states

        @constraint(m, washington != oregon)
        @constraint(m, washington != idaho)
        @constraint(m, oregon != idaho)
        @constraint(m, oregon != nevada)
        @constraint(m, oregon != california)
        @constraint(m, california != nevada)
        @constraint(m, california != arizona)
        @constraint(m, nevada != idaho)
        @constraint(m, nevada != utah)
        @constraint(m, nevada != arizona)
        @constraint(m, idaho != montana)
        @constraint(m, idaho != wyoming)
        @constraint(m, idaho != utah)
        @constraint(m, utah != wyoming)
        @constraint(m, utah != colorado)
        @constraint(m, utah != new_mexico)
        @constraint(m, utah != arizona)
        @constraint(m, arizona != colorado)
        @constraint(m, arizona != new_mexico)
        @constraint(m, montana != north_dakota)
        @constraint(m, montana != south_dakota)
        @constraint(m, montana != wyoming)
        @constraint(m, wyoming != south_dakota)
        @constraint(m, wyoming != nebraska)
        @constraint(m, wyoming != colorado)
        @constraint(m, colorado != nebraska)
        @constraint(m, colorado != kansas)
        @constraint(m, colorado != oklahoma)
        @constraint(m, colorado != new_mexico)
        @constraint(m, new_mexico != oklahoma)
        @constraint(m, new_mexico != texas)
        @constraint(m, north_dakota != minnesota)
        @constraint(m, north_dakota != south_dakota)
        @constraint(m, south_dakota != minnesota)
        @constraint(m, south_dakota != iowa)
        @constraint(m, south_dakota != nebraska)
        @constraint(m, nebraska != iowa)
        @constraint(m, nebraska != missouri)
        @constraint(m, nebraska != kansas)
        @constraint(m, kansas != missouri)
        @constraint(m, kansas != oklahoma)
        @constraint(m, oklahoma != arkansas)
        @constraint(m, oklahoma != texas)
        @constraint(m, texas != arkansas)
        @constraint(m, texas != louisiana)
        @constraint(m, minnesota != wisconsin)
        @constraint(m, minnesota != iowa)
        @constraint(m, iowa != wisconsin)
        @constraint(m, iowa != illinois)
        @constraint(m, iowa != missouri)
        @constraint(m, missouri != illinois)
        @constraint(m, missouri != kentucky)
        @constraint(m, missouri != tennessee)
        @constraint(m, missouri != arkansas)
        @constraint(m, arkansas != tennessee)
        @constraint(m, arkansas != mississippi)
        @constraint(m, arkansas != louisiana)
        @constraint(m, louisiana != mississippi)
        @constraint(m, wisconsin != illinois)
        @constraint(m, wisconsin != michigan)
        @constraint(m, illinois != indiana)
        @constraint(m, illinois != kentucky)
        @constraint(m, kentucky != indiana)
        @constraint(m, kentucky != ohio)
        @constraint(m, kentucky != west_virginia)
        @constraint(m, kentucky != virginia)
        @constraint(m, kentucky != tennessee)
        @constraint(m, tennessee != virginia)
        @constraint(m, tennessee != north_carolina)
        @constraint(m, tennessee != georgia)
        @constraint(m, tennessee != alabama)
        @constraint(m, tennessee != mississippi)
        @constraint(m, mississippi != alabama)
        @constraint(m, michigan != indiana)
        @constraint(m, michigan != ohio)
        @constraint(m, indiana != ohio)
        @constraint(m, alabama != georgia)
        @constraint(m, alabama != florida)
        @constraint(m, ohio != pennsylvania)
        @constraint(m, ohio != west_virginia)
        @constraint(m, maine != new_hampshire)
        @constraint(m, new_hampshire != vermont)
        @constraint(m, new_hampshire != massachusetts)
        @constraint(m, vermont != massachusetts)
        @constraint(m, vermont != new_york)
        @constraint(m, massachusetts != new_york)
        @constraint(m, massachusetts != rhode_island)
        @constraint(m, massachusetts != connecticut)
        @constraint(m, connecticut != rhode_island)
        @constraint(m, new_york != connecticut)
        @constraint(m, new_york != pennsylvania)
        @constraint(m, new_york != new_jersey)
        @constraint(m, pennsylvania != new_jersey)
        @constraint(m, pennsylvania != delaware)
        @constraint(m, pennsylvania != maryland)
        @constraint(m, pennsylvania != west_virginia)
        @constraint(m, new_jersey != delaware)
        @constraint(m, maryland != washington_dc)
        @constraint(m, maryland != west_virginia)
        @constraint(m, maryland != virginia)
        @constraint(m, washington_dc != virginia)
        @constraint(m, west_virginia != virginia)
        @constraint(m, virginia != north_carolina)
        @constraint(m, north_carolina != south_carolina)
        @constraint(m, south_carolina != georgia)
        @constraint(m, georgia != florida)

        # test strictly less than
        @constraint(m, max_color .< states+1)

        @objective(m, Max, max_color)

        optimize!(m)

        status = JuMP.termination_status(m)

        com = CS.get_inner_model(m)

        @test status == MOI.OPTIMAL
        @test com.best_sol == 17
        @test minimum([JuMP.value(var) for var in states]) == 17 == JuMP.value(max_color)
    end

    @testset "small graph coloring correctness test" begin
        cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "logging" => [],
            "lp_optimizer" => cbc_optimizer,
            "keep_logs" => true,
        ))

        n = 10
        @variable(m, 1 <= c[1:n] <= n, Int)
        @variable(m, 1 <= max_color <= n, Int)

        edge_list = [
            (1, 2),
            (1, 5),
            (1, 7),
            (1, 8),
            (2, 4),
            (2, 7),
            (2, 8),
            (2, 9),
            (3, 7),
            (3, 9),
            (3, 10),
            (4, 9),
            (4, 10),
            (5, 10),
            (6, 7),
            (6, 9),
            (7, 8),
            (7, 10),
        ]

        for edge in edge_list
            @constraint(m, c[edge[1]] != c[edge[2]])
        end

        @constraint(m, max_color .>= c)
        @objective(m, Min, max_color)

        optimize!(m)

        status = JuMP.termination_status(m)
        objval = JuMP.objective_value(m)
        @test status == MOI.OPTIMAL
        @test objval ≈ 4
        colors = JuMP.value.(c)
        for edge in edge_list
            @test !(colors[edge[1]] ≈ colors[edge[2]])
        end
    end
end


text:
function gelman_rubin(samplers)
    # parameters in db must be scalars, vectors or matrices

    m = length(samplers)
    n = length(samplers[1])

    assert(all((map(length,samplers) .- n).==0)) # ... for now

    paramnames = names(samplers[1].db[1])
    dnostics = Dict()
    for p in paramnames
        p == :energy && continue
        ex = getfield(samplers[1].db[1],p)
        extype = typeof(ex)

        if extype <: Number
            # vector of vectors
            posts = [convert(Vector{extype}, map(y->getfield(y,p), x.db)) for x in samplers]
        elseif ndims(ex) == 1
            # vector of matrices
            posts = [vcat(map(y->getfield(y,p)', x.db)...) for x in samplers]
        elseif ndims(ex) == 2
            # vector of matrices
            posts = [vcat(map(y->vec(getfield(y,p))', x.db)...) for x in samplers]
        end

        # Dims: [chain] X [iteration X params]

        vars = vcat(map(x->var(x,1), posts)...)
        means = vcat(map(x->mean(x,1), posts)...)
        # Dims: [chain X params]
        W = mean(vars,1)
        B_jk = var(means,1)
        # Dims: [params]
        dnostics[string(p)] = vec(sqrt(((n-1)/n .* W .+ B_jk) ./ W))
    end
    return dnostics
end

function logsum(x::AbstractArray)
    if length(x) == 1
        return x[1]
    else
        a = maximum(x)
        return a .+ log(sum_kbn(exp(x.-a)))
    end
end

function logsum(x::AbstractArray, w::AbstractArray)
    if length(x) == 1 && w[1] == 1.0
        return x[1]
    else
        a = maximum(x)
        return a .+ log(sum_kbn(w.*exp(x.-a)))
    end
end

function logsum(x::Float64, y::Float64)
    x > y ? x+log1p(exp(y-x)) : y+log1p(exp(x-y))
end

function get_bbox(data; factor=1.5, positive=true)
    @show D = size(data,2)

    #maxs = vec(maximum(data, 1))
    maxs = vec(mapslices(x->quantile(x,0.75), data, 2))
    maxs .+= 3*maximum(maxs) / 4

    if positive
        mins = zeros(D)
    else
        #mins = vec(minimum(data,1))
        maxs = vec(mapslices(x->quantile(x,0.25), data, 2))
    end

    spreads = maxs .- mins

    if !positive
        mins .-= spreads .* factor
    end
    maxs .+= spreads .* factor
    mins, maxs
end

function get_grid(data; factor=1.5, N=30, positive=true)
    mins, maxs = get_bbox(data, factor=factor, positive=positive)
    mins, maxs, gen_grid(mins, maxs, N)
end

gen_unit_grid(mins,maxs) = gen_grid(mins,maxs,typemax(Int))

function gen_grid(mins, maxs, N=30)
    D = length(mins)
    stepsizes = ceil(float(maxs .- mins) ./N)
    ranges = [mins[i]:stepsizes[i]:maxs[i] for i=1:D]
    grid = hcat(map(collect, product(ranges...))...)
    return ranges, grid
end


text:
#this file holds functions specific to acquiring the NCCS data

function getNCCSFile(url::String, fileName::String, outPath::String;
  outName::String = fileName, validate::Bool=true, update::Bool=false)::Nothing

  if update && isfile("$outPath\\$outName.gz")# if the file exists
    if validate #validate if desired
      try #validate the file
        validateSingleTypeCSV("$outPath\\$outName.gz")
        println("Download Validated: $outName")
        return nothing
      end
    else
      return nothing
    end
  end

  #download the file
  #oStream::GZipStream = GZip.open("$outPath\\$outName.gz", "w")
  oStream::GZipStream = GZip.open("$outPath\\$outName.gz", "w")
  #weird error if using the transcode codec

  #acquire the data
  HTTP.open(io::IO->
    while !eof(io)
      write(oStream, readavailable(io))
    end,
    "GET", "$url/$fileName", retries=4, connection_limit=16,
    pipeline_limit=32, readtimeout=60)
  GZip.close(oStream)

  #check if its a valid file
  if validate
    try
      validateSingleTypeCSV("$outPath\\$outName.gz")
      println("Download Validated: $outName")
    catch err
      warn("DL Error with $outName. \nErrorMessage:\n$err")
    end
  end

  return nothing
end

#vectorized version (useful for parallelization)
function getNCCSFiles(urls::Vector{String}, fileNames::Vector{String},
  outPath::String, mode::Symbol = :serial;
  outNames::Vector{String} = fileNames, delay::Float64 = 0.0,
  validate::Bool = true, update::Bool=false)::Nothing

  #short-circuit for parallel operation
  if mode == :parallel
    return getNCCSFilesPar(urls, fileNames, outPath, outNames=outNames,
    delay=delay, validate=validate, update=update)
  end

  for i::Int ∈ 1:length(urls)

    if delay ≥ .001 #allows for a delay to avoid slamming the site
      sleep(delay)
    end

    getNCCSFile(urls[i], fileNames[i], outPath, validate=validate, update=update)
  end

  return nothing
end

function getNCCSFilesPar(urls::Vector{String}, fileNames::Vector{String},
  outPath::String; outNames::Vector{String} = fileNames,
  delay::Float64 = 0.0, validate::Bool = true, update::Bool=false)::Nothing

  numFiles::Int = length(urls)
  pids = workers()
  np::Int = length(pids)

  #divide up the work
  assignments::Vector{Vector{Int}} = ((i::Int)->Vector{Int}()).(1:np)
  for i::Int ∈ 1:numFiles
    push!(assignments[i%np+1], i) #assigns a file to a pid
  end

  #set up a vector for pushing and receiving assignments
  futures::Vector{Future} = Vector{Future}()
  sizehint!(futures, np)

  for i::Int ∈ 1:np
    if length(assignments[i]) > 0
      push!(futures,
        (@spawnat pids[i] getNCCSFiles(urls[assignments[i]],
          fileNames[assignments[i]], outPath,
          outNames=outNames, delay=delay, validate=validate, update=update)))
    end
  end

  #get the work
  (fetch).(futures)

  return nothing
end

function getNCCSFilesIndex(;inPath::String = DATA_PATH, filesIndexName::String = FILES_INDEX_NAME)::DataFrame
  return CSV.read("$inPath\\$filesIndexName.csv", strings=:raw,
    categorical=false)
end

function getNCCSFiles(outPath::String; mode::Symbol=:parallel,
  delay::Float64=0.0, validate::Bool = true, update::Bool=false)

  #read in the metadata
  fileDF::DataFrame = getNCCSFilesIndex()
  numFiles::Int = size(fileDF,1)

  #build the urls
  urls::Vector{String} = ((urlStem::String, urlYear::Int)->
    "$urlStem/$urlYear").(fileDF[:urlStems], fileDF[:urlYears])
  fileNames::Vector{String} =  fileDF[:fileNames]

  return getNCCSFiles(urls, fileNames, outPath, mode, delay=delay,
    validate=validate, update=update)

end

#serializes and compresses an already downloaded NCCS file
function serializeNCCSFile(meta::MetaDictionary; nccsPath::String = NCCS_PATH,
  iStreamer::Function = I_CSV_STREAMER, oStreamer::Function = O_JLS_STREAMER,
  serializeOnly::Bool = true, verbose::Bool = true)::DataFrame

  #these are the column names that we expect
  expectedCols::Vector{Symbol} = collect(keys(meta.dictionary))

  inPath::String = "$nccsPath\\$(meta.dataFileName).csv.gz"
  outPath::String = "$nccsPath\\$(meta.dataFileName).jls.gz"

  if verbose
    println("Serializing $(meta.dataFileName)...")
  end

  iStream = iStreamer(inPath)
  df::DataFrame = CSV.read(iStream, strings=:raw, categorical=false, rows=2)
  close(iStream)

  #get the column names and make sure each one is in the list of columns we want
  actualCols::Vector{String} = (s::Symbol->lowercase(string(s))).(names(df))

  #form the lookup dictionary for the type. Takes a column name and gets the type, with a default of String
  typeDict::Dict = Dict(
    actualCols[i] => (
    haskey(meta.dictionary, actualCols[i]) ? meta.dictionary[Symbol(actualCols[i])]: Union{String, Missing})
    for i::Int ∈ 1:(length(actualCols)))

  #now read in the file
  iStream = iStreamer(inPath)
  #CSV.validate(iStream, header=actualCols, datarow=2, types=typeDict)
  df = CSV.read(iStream, header=actualCols, datarow=2, types=typeDict)
  close(iStream)

  #remove columns we don't expect or want
  try
    df = df[:,expectedCols]
  catch err
    error("$err\nmeta: $(meta)")
  end

  oStream = oStreamer(outPath)
  serialize(oStream, df)
  close(oStream)

  if verbose
    println("Serialization of $(meta.dataFileName) complete.")
  end

  if serializeOnly
    df = DataFrame([[nothing]]) #free the memory with an empty dataframe
  end

  return df

end

#loads an already downloaded NCCS file
function loadNCCSFile(meta::MetaDictionary; nccsPath::String = NCCS_PATH,
  refresh::Bool = true, iStreamer::Function = I_JLS_STREAMER, serializeOnly::Bool = true)

  filePath::String = "$nccsPath\\$(meta.dataFileName).jls.gz"

  #only re-process the file if necessesary
  if refresh || !isfile(filePath)
    df::DataFrame = serializeNCCSFile(meta, serializeOnly=serializeOnly)
  elseif !serializeOnly
    iStream = iStreamer(filePath)
    df = deserialize(iStream)
    close(iStream)
  else #do nothing
    df = DataFrame([[nothing]]) #an empty dataframe
  end

  return df
end

#loads and/or serializes a vector of corresponding files
function loadNCCSFiles(metadata::Vector{MetaDictionary}; mode::Symbol = :serial,
  refresh::Bool = true, returnFiles::Bool = false, serializeOnly::Bool = true)::Vector{DataFrame}

  if mode == :parallel #relaunch the dedicated parallel function
    return loadNCCSFilesPar(metadata::Vector{MetaDictionary}, refresh=refresh, serializeOnly=serializeOnly)
  end

  if serializeOnly
    (meta::MetaDictionary->loadNCCSFile(meta,refresh=refresh, serializeOnly=serializeOnly)).(metadata)
    return Vector{DataFrame}()
  else
    return (meta::MetaDictionary->loadNCCSFile(meta,refresh=refresh)).(metadata)
  end
end

#loads many NCCS serial files in parallel
function loadNCCSFilesPar(metadata::Vector{MetaDictionary}; refresh::Bool = true,
  serializeOnly::Bool=true)::Vector{DataFrame}

  numFiles::Int = length(metadata)
  pids = workers()
  np::Int = length(pids)

  #divide up the work
  assignments::Vector{Vector{Int}} = ((i::Int)->Vector{Int}()).(1:np)
  for i::Int ∈ 1:numFiles
    push!(assignments[i%np+1], i) #assigns a file to a pid
  end

  #set up a vector for pushing and receiving assignments
  futures::Vector{Future} = Vector{Future}()
  sizehint!(futures, np)

  for i::Int ∈ 1:np
    if length(assignments[i]) > 0
      push!(futures,
        (@spawnat pids[i] loadNCCSFiles(metadata[assignments[i]], refresh=refresh, serializeOnly=serializeOnly)))
    end
  end

  #get the work
  return vcat((fetch).(futures)...)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# common facilities

typealias DenseVector{T} DenseArray{T,1}
typealias DenseMatrix{T} DenseArray{T,2}
typealias DenseVecOrMat{T} Union(DenseVector{T},DenseMatrix{T})

typealias NumericArray{T<:Number,N} DenseArray{T,N}
typealias DenseArrOrNum{T<:Number} Union(DenseArray{T}, T)
typealias ArrayOrNum{T<:Number} Union(Array{T}, T)

typealias NumericVector{T<:Number} DenseVector{T}
typealias NumericMatrix{T<:Number} DenseMatrix{T}

typealias ContiguousArrOrNum{T<:Number} Union(ContiguousArray{T}, T)
typealias ContiguousNumericArray{T<:Number} ContiguousArray{T}
typealias ContiguousRealArray{T<:Real} ContiguousArray{T}

typealias DimSpec Union(Int,Dims,Vector{Int})
typealias BlasFP Union(Float32,Float64)

parent(a::Number) = a
stride(a::Number, d::Int) = 1
ellipview(a::Number, i::Int) = a

getvalue(a::Number, i::Integer) = a
getvalue(a::DenseArray, i::Integer) = a[i]

gt_or_nan(s::Number, x) = (s > x)
lt_or_nan(s::Number, x) = (s < x)

gt_or_nan(s::FloatingPoint, x) = (s > x || s != s)
lt_or_nan(s::FloatingPoint, x) = (s < x || s != s)

type _Max <: Functor{2} end
evaluate{T<:Real}(::_Max, x::T, y::T) = ifelse(y > x, y, x)

type NonnegMax <: Functor{2} end
evaluate{T<:Real}(::NonnegMax, x::T, y::T) = ifelse(y > x, y, x)

type _Min <: Functor{2} end
evaluate{T<:Real}(::_Min, x::T, y::T) = ifelse(y < x, y, x)

result_type{T<:Real}(::Union(_Max,_Min,NonnegMax), ::Type{T}, ::Type{T}) = T


text:
nanmean(x::AbstractVector)::Float64 = mean(map(n -> isnan(n) ? 0 : n, x))
nanmean(x::AbstractMatrix)::Vector{Float64} = map(nanmean, eachrow(x))
absnanmean(x::AbstractVector)::Float64 = nanmean(map(abs, x))
absnanmean(x::AbstractMatrix)::Vector{Float64} = map(absnanmean, eachrow(x))

# The code for `calc_drop` and `drop_cor` below is inspired by the code by Brian Pietracatella, found at
# https://towardsdatascience.com/are-you-dropping-too-many-correlated-features-d1c96654abe6

function calc_drop(res::DataFrame)::Vector{Symbol}
    # Get all variables that have some high correlation
    all_corr_vars = unique([res[!, :v1]; res[!, :v2]])
    # Get all variables that are possibly dropped
    poss_drop = unique(res[!, :drop])
    # Get variables that are definitely kept
    keep = setdiff(all_corr_vars, poss_drop)
    # Drop variables that are highly correlated to a variable that is definitely kept
    p = res[in(keep).(res[!, :v1]) .| in(keep).(res[!, :v2]), [:v1, :v2]]
    q = unique([p[!, :v1]; p[!, :v2]])
    drop = setdiff(q, keep)
    # Remove dropped variables from the list of possibilities
    poss_drop = setdiff(poss_drop, drop)
    # Get variables that (correlate highly to variables that) are still possibly dropped
    m = res[in(poss_drop).(res[!, :v1]) .| in(poss_drop).(res[!, :v2]), [:v1, :v2, :drop]]
    # Drop variables that are not (correlated to variables that are) already dropped
    more_drop = m[.!in(drop).(m[!, :v1]) .& .!in(drop).(m[!, :v2]), :drop]
    return [drop; more_drop]
end

function drop_cor(cormat::Matrix{Float64}, xs::Vector{Symbol}, t::Float64)::Vector{Symbol}
    msize = size(cormat, 1)
    (msize == size(cormat, 2)) || throw("cormat is not a square matrix")
    (msize == length(xs)) || throw("length of xs does not match size of cormat")

    avgcor = absnanmean(cormat)

    res = DataFrame(v1 = Symbol[], v2 = Symbol[], v1mean = Float64[], v2mean = Float64[], corr = Float64[], drop = Symbol[])

    @inbounds for row in 1:(msize - 1)
        for col in (row + 1):msize
            if abs(cormat[row, col]) > t
                rowsym = xs[row]
                colsym = xs[col]
                drop = avgcor[row] > avgcor[col] ? rowsym : colsym
                entry = [
                    rowsym, colsym,
                    avgcor[row], avgcor[col],
                    cormat[row, col], drop
                ]
                push!(res, entry)
            end
        end
    end

    return calc_drop(res)

end

max_stars_repo_name: KasperNooteboom/thesis-rvfl-fs

max_issues_repo_name: KasperNooteboom/thesis-rvfl-fs

max_forks_repo_name: KasperNooteboom/thesis-rvfl-fs


text:
#
# Datasets
#

# Use an Enum here when it doesn't make us allocate
const LC_COMPACT_STORAGE = 0x00
const LC_CONTIGUOUS_STORAGE = 0x01
const LC_CHUNKED_STORAGE = 0x02

function read_dataset(f::JLDFile, offset::RelOffset)
    if haskey(f.jloffset, offset)
        # Stored as WeakRefs and may no longer exist
        val = f.jloffset[offset].value
        val !== nothing && return val
    end

    io = f.io
    seek(io, fileoffset(f, offset))
    cio = begin_checksum(io)
    sz = read_obj_start(cio)
    pmax = position(cio) + sz

    # Messages
    dataspace = ReadDataspace()
    attrs = EMPTY_READ_ATTRIBUTES
    datatype_class::UInt8 = 0
    datatype_offset::Int64 = 0
    data_offset::Int64 = 0
    data_length::Int = 0
    while position(cio) < pmax
        msg = read(cio, HeaderMessage)
        endpos = position(cio) + msg.size
        if msg.msg_type == HM_DATASPACE
            dataspace = read_dataspace_message(io)
        elseif msg.msg_type == HM_DATATYPE
            datatype_class, datatype_offset = read_datatype_message(io, f, (msg.flags & 2) == 2)
        elseif msg.msg_type == HM_FILL_VALUE
            (read(cio, UInt8) == 3 && read(cio, UInt8) == 0x09) || throw(UnsupportedFeatureException())
        elseif msg.msg_type == HM_DATA_LAYOUT
            read(cio, UInt8) == 3 || throw(UnsupportedVersionException())
            storage_type = read(cio, UInt8)
            if storage_type == LC_COMPACT_STORAGE
                data_length = read(cio, UInt16)
                data_offset = position(cio)
            elseif storage_type == LC_CONTIGUOUS_STORAGE
                data_offset = fileoffset(f, read(cio, RelOffset))
                data_length = read(cio, Length)
            else
                throw(UnsupportedFeatureException())
            end
        elseif msg.msg_type == HM_ATTRIBUTE
            if attrs === EMPTY_READ_ATTRIBUTES
                attrs = ReadAttribute[read_attribute(cio, f)]
            else
                push!(attrs, read_attribute(cio, f))
            end
        elseif (msg.flags & 2^3) != 0
            throw(UnsupportedFeatureException())
        elseif msg.msg_type == HM_NIL
            break
        end
        seek(cio, endpos)
    end
    seek(cio, pmax)

    # Checksum
    end_checksum(cio) == read(io, UInt32) || throw(InvalidDataException())

    # TODO verify that data length matches
    val = read_data(f, dataspace, datatype_class, datatype_offset, data_offset,
                    offset, attrs)
    val
end

# Read data from an attribute
read_attr_data(f::JLDFile, attr::ReadAttribute) =
    read_data(f, attr.dataspace, attr.datatype_class, attr.datatype_offset,
              attr.data_offset)

# Read data from an attribute, assuming a specific HDF5 datatype and
# ReadRepresentation. If the HDF5 datatype does not match, throws an
# UnsupportedFeatureException. This allows better type stability while
# simultaneously validating the data.
function read_attr_data(f::JLDFile, attr::ReadAttribute, expected_datatype::H5Datatype, rr::ReadRepresentation)
    io = f.io
    if attr.datatype_class == class(expected_datatype)
        seek(io, attr.datatype_offset)
        dt = read(io, typeof(expected_datatype))
        if dt == expected_datatype
            seek(f.io, attr.data_offset)
            BOXED_READ_DATASPACE[] = (attr.dataspace, NULL_REFERENCE)
            return read_data(f, rr)
        end
    end
    throw(UnsupportedFeatureException())
end

# Read data from a file. If datatype_class is typemax(UInt8), the
# datatype is assumed to be committed, and datatype_offset points to
# the offset of the committed datatype's header. Otherwise,
# datatype_offset points to the offset of the datatype attribute.
function read_data(f::JLDFile, dataspace::ReadDataspace,
                   datatype_class::UInt8, datatype_offset::Int64,
                   data_offset::Int64, header_offset::RelOffset=NULL_REFERENCE,
                   attributes::Union{Vector{ReadAttribute},Void}=nothing)
    # See if there is a julia type attribute
    io = f.io
    if datatype_class == typemax(UInt8) # Committed datatype
        rr = jltype(f, f.datatype_locations[h5offset(f, datatype_offset)])
        seek(io, data_offset)
        BOXED_READ_DATASPACE[] = (dataspace, header_offset)
        read_data(f, rr, attributes)
    else
        seek(io, datatype_offset)
        @read_datatype io datatype_class dt begin
            rr = jltype(f, dt)
            seek(io, data_offset)
            BOXED_READ_DATASPACE[] = (dataspace, header_offset)
            read_data(f, rr, attributes)
        end
    end
end

# We can avoid a box by putting the dataspace here instead of passing
# it to read_data.  That reduces the allocation footprint, but doesn't
# really seem to help with performance.
const BOXED_READ_DATASPACE = Ref{Tuple{ReadDataspace,RelOffset}}()

# Most types can only be scalars or arrays
function read_data(f::JLDFile{MmapIO}, rr,
                   attributes::Union{Vector{ReadAttribute},Void}=nothing)
    dataspace, header_offset = BOXED_READ_DATASPACE[]
    if dataspace.dataspace_type == DS_SCALAR
        io = f.io
        inptr = io.curptr
        obj = jlconvert(rr, f, inptr, header_offset)
        io.curptr = inptr + sizeof(rr)
        obj
    elseif dataspace.dataspace_type == DS_SIMPLE
        read_array(f, dataspace, rr, attributes, header_offset)
    else
        throw(UnsupportedFeatureException())
    end
end

# Reference arrays can only be arrays or null dataspace (for Union{} case)
function read_data(f::JLDFile{MmapIO}, rr::ReadRepresentation{Any,RelOffset},
                   attributes::Vector{ReadAttribute})
    dataspace, header_offset = BOXED_READ_DATASPACE[]
    if dataspace.dataspace_type == DS_SIMPLE
        # Since this is an array of references, there should be an attribute
        # informing us of the type
        io = f.io
        startpos = position(io)
        for x in attributes
            if x.name == :julia_type
                T = read_attr_data(f, x)
                if isa(T, UnknownType)
                    str = typestring(T)
                    warn("type $(str) does not exist in workspace; interpreting Array{$str} as Array{Any}")
                    T = Any
                end
                seek(io, startpos)
                return read_array(f, dataspace, ReadRepresentation{T,RelOffset}(),
                                  attributes, header_offset)
            end
        end
    elseif dataspace.dataspace_type == DS_NULL
        return read_empty(ReadRepresentation{Union{},nothing}(), f,
                          attributes[find_dimensions_attr(attributes)],
                          header_offset)
    end
    throw(UnsupportedFeatureException())
end

# Types with no payload can only be null dataspace
function read_data{T}(f::JLDFile{MmapIO}, rr::ReadRepresentation{T,nothing},
                      attributes::Vector{ReadAttribute})
    dataspace, header_offset = BOXED_READ_DATASPACE[]
    dataspace.dataspace_type == DS_NULL || throw(UnsupportedFeatureException())

    dimensions_attr_index = find_dimensions_attr(attributes)
    if dimensions_attr_index == 0
        jlconvert(rr, f, Ptr{Void}(0), header_offset)
    else
        # dimensions attribute => array of empty type
        read_empty(rr, f, attributes[dimensions_attr_index], header_offset)
    end
end

function read_data{T,S}(f::JLDFile{MmapIO}, rr::ReadRepresentation{T,CustomSerialization{S,nothing}},
                      attributes::Vector{ReadAttribute})
    dataspace, header_offset = BOXED_READ_DATASPACE[]
    dataspace.dataspace_type == DS_NULL || throw(UnsupportedFeatureException())

    dimensions_attr_index = find_dimensions_attr(attributes)
    if dimensions_attr_index == 0
        jlconvert(rr, f, Ptr{Void}(0), header_offset)
    else
        # dimensions attribute => array of empty type
        read_empty(rr, f, attributes[dimensions_attr_index], header_offset)
    end
end

function find_dimensions_attr(attributes::Vector{ReadAttribute})
    dimensions_attr_index = 0
    julia_type_attr_index = 0
    for i = 1:length(attributes)
        x = attributes[i]
        if x.name == :dimensions
            dimensions_attr_index = i
        end
    end
    dimensions_attr_index
end

function read_empty{T}(rr::ReadRepresentation{T}, f::JLDFile{MmapIO},
                    dimensions_attr::ReadAttribute, header_offset::RelOffset)
    io = f.io
    inptr = io.curptr
    seek(io, dimensions_attr.dataspace.dimensions_offset)
    ndims = read(io, Int)
    seek(io, dimensions_attr.data_offset)
    v = construct_array(io, T, ndims)
    if isleaftype(T)
        for i = 1:length(v)
            @inbounds v[i] = jlconvert(rr, f, Ptr{Void}(0), header_offset)
        end
    end
    header_offset !== NULL_REFERENCE && (f.jloffset[header_offset] = WeakRef(v))
    io.curptr = inptr
    v
end

get_ndims_offset(f::JLDFile, dataspace::ReadDataspace, attributes::Void) =
    (dataspace.dimensionality, dataspace.dimensions_offset)

function get_ndims_offset(f::JLDFile, dataspace::ReadDataspace, attributes::Vector{ReadAttribute})
    ndims = dataspace.dimensionality
    offset = dataspace.dimensions_offset
    if !isempty(attributes)
        for x in attributes
            if x.name == :dimensions
                (x.dataspace.dataspace_type == DS_SIMPLE &&
                 x.dataspace.dimensionality == 1) || throw(InvalidDataException())
                seek(f.io, x.dataspace.dimensions_offset)
                ndims = UInt8(read(f.io, Length))
                offset = x.data_offset
            end
        end
    end
    (ndims, offset)
end

# Construct array by reading ndims dimensions from io
# Assumes io has already been seeked t oteh correct position
function construct_array{T}(io::IO, ::Type{T}, ndims::Int)
    if ndims == 1
        n = read(io, Int)
        Vector{T}(n)
    elseif ndims == 2
        d2 = read(io, Int)
        d1 = read(io, Int)
        Matrix{T}(d1, d2)
    elseif ndims == 3
        d3 = read(io, Int)
        d2 = read(io, Int)
        d1 = read(io, Int)
        Array{T,3}(d1, d2, d3)
    else
        ds = reverse!(read(io, Int, ndims))
        Array{T}(tuple(ds...))
    end
end

function read_array{T,RR}(f::JLDFile, dataspace::ReadDataspace,
                          rr::ReadRepresentation{T,RR},
                          attributes::Union{Vector{ReadAttribute},Void},
                          header_offset::RelOffset)
    io = f.io
    inptr = io.curptr
    ndims, offset = get_ndims_offset(f, dataspace, attributes)
    seek(io, offset)
    v = construct_array(io, T, Int(ndims))
    header_offset !== NULL_REFERENCE && (f.jloffset[header_offset] = WeakRef(v))
    n = length(v)

    if RR === T && isbits(T)
        nb = n*sizeof(RR)
        if nb > 1048576 # TODO tweak
            # It turns out that regular IO is faster here (at least on OS X)
            mmapio = f.io
            regulario = mmapio.f
            seek(regulario, inptr - pointer(f.io.arr))
            if ccall(:ios_readall, UInt,
                     (Ptr{Void}, Ptr{Void}, UInt), regulario.ios, v, nb) < nb
                throw(EOFError())
            end
        else
            unsafe_copy!(pointer(v), convert(Ptr{T}, inptr), Int(n))
        end
    # Would this actually help with performance?
    # elseif isbits(T)
    #     @simd for i = 1:n
    #         jlconvert!(outptr, rr, f, inptr)
    #         inptr += sizeof(RR)
    #         outptr += sizeof(T)
    #     end
    else
        @simd for i = 1:n
            if !jlconvert_canbeuninitialized(rr) || jlconvert_isinitialized(rr, inptr)
                @inbounds v[i] = jlconvert(rr, f, inptr, NULL_REFERENCE)
            end
            inptr += sizeof(RR)
        end
    end

    io.curptr = inptr + sizeof(RR) * n
    v
end

function payload_size(dataspace::WriteDataspace, datatype::H5Datatype, datasz::Int, layout_class::UInt8)
    sz = sizeof(dataspace) + sizeof(datatype) + 2 + (4 + length(dataspace.attributes))*sizeof(HeaderMessage) + 2
    for attr in dataspace.attributes
        sz += sizeof(attr)
    end
    if layout_class == LC_COMPACT_STORAGE
        sz + 2 + datasz
    else
        sz + sizeof(RelOffset) + sizeof(Length)
    end
end

# Might need to do something else someday for non-mmapped IO
function write_data{S}(f::JLDFile, data, odr::S, wsession::JLDWriteSession)
    io = f.io
    ensureroom(io, sizeof(odr))
    arr = io.arr
    cp = io.curptr
    h5convert!(cp, odr, f, data, wsession)
    io.curptr = cp + sizeof(odr)
    arr # Keep old array rooted until the end
end

# Like isdefined, but assumes arr is a pointer array, and can be inlined
unsafe_isdefined(arr::Array, i::Int) =
    unsafe_load(Ptr{Ptr{Void}}(pointer(arr)+(i-1)*sizeof(Ptr{Void}))) != Ptr{Void}(0)

function write_data{T,S}(f::JLDFile, data::Array{T}, odr::S, wsession::JLDWriteSession)
    io = f.io
    ensureroom(io, sizeof(odr) * length(data))
    arr = io.arr
    cp = io.curptr
    ep = io.endptr
    @simd for i = 1:length(data)
        if (isleaftype(T) && isbits(T)) || unsafe_isdefined(data, i)
            # For now, just don't write anything unless the field is defined
            @inbounds h5convert!(cp, odr, f, data[i], wsession)
        else
            @inbounds h5convert_uninitialized!(cp, odr)
        end
        cp += sizeof(odr)
    end
    io.curptr = cp
    arr # Keep old array rooted until the end
end

function write_dataset{S}(f::JLDFile, dataspace::WriteDataspace, datatype::H5Datatype, odr::S, data, wsession::JLDWriteSession)
    io = f.io
    datasz = sizeof(odr) * numel(dataspace)
    layout_class = datasz < 8192 ? LC_COMPACT_STORAGE : LC_CONTIGUOUS_STORAGE
    psz = payload_size(dataspace, datatype, datasz, layout_class)
    fullsz = sizeof(ObjectStart) + size_size(psz) + psz + 4

    header_offset = f.end_of_data
    seek(io, header_offset)
    f.end_of_data = header_offset + fullsz + (layout_class == LC_CONTIGUOUS_STORAGE ? datasz : 0)

    if typeof(data).mutable && !isa(wsession, JLDWriteSession{Union{}})
        wsession.h5offset[object_id(data)] = h5offset(f, header_offset)
        push!(wsession.objects, data)
    end

    cio = begin_checksum(io, fullsz - 4)

    write(cio, ObjectStart(size_flag(psz)))
    write_size(cio, psz)

    # Dataspace
    write(cio, HeaderMessage(HM_DATASPACE, sizeof(dataspace), 0))
    write(cio, dataspace)

    # Datatype
    write(cio, HeaderMessage(HM_DATATYPE, sizeof(datatype), 1+2*isa(datatype, CommittedDatatype)))
    write(cio, datatype)

    # Fill value
    write(cio, HeaderMessage(HM_FILL_VALUE, 2, 0))
    write(cio, UInt8(3)) # Version
    write(cio, 0x09)     # Flags

    # Attributes
    for attr in dataspace.attributes
        write(cio, HeaderMessage(HM_ATTRIBUTE, sizeof(attr), 0))
        write_attribute(cio, f, attr, f.datatype_wsession)
    end

    # Data storage layout
    if layout_class == LC_COMPACT_STORAGE
        write(cio, HeaderMessage(HM_DATA_LAYOUT, 4+datasz, 0))
        write(cio, UInt8(3))                  # Version
        write(cio, LC_COMPACT_STORAGE)        # Layout class
        write(cio, UInt16(datasz))            # Size
        if datasz != 0
            write_data(f, data, odr, wsession)
            seek(io, header_offset + fullsz - 4)
        end
        write(io, end_checksum(cio))
    else
        write(cio, HeaderMessage(HM_DATA_LAYOUT, 2+sizeof(RelOffset)+sizeof(Length), 0))
        write(cio, UInt8(3))                            # Version
        write(cio, LC_CONTIGUOUS_STORAGE)               # Layout class
        write(cio, h5offset(f, header_offset + fullsz)) # RelOffset
        write(cio, Length(sizeof(data)))                # Length
        write(io, end_checksum(cio))
        if datasz != 0
            write_data(f, data, odr, wsession)
        end
    end

    h5offset(f, header_offset)
end

@Base.pure ismutabletype(x::DataType) = x.mutable

@inline function write_dataset(f::JLDFile, x, wsession::JLDWriteSession)
    odr = objodr(x)
    write_dataset(f, WriteDataspace(f, x, odr), h5type(f, x), odr, x, wsession)
end

@inline function write_ref_mutable(f::JLDFile, x, wsession::JLDWriteSession)
    offset = get(wsession.h5offset, object_id(x), RelOffset(0))
    offset != RelOffset(0) ? offset : write_dataset(f, x, wsession)::RelOffset
end

@inline write_ref_mutable(f::JLDFile, x, wsession::JLDWriteSession{Union{}}) =
    write_dataset(f, x, wsession)

function write_ref(f::JLDFile, x, wsession::JLDWriteSession)
    if ismutabletype(typeof(x))
        write_ref_mutable(f, x, wsession)::RelOffset
    else
        write_dataset(f, x, wsession)::RelOffset
    end
end
write_ref(f::JLDFile, x::RelOffset, wsession::JLDWriteSession) = x


text:
using DelimitedFiles
using Plots; gr()

function plot_iterations_linear(
                                a::Vector{Float64},
                                b::Vector{Float64},
                                filename::String,
                                imagename::String;
                                p_step::Float64=1.0e-2,
                                p_levels::Int64=30,
                                show_sample_set::Bool=false
                                )

    # Reads the execution data
    data = readdlm(filename)

    # Splits the useful data
    data_c = data[1:(end - 1), 14]
    data_g = data[1:(end - 1), 15:16]
    data_xbase = data[1:(end - 1), 17:18]
    data_xopt = data[1:(end - 1), 19:20]

    if show_sample_set

        data_Y = data[1:(end - 1), 24:27]

    end

    # Creates 'x' and 'y'
    x = a[1]:p_step:b[1]
    y = a[2]:p_step:b[2]

    # Generic model function
    function model_i(i, x, y)

        return data_c[i] + data_g[i, 1] * x + data_g[i, 2] * y

    end

    for i=1:length(data_c)

        # i-th iteration model

        model(x, y) = model_i(i, x, y)
        plot_model = contour( x, y, model, levels = p_levels )

        if show_sample_set

            Y_x = [ data_Y[i, 1] + data_xbase[i, 1], data_Y[i, 3] + data_xbase[i, 1] ]
            Y_y = [ data_Y[i, 2] + data_xbase[i, 2], data_Y[i, 4] + data_xbase[i, 2] ]

            scatter!( plot_model, Y_x, Y_y, color = :gray40, markersize = 6, marker = :circle,  lab = "sample point" )

        end

        scatter!( plot_model, [ data_xbase[i, 1] ], [ data_xbase[i, 2] ], color = :green3, markersize = 12, marker = :circle,  lab = "xbase" )
        scatter!( plot_model, [ data_xopt[i, 1] ] , [ data_xopt[i, 2] ], color = :red  , markersize = 9,  marker = :diamond, lab = "xopt" )

        plot( plot_model, size = (900, 600))
        figname = string( imagename, "$(i)" )
        png( figname )

    end

end

# -----------------------------------------------
# Problem definition
# -----------------------------------------------

a = [0.0, 0.0]
b = [5.0, 5.0]
path_to_source = "./data_files/simple_runtest_01.dat"
path_to_image = "./images/test_01/fig_model_"

# -----------------------------------------------
# Function call
# -----------------------------------------------

plot_iterations_linear( a, b, path_to_source, path_to_image; show_sample_set=true)


text:
using AlphaZero
using ProgressMeter

# include("../games/tictactoe/main.jl")
# using .Tictactoe: Game

include("../games/connect-four/main.jl")
using .ConnectFour: Game

baseline = MinMax.Player{Game}(depth=5, amplify_rewards=true, τ=0.2)
mcts = MCTS.Env{Game}(MCTS.RolloutOracle{Game}(), nworkers=1)
player = MctsPlayer(mcts, niters=1000, τ=ConstSchedule(0.5))

num_games = 250
bar = Progress(num_games)
avgz = AlphaZero.pit(player, baseline, num_games, gamma=1.0,
    reset_every=nothing,
    color_policy=ALTERNATE_COLORS,
    flip_probability=0.) do i, z, t
  #AlphaZero.debug_trace(t)
  next!(bar)
end
println("Average reward: $avgz")

max_stars_repo_name: btheodorou/Bridge-NN-Research

max_issues_repo_name: btheodorou/Bridge-NN-Research

max_forks_repo_name: btheodorou/Bridge-NN-Research


text:
# sources are applied additively
@generated function atmos_source!(
    stuple::Tuple,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    N = fieldcount(stuple)
    return quote
        Base.Cartesian.@nexprs $N i -> atmos_source!(
            stuple[i],
            atmos,
            source,
            state,
            diffusive,
            aux,
            t,
            direction,
        )
        return nothing
    end
end

function atmos_source!(
    ::Gravity,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    # Migrated to Σsources
end

function atmos_source!(
    ::Coriolis,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    # Migrated to Σsources
end

function atmos_source!(
    subsidence::Subsidence,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    # Migrated to Σsources
end

function atmos_source!(
    s::GeostrophicForcing,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    # Migrated to Σsources
end

function atmos_source!(
    s::RayleighSponge,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    # Migrated to Σsources
end

function atmos_source!(
    ::CreateClouds,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    # Migrated to Σsources
end

function atmos_source!(
    s::RemovePrecipitation,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    # Migrated to Σsources
end

function atmos_source!(
    ::WarmRain_1M,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    # Migrated to Σsources
end

function atmos_source!(
    ::RainSnow_1M,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    # Migrated to Σsources
end

max_forks_repo_name: jm-c/CLIMA


text:
using FileIO, ImageIO, PyPlot

## lena reg intro
img, imgw, flow = gen_init(flipped=true, flow_args=[30])

imgshow(img, figtitle="Target", origin_bot=true)
savefig("../plots/intro_lena_orig.png")

imgshow(imgw, figtitle="Source", origin_bot=true)
savefig("../plots/intro_lena_warped.png")
#
# img_viewed = viewimg(imgw)
# save("../plots/tmp.png", img_viewed)

imgoverlay(img, imgw, figtitle="Blending of Target and Source", origin_bot=true)
savefig("../plots/intro_lena_bl_or_movement.png")
overlay = imgoverlay_v2(img, imgw, flipped=true)
save("../plots/intro_lena_bl_or_movement.png", overlay)

showflow(flow.*(-1), figtitle="Displacement Field")
savefig("../plots/intro_lena_flow.png")

## lena reg intro with red rectangle
using PyCall
using PyPlot: savefig
@pyimport matplotlib.patches as patches

function add_rect(ax)
    rect = patches.Rectangle((90,140),30,30,linewidth=1,edgecolor="r",facecolor="none")
    ax.add_patch(rect)
    return gcf()
end

img, imgw, flow = gen_init(:lena, :tiled, flipped=true)

imgshow(img, origin_bot=true)
imgshow(imgw, origin_bot=true)
showflow(flow)

imgshow(img, figtitle="Image 1", origin_bot=true)
add_rect(gca())
savefig("../plots/rect_intro_lena_orig.png")

imgshow(imgw, figtitle="Image 2", origin_bot=true)
add_rect(gca())
savefig("../plots/rect_intro_lena_warped.png")

imgoverlay(img, imgw, figtitle="Blending of Target and Source", origin_bot=true)
add_rect(gca())
savefig("../plots/rect_intro_lena_bl_or_movement.png")

showflow(flow, figtitle="Displacement Field")
add_rect(gca())
savefig("../plots/rect_intro_lena_flow.png")

# rect disp
flow_rect = flow[140:170, 90:120]
showflow(flow_rect, figtitle="Displacement Field of Rectangle")
rect = patches.Rectangle((1,1),30,30,linewidth=1,edgecolor="r",facecolor="none")
ax = gca()
ax.add_patch(rect)
gcf()
savefig("../plots/rect_intro_lena_flow_zoom.png")

## lap const displacement
img, imgw, flow = gen_init(:lena, :uniform, flow_args=[2+1im, 20], flipped=true)

imgshow(img, figtitle="Target", origin_bot=true)
savefig("../plots/lap_const_lena_orig.png")

imgshow(imgw, figtitle="Source", origin_bot=true)
savefig("../plots/lap_const_lena_warped.png")

imgoverlay(img, imgw, figtitle="Blending of Target and Source")
savefig("../plots/lap_const_lena_bl_or_movement.png")

overlay = imgoverlay_v2(img, imgw, flipped=true)
save("../plots/lap_const_lena_bl_or_movement.png", overlay)

showflow(flow.*(-1), figtitle="Constant Displacement Field")
savefig("../plots/lap_const_lena_flow.png")

smooth = ImageFiltering.Kernel.gaussian(1)

## Make images for testing

img, imgw, flow = gen_init(:lena, flow_args=[20, 60])

showflow(flow)

imgshow(imgw)

using Images, FileIO, ImageIO

save("../testimages/smoothly_varying/target.png", img)
save("../testimages/smoothly_varying/source.png", imgw)

## make tests images:
using Printf
#params:
num_tests = 3
basedir = "../testimages/"
folder = "settings_1/"
flow_args_smoothly_var = [20, 60]
flow_args_const = [20, 300]
img_select = :chess

if img_select == :chess
    folder_name = "chess_"
elseif img_select == :lena
    folder_name = "lena_"
end

for test in 1:num_tests
    name = folder_name * string(test) * "/"
    base_path = basedir*folder*name

    for setting in ["smoothly_varying", "constant"]
        if setting == "smoothly_varying"
            img, imgw, flow = gen_init(img_select, flow_args=flow_args_smoothly_var)
        elseif setting == "constant"
            img, imgw, flow = gen_init(img_select, flow_args=flow_args_const)
        end
        path = base_path * setting * "/"
        mkpath(path);

        # save the max displacement into a file
        f = open(path * "max_disp.txt", "w")
        magnitudes = map(x -> LAP_julia.vec_len(x), flow)
        max_mag = maximum(filter(!isnan, magnitudes))
        max_disp_str = @sprintf("%0.2f", max_mag)
        write(f, max_disp_str)
        close(f)

        save(path * "/target.png", img)
        save(path * "/source.png", imgw)
    end
end

## evaluate tests

methods = ["bunwarpj", "pflap_matlab"]

img, imgw, flow = gen_init(flow_args=[20, 120])
showflow(flow)

u_est, source_reg = pflap(img, imgw)

LAP_julia.mean(abs.(img .- imgw))
LAP_julia.mean(abs.(img .- source_reg))

dist_euclid = sqrt(sum((img .- imgw).^2))/length(img)
dist_euclid = sqrt(sum((img .- source_reg).^2))/length(img)

ssd(img, imgw)
ssd(img, source_reg)
ssd(img, source_reg_p)

imgshow(source_reg)

LAP_julia.mean(source_reg)

u_est_p, source_reg_p = sparse_pflap(img, imgw)

showflow(u_est_p)

imgshow(source_reg_p)
imgshow(img)

LAP_julia.mean(abs.(img .- source_reg_p))
dist_euclid = sqrt(sum(img .- imgw).^2)/length(img)

u_est, source_reg_p_ed, figs, Δ_u = sparse_pflap(img, imgw)

imgshow(source_reg_p_ed)
imgshow(img)

figs[1,5]

## single_lap

u_est = single_lap(img, imgw, 20, [41,41])

showflow(u_est)
LAP_julia.inpaint_nans!(u_est)
u_est_sm = LAP_julia.smooth_with_gaussian!(u_est, 20)

showflow(u_est_sm)

vec, edge = LAP_julia.gradient_magnitude(img)

imgshow(edge, figtitle="Edge Image", origin_bot=true)
savefig("../plots/sparse_lap_edge_image.png")

inds, mag = find_edge_points(img, spacing=10, point_count=300, debug=true)

imgshow(mag .== 0)

imggg = imgoverlay_v2(edge, mag .== 0, flipped=true)

imgshow(edge, figtitle="Edge Image", origin_bot=true)
fig = addpoints(inds)

savefig("../plots/sparse_lap_edge_image_with_points.png")

## sparse lap

img, imgw, flow = gen_init(:lena, flow_args=[10, 150])
showflow(flow, figtitle="Ground Truth Displacement")
savefig("../plots/sparse_lap_ground_truth_flow.png")

fhs = 15
window_size = [31, 31]

mask = parent(padarray(trues(size(img).-(2*fhs, 2*fhs)), Fill(false, (fhs, fhs), (fhs, fhs))))
inds = find_edge_points(img, mask=mask)
points = LAP_julia.inds_to_points(inds)
new_estim = single_lap_at_points(img, imgw, fhs, window_size, points, 3)
showflow(new_estim, disp_type=:sparse, figtitle="Estimated Displacement Vectors")
savefig("../plots/sparse_lap_estimate_image.png")

full_new_estim = interpolate_flow(new_estim, inds)
showflow(full_new_estim, figtitle="Interpolated Displacement Vectors")
savefig("../plots/sparse_lap_estimate_interpolated_image.png")

imgshow(edge, figtitle="Edge Image with Points")
PyPlot.scatter([ind[2] for ind in inds], [ind[1] for ind in inds], marker = :x); gcf()
savefig("../plots/sparse_lap_edge_image_with_points.png")

img, imgw, flow = gen_init(flipped=true)

inds = find_edge_points(img, point_count= 300)

asd = zeros(size(img))

for ind in inds
    fill_square!(asd, Tuple(ind), 3, 1)
end

sum(asd)/length(img)

ver = imgoverlay_v2(asd, img, flipped=true)
save("../plots/convolution_area_small.png", ver)

imgshow(asd, figtitle="Convolution Area", origin_bot = true)
savefig("../plots/convolution_area.png")

function fill_square!(mag, pos, spacing, marker)
    mag_size = size(mag)
    for k = (pos[1]-spacing <= 1 ? 1 : pos[1]-spacing) : (pos[1]+spacing >= mag_size[1] ? mag_size[1] : pos[1]+spacing),
        l = (pos[2]-spacing <= 1 ? 1 : pos[2]-spacing) : (pos[2]+spacing >= mag_size[2] ? mag_size[2] : pos[2]+spacing)
        mag[k, l] = marker
    end
end

N = 20
A = rand(5,5,N)
b = rand(5,N)

@time LAP_julia.multi_mat_div_gem(A, b)
@time LAP_julia.multi_mat_div_qr(A, b)


text:
abstract type Basis{T} <: LazyQuasiMatrix{T} end
abstract type Weight{T} <: LazyQuasiVector{T} end

const WeightedBasis{T, A<:AbstractQuasiVector, B<:Basis} = BroadcastQuasiMatrix{T,typeof(*),<:Tuple{A,B}}

struct WeightLayout <: MemoryLayout end
abstract type AbstractBasisLayout <: MemoryLayout end
struct BasisLayout <: AbstractBasisLayout end
struct SubBasisLayout <: AbstractBasisLayout end
struct MappedBasisLayout <: AbstractBasisLayout end

abstract type AbstractAdjointBasisLayout <: MemoryLayout end
struct AdjointBasisLayout <: AbstractAdjointBasisLayout end
struct AdjointSubBasisLayout <: AbstractAdjointBasisLayout end
struct AdjointMappedBasisLayout <: AbstractAdjointBasisLayout end

MemoryLayout(::Type{<:Basis}) = BasisLayout()
MemoryLayout(::Type{<:Weight}) = WeightLayout()

adjointlayout(::Type, ::BasisLayout) = AdjointBasisLayout()
adjointlayout(::Type, ::SubBasisLayout) = AdjointSubBasisLayout()
adjointlayout(::Type, ::MappedBasisLayout) = AdjointMappedBasisLayout()
broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::BasisLayout) = BasisLayout()
broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::SubBasisLayout) = SubBasisLayout()

combine_mul_styles(::AbstractBasisLayout) = LazyQuasiArrayApplyStyle()
combine_mul_styles(::AbstractAdjointBasisLayout) = LazyQuasiArrayApplyStyle()

ApplyStyle(::typeof(pinv), ::Type{<:Basis}) = LazyQuasiArrayApplyStyle()
pinv(J::Basis) = apply(pinv,J)

_multup(a::Tuple) = Mul(a...)
_multup(a) = a

function ==(A::Basis, B::Basis)
    axes(A) == axes(B) && throw(ArgumentError("Override == to compare bases of type $(typeof(A)) and $(typeof(B))"))
    false
end

ApplyStyle(::typeof(\), ::Type{<:Basis}, ::Type{<:AbstractQuasiMatrix}) = LdivApplyStyle()
ApplyStyle(::typeof(\), ::Type{<:Basis}, ::Type{<:AbstractQuasiVector}) = LdivApplyStyle()
ApplyStyle(::typeof(\), ::Type{<:SubQuasiArray{<:Any,2,<:Basis}}, ::Type{<:AbstractQuasiMatrix}) = LdivApplyStyle()
ApplyStyle(::typeof(\), ::Type{<:SubQuasiArray{<:Any,2,<:Basis}}, ::Type{<:AbstractQuasiVector}) = LdivApplyStyle()

copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(+)}}) = +(broadcast(\,Ref(L.A),arguments(L.B))...)
copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(+)},<:Any,<:AbstractQuasiVector}) =
    +(broadcast(\,Ref(L.A),arguments(L.B))...)

function copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(-)}})
    a,b = arguments(L.B)
    (L.A\a)-(L.A\b)
end

function copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(-)},<:Any,<:AbstractQuasiVector})
    a,b = arguments(L.B)
    (L.A\a)-(L.A\b)
end

function copy(P::Ldiv{BasisLayout,BasisLayout})
    A, B = P.A, P.B
    A == B || throw(ArgumentError("Override materialize for $(typeof(A)) \\ $(typeof(B))"))
    Eye(size(A,2))
end
function copy(P::Ldiv{SubBasisLayout,SubBasisLayout})
    A, B = P.A, P.B
    (parent(A) == parent(B) && parentindices(A) == parentindices(B)) ||
        throw(ArgumentError("Override materialize for $(typeof(A)) \\ $(typeof(B))"))
    Eye(size(A,2))
end

function copy(P::Ldiv{MappedBasisLayout,MappedBasisLayout})
    A, B = P.A, P.B
    demap(A)\demap(B)
end
# function copy(P::Ldiv{MappedBasisLayout,SubBasisLayout})
#     A, B = P.A, P.B
#     # use lazy_getindex to avoid sparse arrays
#     lazy_getindex(parent(A)\parent(B),:,parentindices(B)[2])
# end

for Bas1 in (:Basis, :WeightedBasis), Bas2 in (:Basis, :WeightedBasis)
    @eval ==(A::SubQuasiArray{<:Any,2,<:$Bas1}, B::SubQuasiArray{<:Any,2,<:$Bas2}) =
        all(parentindices(A) == parentindices(B)) && parent(A) == parent(B)
end

# expansion
grid(P) = error("Overload Grid")
function transform(L)
    p = grid(L)
    p,L[p,:]
end

function copy(L::Ldiv{<:AbstractBasisLayout,<:Any,<:Any,<:AbstractQuasiVector})
    p,T = transform(L.A)
    T \ L.B[p]
end

copy(L::Ldiv{<:AbstractBasisLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) =
    copy(Ldiv{LazyLayout,ApplyLayout{typeof(*)}}(L.A, L.B))

function copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(*)},<:AbstractQuasiMatrix,<:AbstractQuasiVector})
    p,T = transform(L.A)
     T \ L.B[p]
end

## materialize views

# materialize(S::SubQuasiArray{<:Any,2,<:ApplyQuasiArray{<:Any,2,typeof(*),<:Tuple{<:Basis,<:Any}}}) =
#     *(arguments(S)...)

# Differentiation of sub-arrays
function copy(M::QMul2{<:Derivative,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:Inclusion,<:Any}}})
    A, B = M.args
    P = parent(B)
    (Derivative(axes(P,1))*P)[parentindices(B)...]
end

function copy(M::QMul2{<:Derivative,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AffineQuasiVector,<:Any}}})
    A, B = M.args
    P = parent(B)
    kr,jr = parentindices(B)
    (Derivative(axes(P,1))*P*kr.A)[kr,jr]
end

function copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(*)},<:AbstractQuasiMatrix})
    args = arguments(L.B)
    # this is a temporary hack
    if args isa Tuple{AbstractQuasiMatrix,Number}
        (L.A \  first(args))*last(args)
    elseif args isa Tuple{Number,AbstractQuasiMatrix}
        first(args)*(L.A \ last(args))
    else
        error("Not implemented")
    end
end

# we represent as a Mul with a banded matrix
sublayout(::AbstractBasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractUnitRange}}) = SubBasisLayout()
sublayout(::BasisLayout, ::Type{<:Tuple{<:AffineQuasiVector,<:AbstractUnitRange}}) = MappedBasisLayout()

demap(x) = x
demap(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Any,<:Slice}}) = parent(V)
function demap(V::SubQuasiArray{<:Any,2})
    kr, jr = parentindices(V)
    demap(parent(V)[kr,:])[:,jr]
end

##
# SubLayout behaves like ApplyLayout{typeof(*)}

combine_mul_styles(::SubBasisLayout) = combine_mul_styles(ApplyLayout{typeof(*)}())
_arguments(::SubBasisLayout, A) = _arguments(ApplyLayout{typeof(*)}(), A)
call(::SubBasisLayout, ::SubQuasiArray) = *

combine_mul_styles(::AdjointSubBasisLayout) = combine_mul_styles(ApplyLayout{typeof(*)}())
_arguments(::AdjointSubBasisLayout, A) = _arguments(ApplyLayout{typeof(*)}(), A)
arguments(::AdjointSubBasisLayout, A) = arguments(ApplyLayout{typeof(*)}(), A)
call(::AdjointSubBasisLayout, ::SubQuasiArray) = *

function arguments(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Inclusion,<:AbstractUnitRange}})
    A = parent(V)
    _,jr = parentindices(V)
    first(jr) ≥ 1 || throw(BoundsError())
    P = _BandedMatrix(Ones{Int}(1,length(jr)), axes(A,2), first(jr)-1,1-first(jr))
    A,P
end

copy(L::Ldiv{BasisLayout,SubBasisLayout}) = apply(\, L.A, ApplyQuasiArray(L.B))

include("splines.jl")