text:
module FormattingFuncs

export format_integer, format_float, format_number, format_numbers

using DataFrames
using Printf

"""
Format an integer number as a string with a `thousands separator` that can be specified
(default `.`).
"""
function format_integer(number::Integer;
                        thousandsSeparator::String=".")::String
    # Format number as string for easier processing
    negative_number::Bool = number < 0
    number_string::String = @sprintf("%.0f", negative_number ? -1 * number : number)
    # start = start position of processing, stop = end position of processing
    start::Int64 = 1 + length(number_string) % 3
    stop::Int64 = length(number_string) - length(number_string) % 3
    # Get the various parts of the number that need to be separated
    results::Vector{String} = Vector{String}()
    if number_string[1:start-1] != ""
        push!(results, number_string[1:start-1])
    end
    results = [results ; [number_string[i:i+2] for i in start:3:stop]]
    # Return the formatted integer as a string
    if negative_number
        return "-" * join(results, thousandsSeparator)
    end
    return join(results, thousandsSeparator)
end

"""
Format a floating point number with a `decimal separator` and a `thousands separator`,
both of which can be specified. Defaults: decimal `,`, thousands `.`.
"""
function format_float(number::Float64;
                      thousandsSeparator::String=".",
                      decimalSeparator::String=",",
                      decimalPlaces::Int64=2)::String
    # Get the fractional and integral parts of the specified number as strings
    fpart, ipart = modf(round(number, digits=decimalPlaces))
    negative_number::Bool = number < 0
    ipart_string::String = format_integer(Int(ipart), thousandsSeparator=thousandsSeparator)
    fpart_string::String = @sprintf("%.50f", negative_number ? -1 * fpart : fpart)
    # Determine length of fractional part
    fpart_length = minimum([2 + decimalPlaces, length(fpart_string)])
    # Return formatted float number as a string
    if decimalPlaces == 0
        return ipart_string
    end
    return ipart_string * decimalSeparator * fpart_string[3:fpart_length]
end

"""
Format an integer or a floating point number as a string with a `decimal separator` and a
`thousands separator`, both of which can be specified. Defaults: decimal `,`, thousands `.`.

For floating point numbers the number of `decimal places` can be specified (default `2`).

For integer numbers there is an option to `show decimals` (default `false`).
"""
function format_number(number::T;
                       thousandsSeparator::String=".",
                       decimalSeparator::String=",",
                       decimalPlaces::Int64=2,
                       intShowDecimals::Bool=false)::String where T<:Real
    # Determine number type and call appropriate formatting function
    if isa(number, Integer)
        if intShowDecimals
            return format_float(Float64(number), thousandsSeparator=thousandsSeparator,
                                decimalSeparator=decimalSeparator, decimalPlaces=decimalPlaces)
        else
            return format_integer(Int(number), thousandsSeparator=thousandsSeparator)
        end
    elseif isa(number, Float16) || isa(number, Float32) || isa(number, Float64)
        return format_float(Float64(number), thousandsSeparator=thousandsSeparator,
                            decimalSeparator=decimalSeparator, decimalPlaces=decimalPlaces)
    else
        return string(number)
    end
end

"""
Format a vector of numbers (integer and / or floating point) as strings with a `decimal separator`
and a `thousands separator`, both of which can be specified. Defaults: decimal `,`, thousands `.`.

When a `DataFrame` is passed, all number columns will be formatted as strings: the same separators
as mentioned for vectors of numbers apply.

For floating point numbers the number of `decimal places` can be specified (default `2`).

For integer numbers there is an option to `show decimals` (default `false`).
"""
function format_numbers(v::Vector{T};
                        thousandsSeparator::String=".",
                        decimalSeparator::String=",",
                        decimalPlaces::Int64=2,
                        intShowDecimals::Bool=false)::Vector{String} where T<:Real
    return format_number.(v, thousandsSeparator=thousandsSeparator, decimalSeparator=decimalSeparator,
                          decimalPlaces=decimalPlaces, intShowDecimals=intShowDecimals)
end

function format_numbers(df::DataFrame;
                        thousandsSeparator::String=".",
                        decimalSeparator::String=",",
                        decimalPlaces::Int64=2,
                        intShowDecimals::Bool=false)::DataFrame
    dfc = copy(df)
    # Format number colums (type -> String)
    number_cols = propertynames(select(dfc, findall(col -> all(value -> value isa Union{Missing, Number}, col), eachcol(dfc))))
    foreach(col -> dfc[!, col] = format_numbers(dfc[!, col], thousandsSeparator=thousandsSeparator, decimalSeparator=decimalSeparator,
                                                decimalPlaces=decimalPlaces, intShowDecimals=intShowDecimals), number_cols)
    return dfc
end

end

max_stars_repo_name: bogersw/julia-snippets

max_issues_repo_name: bogersw/julia-snippets

max_forks_repo_name: bogersw/julia-snippets


text:
module Architectures

export AbstractArchitecture, AbstractMultiArchitecture
export CPU, GPU, CUDAGPU, ROCMGPU
export device, device_event, architecture, array_type, arch_array

using CUDA
using AMDGPU
using KernelAbstractions
using CUDAKernels
using ROCKernels
using Adapt
using OffsetArrays

"""
    AbstractArchitecture

Abstract supertype for architectures supported by Oceananigans.
"""
abstract type AbstractArchitecture end

"""
    AbstractMultiArchitecture

Abstract supertype for Distributed architectures supported by Oceananigans.
"""
abstract type AbstractMultiArchitecture <: AbstractArchitecture end

"""
    CPU <: AbstractArchitecture

Run Oceananigans on one CPU node. Uses multiple threads if the environment
variable `JULIA_NUM_THREADS` is set.
"""
struct CPU <: AbstractArchitecture end

"""
    GPU <: AbstractArchitecture

Run Oceananigans on a single NVIDIA CUDA GPU.
"""
struct GPU{D} <: AbstractArchitecture
    device :: D
end

const CUDAGPU = GPU{<:CUDAKernels.CUDADevice}
const ROCMGPU = GPU{<:ROCKernels.ROCDevice}

# Convenience, non-public constructors (may be better to remove these eventually for code clarity)
CUDAGPU() = GPU(CUDAKernels.CUDADevice())
ROCMGPU() = GPU(ROCKernels.ROCDevice())

#####
##### These methods are extended in Distributed.jl
#####

device(::CPU) = KernelAbstractions.CPU()
device(::CUDAGPU) = CUDAKernels.CUDADevice()
device(::ROCMGPU) = ROCKernels.ROCDevice()

architecture() = nothing
architecture(::Number) = nothing
architecture(::Array) = CPU()
architecture(::CuArray) = CUDAGPU()
architecture(::ROCArray) = ROCMGPU()
architecture(a::SubArray) = architecture(parent(a))
architecture(a::OffsetArray) = architecture(parent(a))

"""
    child_architecture(arch)

Return `arch`itecture of child processes.
On single-process, non-distributed systems, return `arch`.
"""
child_architecture(arch) = arch

array_type(::CPU) = Array
array_type(::CUDAGPU) = CuArray
array_type(::ROCMGPU) = ROCArray

arch_array(::CPU, a::Array)   = a
arch_array(::CPU, a::CuArray) = Array(a)
arch_array(::CPU, a::ROCArray) = Array(a)
arch_array(::CUDAGPU, a::Array)   = CuArray(a)
arch_array(::CUDAGPU, a::CuArray) = a
arch_array(::ROCMGPU, a::Array) = ROCArray(a)
arch_array(::ROCMGPU, a::ROCArray) = a

arch_array(arch, a::AbstractRange) = a
arch_array(arch, a::OffsetArray) = OffsetArray(arch_array(arch, a.parent), a.offsets...)
arch_array(arch, ::Nothing) = nothing
arch_array(arch, a::Number) = a

device_event(arch) = Event(device(arch))

end # module

ext: jl


text:
"""
Tools for Pauli strings and binary symplectic vectors / matrices.
"""
module PauliTools

# exports
export bsp, pack, unpack, to_bsf, to_pauli, weight

@doc raw"""
    bsp(
        A::AbstractVecOrMat{Bool}, B::AbstractVecOrMat{Bool}
    ) -> Union{Bool,BitVector,BitMatrix}

Return the binary symplectic product of `A` with `B`, given in binary symplectic form.

The binary symplectic product ``⊙`` is defined as ``A ⊙ B ≡ A Λ B \bmod 2`` where
``Λ = \left[\begin{smallmatrix} 0 & I \\ I & 0 \end{smallmatrix}\right]``.

# Examples
```jldoctest
julia> a = BitVector([1, 0, 0, 0]);  # XI

julia> b = BitVector([0, 0, 1, 0]);  # ZI

julia> bsp(a', b)
true
```
```jldoctest
julia> stabilizers = BitMatrix(  # 5-qubit stabilizers
       [1 0 0 1 0 0 1 1 0 0     # XZZXI
        0 1 0 0 1 0 0 1 1 0     # IXZZX
        1 0 1 0 0 0 0 0 1 1     # XIXZZ
        0 1 0 1 0 1 0 0 0 1]);  # ZXIXZ

julia> error = BitVector([0, 0, 1, 1, 0, 0, 1, 0, 1, 0]);  # IZXYI

julia> bsp(stabilizers, error)
4-element BitVector:
 0
 1
 1
 0
```
"""
function bsp(a::AbstractVecOrMat{Bool}, b::AbstractVecOrMat{Bool})
    # circshift b by half its 1st dimension to emulate symplectic product
    return isodd.(a * circshift(b, size(b, 1)/2))  # mod elements to base 2
end

"""
    pack(bsf::AbstractVector{Bool}) -> Tuple{String,Int}

Pack a binary vector into a concise representation, typically for log output. See also
[`unpack`](@ref).

# Examples
```jldoctest
julia> a = BitVector([1, 0, 1, 0, 1, 1]);  # XZY

julia> b = pack(a)  # (hex_value, length)
("2b", 6)
julia> unpack(b) == a
true
```
"""
function pack(bsf::AbstractVector{Bool})
    bit_str = join(d ? '1' : '0' for d in bsf)  # convert to bit-string
    val = parse(BigInt, bit_str; base=2)  # parse to val (BigInt to support long vectors)
    return string(val, base=16), length(bsf)  # hex_str, len
end

"""
    unpack(packed_bsf::Tuple{String,Int}) -> BitVector

Unpack a binary vector from a concise representation, typically from log output. See also
[`pack`](@ref).

# Examples
```jldoctest
julia> a = ("2b", 6);  # (hex_value, length)

julia> b = unpack(a)  # XZY
6-element BitVector:
 1
 0
 1
 0
 1
 1
julia> pack(b) == a
true
```
"""
function unpack(packed_bsf::Tuple{String,Int})
    hex_str, len = packed_bsf
    val = parse(BigInt, hex_str; base=16)  # parse to val (BigInt to support long vectors)
    return reverse!(BitVector(digits(val; base=2, pad=len)))  # binary-digits as BitVector
end

"""
    to_bsf(
        pauli::Union{AbstractString,AbstractVector{<:AbstractString}}
    ) -> Union{BitVector,BitMatrix}

Convert the Pauli string operator(s) to binary symplectic form.

A single Pauli string is converted to a vector. A vector of Pauli strings is converted to a
matrix where each row corresponds to a Pauli.

# Examples
```jldoctest
julia> to_bsf("XIZIY")
10-element BitVector:
 1
 0
 0
 0
 1
 0
 0
 1
 0
 1
```
```jldoctest
julia> to_bsf(["XIZIY", "IXZYI"])
2×10 BitMatrix:
 1  0  0  0  1  0  0  1  0  1
 0  1  0  1  0  0  0  1  1  0
```
"""
function to_bsf(pauli::AbstractString)
    p = collect(pauli)
    ys = p .== 'Y'
    return vcat((p .== 'X') .| ys, (p .== 'Z') .| ys)
end
function to_bsf(paulis::AbstractVector{<:AbstractString})
    # for each pauli string, convert to bsf vector, transpose; then concatenate as rows
    return vcat((transpose(to_bsf(p)) for p in paulis)...)
end

"""
    to_pauli(bsf::AbstractVecOrMat{Bool}) -> Union{String,Vector{String}}

Convert the binary symplectic form to Pauli string operator(s).

A vector is converted to a single Pauli string. A matrix is converted row-by-row to a
collection of Pauli strings.

# Examples
```jldoctest
julia> to_pauli(BitVector([1, 0, 0, 0, 1, 0, 0, 1, 0, 1]))
"XIZIY"
```
```jldoctest
julia> to_pauli(BitMatrix([1 0 0 0 1 0 0 1 0 1; 0 1 0 1 0 0 0 1 1 0]))
2-element Vector{String}:
 "XIZIY"
 "IXZYI"
```
"""
function to_pauli(bsf::AbstractVector{Bool})
    l = length(bsf)  # bsf = (1 0 0 0 1 | 0 0 1 0 1)
    x, z = bsf[1:l÷2], bsf[l÷2+1:end]  # x = (1 0 0 0 1), z = (0 0 1 0 1)
    return String(map(i -> "IXZY"[i+1], x+2z))  # x+2z = (1 0 2 0 3) -> "XIZIY"
end
function to_pauli(bsfs::AbstractMatrix{Bool})
    return String[to_pauli(bsf) for bsf in eachrow(bsfs)]
end

"""
    weight(bsf::AbstractVecOrMat{Bool}) -> Int

Return the weight of the binary symplectic form.

# Examples
```jldoctest
julia> weight(BitVector([1, 0, 0, 0, 1, 0, 0, 1, 0, 1]))
3
```
```jldoctest
julia> weight(BitMatrix([1 0 0 0 1 0 0 1 0 1; 1 1 1 1 1 0 0 0 0 0]))
8
```
"""
function weight(bsf::AbstractVector{Bool})
    l = length(bsf)  # bsf = (1 0 0 | 1 1 0)
    return count(>(0), bsf[1:l÷2] + bsf[l÷2+1:end])  # count(>(0), (2 1 0))
end
function weight(bsfs::AbstractMatrix{Bool})
    l = size(bsfs, 2)  # bsf = (1 0 0 | 1 1 0; 1 1 1 | 0 1 0)
    return count(>(0), bsfs[:,1:l÷2] + bsfs[:,l÷2+1:end])  # count(>(0), (2 1 0; 1 2 1))
end

"""
    weight(pauli::Union{AbstractString,AbstractVector{<:AbstractString}}) -> Int

Return the weight of the Pauli string operator(s).

# Examples
```jldoctest
julia> weight("XIZIY")
3
```
```jldoctest
julia> weight(["XIZIY", "XXXXX"])
8
```
"""
function weight(pauli::AbstractString)
    return count(!=('I'), pauli)  # count(!=('I'), "IXYIZ")
end
function weight(paulis::AbstractVector{<:AbstractString})
    return sum(weight, paulis)
end

end


text:
export AbstractSymmetry
# export AbstractUnitSymmetry

"""
    AbstractSymmetry{OperationType<:AbstractSpaceSymmetryOperation}

Abstract symmetry whose elements are of type `OperationType`.
"""
abstract type AbstractSymmetry{OperationType<:AbstractSpaceSymmetryOperation} end
# abstract type AbstractUnitSymmetry{OperationType<:AbstractSpaceSymmetryOperation} <: AbstractSymmetry{OperationType} end


text:
std = 1.0

@gen (static) function kernel(t::Int, x_prev::Float64, (grad)(alpha::Float64), (grad)(beta::Float64))
    x = @trace(normal(x_prev * alpha + beta, std), :x)
    return x
end

foo = Unfold(kernel)

@load_generated_functions()

@testset "unfold combinator" begin

    @testset "Julia call" begin
        @test length(foo(5, 0., 1.0, 1.0)) == 5
    end

    @testset "simulate" begin
        x_init = 0.1
        alpha = 0.2
        beta = 0.3
        x1 = 1.1
        x2 = 1.2
        x3 = 1.3
        trace = simulate(foo, (3, x_init, alpha, beta))
        x1 = trace[1 => :x]
        x2 = trace[2 => :x]
        x3 = trace[3 => :x]
        choices = get_choices(trace)
        @test length(collect(get_values_shallow(choices))) == 0
        @test length(collect(get_submaps_shallow(choices))) == 3
        expected_score = (logpdf(normal, x1, x_init * alpha + beta, std)
                          + logpdf(normal, x2, x1 * alpha + beta, std)
                          + logpdf(normal, x3, x2 * alpha + beta, std))
        @test isapprox(get_score(trace), expected_score)
        retval = get_retval(trace)
        @test length(retval) == 3
        @test retval[1] == x1
        @test retval[2] == x2
        @test retval[3] == x3
    end

    @testset "generate" begin
        x_init = 0.1
        alpha = 0.2
        beta = 0.3
        x1 = 1.1
        x2 = 1.2
        x3 = 1.3
        constraints = choicemap()
        constraints[1 => :x] = x1
        constraints[3 => :x] = x3
        (trace, weight) = generate(foo, (3, x_init, alpha, beta), constraints)
        choices = get_choices(trace)
        @test choices[1 => :x] == x1
        @test choices[3 => :x] == x3
        @test length(collect(get_values_shallow(choices))) == 0
        @test length(collect(get_submaps_shallow(choices))) == 3
        x2 = choices[2 => :x]
        expected_weight = (logpdf(normal, x1, x_init * alpha + beta, std)
                           + logpdf(normal, x3, x2 * alpha + beta, std))
        expected_score = (logpdf(normal, x1, x_init * alpha + beta, std)
                          + logpdf(normal, x2, x1 * alpha + beta, std)
                          + logpdf(normal, x3, x2 * alpha + beta, std))
        @test isapprox(weight, expected_weight)
        @test isapprox(get_score(trace), expected_score)
        retval = get_retval(trace)
        @test length(retval) == 3
        @test retval[1] == x1
        @test retval[2] == x2
        @test retval[3] == x3
    end

    @testset "propose" begin
        x_init = 0.1
        alpha = 0.2
        beta = 0.3
        (choices, weight, retval) = propose(foo, (3, x_init, alpha, beta))
        @test length(collect(get_values_shallow(choices))) == 0
        @test length(collect(get_submaps_shallow(choices))) == 3
        x1 = choices[1 => :x]
        x2 = choices[2 => :x]
        x3 = choices[3 => :x]
        expected_weight = (logpdf(normal, x1, x_init * alpha + beta, std)
                           + logpdf(normal, x2, x1 * alpha + beta, std)
                           + logpdf(normal, x3, x2 * alpha + beta, std))
        @test isapprox(weight, expected_weight)
        @test length(retval) == 3
        @test retval[1] == x1
        @test retval[2] == x2
        @test retval[3] == x3
    end

    @testset "assess" begin
        x_init = 0.1
        alpha = 0.2
        beta = 0.3
        x1 = 1.1
        x2 = 1.2
        x3 = 1.3
        choices = choicemap()
        choices[1 => :x] = x1
        choices[2 => :x] = x2
        choices[3 => :x] = x3
        (weight, retval) = assess(foo, (3, x_init, alpha, beta), choices)
        expected_weight = (logpdf(normal, x1, x_init * alpha + beta, std)
                           + logpdf(normal, x2, x1 * alpha + beta, std)
                           + logpdf(normal, x3, x2 * alpha + beta, std))
        @test isapprox(weight, expected_weight)
        @test length(retval) == 3
        @test retval[1] == x1
        @test retval[2] == x2
        @test retval[3] == x3
    end

    @testset "update" begin
        x_init = 0.1
        alpha = 0.2
        beta = 0.3
        x1 = 1.1
        x2 = 1.2

        function get_initial_trace()
            constraints = choicemap()
            constraints[1 => :x] = x1
            constraints[2 => :x] = x2
            (trace, _) = generate(foo, (2, x_init, alpha, beta), constraints)
            trace
        end

        # unknown change to args, increasing length from 2 to 3 and change 2 and change params
        trace = get_initial_trace()
        x2_new = 1.3
        x3_new = 1.4
        alpha_new = 0.5
        constraints = choicemap()
        constraints[2 => :x] = x2_new
        constraints[3 => :x] = x3_new
        (trace, weight, retdiff, discard) = update(trace,
                                                   (3, x_init, alpha_new, beta),
                                                   (UnknownChange(), UnknownChange(), UnknownChange(), UnknownChange()), constraints)
        choices = get_choices(trace)
        @test get_args(trace) == (3, x_init, alpha_new, beta)
        @test choices[1 => :x] == x1
        @test choices[2 => :x] == x2_new
        @test choices[3 => :x] == x3_new
        @test discard[2 => :x] == x2
        expected_score = (logpdf(normal, x1, x_init * alpha_new + beta, std)
                          + logpdf(normal, x2_new, x1 * alpha_new + beta, std)
                          + logpdf(normal, x3_new, x2_new * alpha_new + beta, std))
        @test isapprox(get_score(trace), expected_score)
        expected_weight = (logpdf(normal, x3_new, x2_new * alpha_new + beta, std)
                           + logpdf(normal, x2_new, x1 * alpha_new + beta, std)
                           - logpdf(normal, x2, x1 * alpha + beta, std)
                           + logpdf(normal, x1, x_init * alpha_new + beta, std)
                           - logpdf(normal, x1, x_init * alpha + beta, std))
        @test isapprox(weight, expected_weight)
        retval = get_retval(trace)
        @test length(retval) == 3
        @test retval[1] == x1
        @test retval[2] == x2_new
        @test retval[3] == x3_new
        @test isa(retdiff, VectorDiff)
        @test !haskey(retdiff.updated, 1) # no diff
        @test retdiff.updated[2] == UnknownChange() # x changed
        @test !haskey(retdiff.updated, 3) # new, not retained

        # unknown change to args, decreasing length from 2 to 1 and change 1 and change params
        trace = get_initial_trace()
        x1_new = 1.3
        alpha_new = 0.5
        constraints = choicemap()
        constraints[1 => :x] = x1_new
        (trace, weight, retdiff, discard) = update(trace,
                                                   (1, x_init, alpha_new, beta),
                                                   (UnknownChange(), UnknownChange(), UnknownChange(), UnknownChange()), constraints)
        choices = get_choices(trace)
        @test get_args(trace) == (1, x_init, alpha_new, beta)
        @test !has_value(choices, 2 => :x)
        @test !has_value(choices, 3 => :x)
        @test choices[1 => :x] == x1_new
        @test discard[1 => :x] == x1
        @test discard[2 => :x] == x2
        @test isapprox(get_score(trace), logpdf(normal, x1_new, x_init * alpha_new + beta, std))
        expected_weight = (logpdf(normal, x1_new, x_init * alpha_new + beta, std)
                           - logpdf(normal, x1, x_init * alpha + beta, std)
                           - logpdf(normal, x2, x1 * alpha + beta, std))
        @test isapprox(weight, expected_weight)
        retval = get_retval(trace)
        @test length(retval) == 1
        @test retval[1] == x1_new
        @test isa(retdiff, VectorDiff)
        @test retdiff.updated[1] == UnknownChange()
        @test !haskey(retdiff.updated, 2) # removed

        # increasing length from 2 to 4; constrain 4 and let 3 be generated from prior
        # also change alpha
        trace = get_initial_trace()
        x4_new = 1.3
        alpha_new = 0.5
        constraints = choicemap()
        constraints[4 => :x] = x4_new
        (trace, weight, retdiff, discard) = update(trace,
                                                   (4, x_init, alpha_new, beta),
                                                   (UnknownChange(), UnknownChange(), UnknownChange(), UnknownChange()), constraints)
        choices = get_choices(trace)
        @test isempty(discard)
        @test get_args(trace) == (4, x_init, alpha_new, beta)
        @test choices[1 => :x] == x1
        @test choices[2 => :x] == x2
        x3_new = choices[3 => :x]
        @test choices[4 => :x] == x4_new
        expected_score = (logpdf(normal, x1, x_init * alpha_new + beta, std)
                          + logpdf(normal, x2, x1 * alpha_new + beta, std)
                          + logpdf(normal, x3_new, x2 * alpha_new + beta, std)
                          + logpdf(normal, x4_new, x3_new * alpha_new + beta, std))
        @test isapprox(get_score(trace), expected_score)
        expected_weight = (logpdf(normal, x1, x_init * alpha_new + beta, std)
                           - logpdf(normal, x1, x_init * alpha + beta, std)
                           + logpdf(normal, x2, x1 * alpha_new + beta, std)
                           - logpdf(normal, x2, x1 * alpha + beta, std)
                           + logpdf(normal, x4_new, x3_new * alpha_new + beta, std))
        @test isapprox(weight, expected_weight)
        retval = get_retval(trace)
        @test length(retval) == 4
        @test retval[1] == x1
        @test retval[2] == x2
        @test retval[3] == x3_new
        @test retval[4] == x4_new
        @test isa(retdiff, VectorDiff)
        @test !haskey(retdiff.updated, 1) # no diff
        @test !haskey(retdiff.updated, 2) # no diff
        @test !haskey(retdiff.updated, 3) # new
        @test !haskey(retdiff.updated, 4) # new

        # no change to arguments, change nothing
        trace = get_initial_trace()
        (trace, weight, retdiff, discard) = update(trace,
                                                   (2, x_init, alpha, beta),
                                                   (NoChange(), NoChange(), NoChange(), NoChange()), EmptyChoiceMap())
        choices = get_choices(trace)
        @test get_args(trace) == (2, x_init, alpha, beta)
        @test choices[1 => :x] == x1
        @test choices[2 => :x] == x2
        @test isempty(discard)
        expected_score = (logpdf(normal, x1, x_init * alpha + beta, std)
                          + logpdf(normal, x2, x1 * alpha + beta, std))
        @test isapprox(get_score(trace), expected_score)
        @test isapprox(weight, 0.)
        retval = get_retval(trace)
        @test retdiff == NoChange()

        # no change to arguments, change x2
        trace = get_initial_trace()
        x2_new = 3.3
        constraints = choicemap()
        constraints[2 => :x] = x2_new
        (trace, weight, retdiff, discard) = update(trace,
                                                   (2, x_init, alpha, beta),
                                                   (NoChange(), NoChange(), NoChange(), NoChange()), constraints)
        choices = get_choices(trace)
        @test get_args(trace) == (2, x_init, alpha, beta)
        @test choices[1 => :x] == x1
        @test choices[2 => :x] == x2_new
        @test discard[2 => :x] == x2
        expected_score = (logpdf(normal, x1, x_init * alpha + beta, std)
                          + logpdf(normal, x2_new, x1 * alpha + beta, std))
        @test isapprox(get_score(trace), expected_score)
        expected_weight = (logpdf(normal, x2_new, x1 * alpha + beta, std)
                           - logpdf(normal, x2, x1 * alpha + beta, std))
        @test isapprox(weight, expected_weight)
        retval = get_retval(trace)
        @test length(retval) == 2
        @test retval[1] == x1
        @test retval[2] == x2_new
        @test isa(retdiff, VectorDiff)
        @test !haskey(retdiff.updated, 1) # no change
        @test retdiff.updated[2] == UnknownChange()

        # init_changed=true, params_changed=false, change nothing
        trace = get_initial_trace()
        x_init_new = 0.1
        (trace, weight, retdiff, discard) = update(trace,
                                                   (2, x_init_new, alpha, beta),
                                                   (NoChange(), UnknownChange(), NoChange(), NoChange()), EmptyChoiceMap())
        choices = get_choices(trace)
        @test get_args(trace) == (2, x_init_new, alpha, beta)
        @test choices[1 => :x] == x1
        @test choices[2 => :x] == x2
        @test isempty(discard)
        expected_score = (logpdf(normal, x1, x_init_new * alpha + beta, std)
                          + logpdf(normal, x2, x1 * alpha + beta, std))
        @test isapprox(get_score(trace), expected_score)
        expected_weight = (logpdf(normal, x1, x_init_new * alpha + beta, std)
                           - logpdf(normal, x1, x_init * alpha + beta, std))
        @test isapprox(weight, expected_weight)
        retval = get_retval(trace)
        @test length(retval) == 2
        @test retval[1] == x1
        @test retval[2] == x2
        @test retdiff == NoChange()

        # init_changed=false, params_changed=true, change nothing
        trace = get_initial_trace()
        alpha_new = 0.5
        (trace, weight, retdiff, discard) = update(trace,
                                                   (2, x_init, alpha_new, beta),
                                                   (NoChange(), NoChange(), UnknownChange(), UnknownChange()), EmptyChoiceMap())
        choices = get_choices(trace)
        @test get_args(trace) == (2, x_init, alpha_new, beta)
        @test choices[1 => :x] == x1
        @test choices[2 => :x] == x2
        @test isempty(discard)
        expected_score = (logpdf(normal, x1, x_init * alpha_new + beta, std)
                          + logpdf(normal, x2, x1 * alpha_new + beta, std))
        @test isapprox(get_score(trace), expected_score)
        expected_weight = (logpdf(normal, x1, x_init * alpha_new + beta, std)
                           - logpdf(normal, x1, x_init * alpha + beta, std)
                           + logpdf(normal, x2, x1 * alpha_new + beta, std)
                           - logpdf(normal, x2, x1 * alpha + beta, std))
        @test isapprox(weight, expected_weight)
        retval = get_retval(trace)
        @test length(retval) == 2
        @test retval[1] == x1
        @test retval[2] == x2
        @test retdiff == NoChange()
    end

    @testset "regenerate" begin
        x_init = 0.1
        alpha = 0.2
        beta = 0.3
        x1 = 1.1
        x2 = 1.2

        function get_initial_trace()
            constraints = choicemap()
            constraints[1 => :x] = x1
            constraints[2 => :x] = x2
            (trace, _) = generate(foo, (2, x_init, alpha, beta), constraints)
            trace
        end

        # unknown change to args, increasing length from 2 to 3 and change 2 and change params
        trace = get_initial_trace()
        alpha_new = 0.5
        selection = select(2 => :x)
        (trace, weight, retdiff) = regenerate(trace,
                                              (3, x_init, alpha_new, beta),
                                              (UnknownChange(), UnknownChange(), UnknownChange(), UnknownChange()), selection)
        choices = get_choices(trace)
        @test get_args(trace) == (3, x_init, alpha_new, beta)
        @test choices[1 => :x] == x1
        x2_new = choices[2 => :x]
        x3_new = choices[3 => :x]
        expected_score = (logpdf(normal, x1, x_init * alpha_new + beta, std)
                          + logpdf(normal, x2_new, x1 * alpha_new + beta, std)
                          + logpdf(normal, x3_new, x2_new * alpha_new + beta, std))
        @test isapprox(get_score(trace), expected_score)
        expected_weight = (logpdf(normal, x1, x_init * alpha_new + beta, std)
                           - logpdf(normal, x1, x_init * alpha + beta, std))
        @test isapprox(weight, expected_weight)
        retval = get_retval(trace)
        @test length(retval) == 3
        @test retval[1] == x1
        @test retval[2] == x2_new
        @test retval[3] == x3_new
        @test isa(retdiff, VectorDiff)
        @test !haskey(retdiff.updated, 1) # no diff
        @test retdiff.updated[2] == UnknownChange() # retval changed
        @test !haskey(retdiff.updated, 3) # new, not retained

        # unknown change to args, decreasing length from 2 to 1 and change 1 and change params
        trace = get_initial_trace()
        alpha_new = 0.5
        selection = select(1 => :x)
        (trace, weight, retdiff) = regenerate(trace,
                                              (1, x_init, alpha_new, beta),
                                              (UnknownChange(), UnknownChange(), UnknownChange(), UnknownChange()), selection)
        choices = get_choices(trace)
        @test get_args(trace) == (1, x_init, alpha_new, beta)
        @test !has_value(choices, 2 => :x)
        @test !has_value(choices, 3 => :x)
        x1_new = choices[1 => :x]
        @test isapprox(get_score(trace), logpdf(normal, x1_new, x_init * alpha_new + beta, std))
        @test isapprox(weight, 0.)
        retval = get_retval(trace)
        @test length(retval) == 1
        @test retval[1] == x1_new
        @test isa(retdiff, VectorDiff)
        @test retdiff.updated[1] == UnknownChange() # retval changed
        @test !haskey(retdiff.updated, 2) # removed

        # no change to args, change nothing
        trace = get_initial_trace()
        (trace, weight, retdiff) = regenerate(trace,
                                              (2, x_init, alpha, beta),
                                              (NoChange(), NoChange(), NoChange(), NoChange()), EmptySelection())
        choices = get_choices(trace)
        @test get_args(trace) == (2, x_init, alpha, beta)
        @test choices[1 => :x] == x1
        @test choices[2 => :x] == x2
        expected_score = (logpdf(normal, x1, x_init * alpha + beta, std)
                          + logpdf(normal, x2, x1 * alpha + beta, std))
        @test isapprox(get_score(trace), expected_score)
        @test isapprox(weight, 0.)
        retval = get_retval(trace)
        @test length(retval) == 2
        @test retval[1] == x1
        @test retval[2] == x2
        @test retdiff == NoChange()

        # no change to args, change x2
        trace = get_initial_trace()
        x2_new = 3.3
        selection = select(2 => :x)
        (trace, weight, retdiff) = regenerate(trace,
                                              (2, x_init, alpha, beta),
                                              (NoChange(), NoChange(), NoChange(), NoChange()), selection)
        choices = get_choices(trace)
        @test get_args(trace) == (2, x_init, alpha, beta)
        @test choices[1 => :x] == x1
        x2_new = choices[2 => :x]
        expected_score = (logpdf(normal, x1, x_init * alpha + beta, std)
                          + logpdf(normal, x2_new, x1 * alpha + beta, std))
        @test isapprox(get_score(trace), expected_score)
        @test isapprox(weight, 0.)
        retval = get_retval(trace)
        @test length(retval) == 2
        @test retval[1] == x1
        @test retval[2] == x2_new
        @test isa(retdiff, VectorDiff)
        @test !haskey(retdiff.updated, 1) # no diff
        @test retdiff.updated[2] == UnknownChange() # retval changed

        # init_changed=true, params_changed=false, change nothing
        trace = get_initial_trace()
        x_init_new = -0.1
        (trace, weight, retdiff) = regenerate(trace,
                                              (2, x_init_new, alpha, beta),
                                              (NoChange(), UnknownChange(), NoChange(), NoChange()), EmptySelection())
        choices = get_choices(trace)
        @test get_args(trace) == (2, x_init_new, alpha, beta)
        @test choices[1 => :x] == x1
        @test choices[2 => :x] == x2
        expected_score = (logpdf(normal, x1, x_init_new * alpha + beta, std)
                          + logpdf(normal, x2, x1 * alpha + beta, std))
        @test isapprox(get_score(trace), expected_score)
        expected_weight = (logpdf(normal, x1, x_init_new * alpha + beta, std)
                           - logpdf(normal, x1, x_init * alpha + beta, std))
        @test isapprox(weight, expected_weight)
        retval = get_retval(trace)
        @test length(retval) == 2
        @test retval[1] == x1
        @test retval[2] == x2
        @test retdiff == NoChange()

        # init_changed=false, params_changed=true, change nothing
        trace = get_initial_trace()
        alpha_new = 0.5
        (trace, weight, retdiff) = regenerate(trace,
                                              (2, x_init, alpha_new, beta),
                                              (NoChange(), NoChange(), UnknownChange(), UnknownChange()), EmptySelection())
        choices = get_choices(trace)
        @test get_args(trace) == (2, x_init, alpha_new, beta)
        @test choices[1 => :x] == x1
        @test choices[2 => :x] == x2
        expected_score = (logpdf(normal, x1, x_init * alpha_new + beta, std)
                          + logpdf(normal, x2, x1 * alpha_new + beta, std))
        @test isapprox(get_score(trace), expected_score)
        expected_weight = (logpdf(normal, x1, x_init * alpha_new + beta, std)
                           - logpdf(normal, x1, x_init * alpha + beta, std)
                           + logpdf(normal, x2, x1 * alpha_new + beta, std)
                           - logpdf(normal, x2, x1 * alpha + beta, std))
        @test isapprox(weight, expected_weight)
        retval = get_retval(trace)
        @test length(retval) == 2
        @test retval[1] == x1
        @test retval[2] == x2
        @test retdiff == NoChange()
    end

    @testset "test 1 - accumulate_param_gradients!" begin

        @gen function kernel(t::Int, x_prev::Float64, (grad)(alpha::Float64), (grad)(beta::Float64))
            @param std::Float64
            x = @trace(normal(x_prev * alpha + beta, std), :x)
            return x
        end

        foo = Unfold(kernel)

        std = 1.
        set_param!(kernel, :std, std)

        x_init = 0.1
        alpha = 0.2
        beta = 0.3
        x1 = 1.1
        x2 = 1.2

        constraints = choicemap()
        constraints[1 => :x] = x1
        constraints[2 => :x] = x2
        (trace, _) = generate(foo, (2, x_init, alpha, beta), constraints)

        zero_param_grad!(kernel, :std)
        input_grads = accumulate_param_gradients!(trace, nothing)
        @test input_grads[1] == nothing # length
        @test input_grads[2] == nothing # inital state
        #@test isapprox(input_grads[3], expected_xs_grad) # alpha
        #@test isapprox(input_grads[4], expected_ys_grad) # beta
        expected_std_grad = (logpdf_grad(normal, x1, x_init * alpha + beta, std)[3]
                             + logpdf_grad(normal, x2, x1 * alpha + beta, std)[3])
        @test isapprox(get_param_grad(kernel, :std), expected_std_grad)
    end

    @gen (grad) function ker(t, (grad)(x::Float64))
        return 2.0 * x
    end

    @testset "test 2 - accumulate_param_gradients!" begin
        uf = Unfold(ker)
        tr = simulate(uf, (10, 1.0))
        k_grads = accumulate_param_gradients!(tr, [i == 10 ? 1.0 : nothing for i in 1 : 10])
        @test isapprox(k_grads[2], 1024.0)
    end

    @gen (grad) function kernel(t::Int,
            (grad)(x_prev::Float64),
            (grad)(alpha::Float64),
            (grad)(beta::Float64))
        x = @trace(normal(x_prev * alpha + beta, 1.0), :x)
        return x
    end

    @testset "choice gradients" begin
        uf = Unfold(kernel)
        alpha, beta = 1.0, 0.0
        tr = simulate(uf, (10, 0.0, alpha, beta))
        chm = get_choices(tr)

        # Test propagation over state.
        for i in 1 : 10
            sel = select(i => :x)
            grads, value_choices, grad_choices = choice_gradients(tr, sel, [nothing for _ in 1 : 10])

            # Explicit edge cases.
            grad_t = logpdf_grad(normal, chm[i => :x], i == 1 ? 0.0 : chm[(i - 1) => :x], 1.0) # if first index, init_state = 0.0
            i == 10 ? grad_tp1 = (0.0, 0.0, 0.0) : grad_tp1 = logpdf_grad(normal, chm[(i + 1) => :x], chm[i => :x], 1.0) # if last index, no more iterations.

            # Propagation should be previous choice value grad + next choice mean grad.
            @test isapprox(grad_t[1] + grad_tp1[2], grad_choices[i => :x])
        end

        # Test propagation with non-nothing ret grad.
        for i in 1 : 10
            sel = select(i => :x)
            grads, value_choices, grad_choices = choice_gradients(tr, sel, [k == i ? 1.0 : nothing for k in 1 : 10])
            grad_t = logpdf_grad(normal, chm[i => :x], i == 1 ? 0.0 : chm[(i - 1) => :x], 1.0)
            i == 10 ? grad_tp1 = (0.0, 0.0, 0.0) : grad_tp1 = logpdf_grad(normal, chm[(i + 1) => :x], chm[i => :x], 1.0)

            # Propagation should be same as above, but including accumulating the retval grad.
            @test isapprox(grad_t[1] + grad_tp1[2] + 1.0, grad_choices[i => :x])
        end
    end
end


text:
module DataType
using MLStyle
using MLStyle.MatchImpl
using MLStyle.Qualification
using MLStyle.Record: as_record
using MLStyle.ExprTools

UNREACHABLE = nothing
export @data
@nospecialize

"""
Codegen for `@data`, e.g.,

```
@data A begin              abstract type A end
    A1(Int, Int)     ≡     struct A1
                               _1 :: Int
end                            _2 :: Int
                           end

@data B{T}<:A begin        abstract type B{T} <: A end
    B1(T, Int)        ≡    struct B1{T}
end                            _1 :: T
                               _2 :: Int
                           end

@data C{T}<:B{T} begin          abstract type C{T} <: B{T} end
    C1{A, B} ::
    (A, B)=>C{Tuple{A, B}} ≡    struct C1{A, B} <: C{Tuple{A, B}}
end                                 _1 :: A
                                    _2 :: B
                                end

@data D{T}<:A begin        abstract type D{T} <: A end
    D1(a::T)          ≡    struct D1{T} <: D{T}
end                            a :: T
                           end
```

Additionally, `@data` macros setup up pattern matching for those datatypes.
```julia

@data S begin
    S1(Int)
end

@match S1(1) begin
    S1(x) => x == 1
end # => true
```
"""
macro data(typ, def_variants)
    esc(data(typ, def_variants, __source__, __module__))
end

macro data(qualifier, typ, def_variants)
    deprecate_qualifier_macro(qualifier, __source__)
    esc(data(typ, def_variants, __source__, __module__))
end

function data(typ::Any, def_variants::Expr, line::LineNumberNode, mod::Module)
    typename, tvars_of_abst = @match typ begin
        :($typename{$(a...)}) => (typename, get_type_parameters_ordered(a))
        :($typename{$(a...)} <: $b) => (typename, get_type_parameters_ordered(a))
        :($typename <: $b) => (typename, Symbol[])
        :($(typename::Symbol)) => (typename, Symbol[])
    end

    ret = Expr(:block, :(abstract type $typ end))
    impl!(ret.args, typename, tvars_of_abst, def_variants, line, mod)
    ret
end

function impl!(
    suite::Vector{Any},
    @nospecialize(abs_t),
    abs_tvars::Vector{Symbol},
    variants::Expr,
    ln::LineNumberNode,
    mod::Module
)
    abst() = isempty(abs_tvars) ? abs_t : :($abs_t{$(abs_tvars...)})
    for each in variants.args
        @switch each begin
            @case ::LineNumberNode
            ln = each
            continue

            @case Do[is_enum = false] &&  (
                 :($case{$(generic_tvars...)}::($(params...),) =>
                          $(ret_ty) where {$(constraints...)}
                  )   ||
                  :($case{$(generic_tvars...)}::($(params...),) => $ret_ty) &&
                  let constraints = []
                  end ||
                  :(
                      $case{$(generic_tvars...)}::$arg_ty =>
                          $ret_ty where {$(constraints...)}
                  ) && let params = [arg_ty]
                  end ||
                  :($case{$(generic_tvars...)}::$arg_ty => $ret_ty) &&
                  let params = [arg_ty], constraints = []
                  end ||
                  :($case::($(params...),) => $(ret_ty) where {$(constraints...)}) &&
                  let generic_tvars = []
                  end ||
                  :($case::($(params...),) => $ret_ty) &&
                  let generic_tvars = [], constraints = []
                  end ||
                  :($case::$arg_ty => $ret_ty where {$(constraints...)}) &&
                  let generic_tvars = [], params = [arg_ty]
                  end ||
                  :($case::$arg_ty => $ret_ty) &&
                  let generic_tvars = [], params = [arg_ty], constraints = []
                  end ||
                  :($case($(params...))) &&
                  let ret_ty = abst(), constraints = [], generic_tvars = abs_tvars
                  end ||
                  :($case{$generic_tvars...}::$ret_ty where {$(constraints...)}) &&
                  if error("Defining generic enum $each is invalid as Julia does not support generic value.")
                  end && let params = []
                  end ||
                  :($case{$generic_tvars...}::$ret_ty) &&
                  if error("Defining generic enum $each is invalid as Julia does not support generic value.")
                  end &&
                  let params = [], constraints = []
                  end ||
                  :($(case::Symbol)::$ret_ty) && Do[is_enum = true] &&
                  let params = [], constraints = [], generic_tvars = []
                  end ||
                  (case::Symbol) && Do[is_enum = true] &&
                  if isempty(abs_tvars) ||
                        error("Defining generic enum $case <: $(abst()) is invalid, as Julia does not support generic value.")
                  end &&
                  Do[is_enum = true] &&
                  let ret_ty = abst(), params = [], constraints = [], generic_tvars = []
                  end)

                ctor_name = is_enum ? Symbol(case, "'s constructor") : case
                expr_field_def = Expr(:block, ln)
                suite_field_def = expr_field_def.args
                field_symbols = Symbol[]
                for i in eachindex(params)
                    each = params[i]
                    @match each begin
                        a::Symbol &&
                             # TODO: check length of symbol?
                            (if Base.isuppercase(string(each)[1]) end &&
                                let field = Symbol("_$i"), ty = a end
                            ||  let field = a, ty = Any end) ||
                        :($field :: $ty) ||
                        (:(::$ty) || ty) && let field = Symbol("_$i") end =>
                            begin
                                argdecl = :($field :: $ty)
                                push!(suite_field_def, argdecl)
                                push!(field_symbols, field)
                                nothing
                            end
                        _ => error("invalid datatype field $each")
                    end
                end

                expr_struct_def = Expr(
                    :struct,
                    false,
                    Expr(:(<:),
                        Expr(:curly, ctor_name, generic_tvars...),
                        ret_ty
                    ),
                    expr_field_def
                )

                expr_infer_constructor = if isempty(constraints)
                     UNREACHABLE
                else
                    call = Expr(:call, ctor_name, suite_field_def[2:end]...)
                    type_arguments = get_type_parameters_ordered(generic_tvars)
                    fn_sig = Expr(:where, call, generic_tvars...)
                    Expr(
                        :function,
                        fn_sig,
                        Expr(:block,
                            Expr(:let,
                                Expr(:block, constraints...),
                                Expr(:call,
                                    :($ctor_name{$(type_arguments...)}),
                                    field_symbols...
                                )
                            )
                        )
                    )
                end

                expr_setup_record = if is_enum
                    err_msg = "Enumeration $case should take 0 argument, type parameter or type argument."
                    Expr(
                        :block,
                        ln,
                        :($MLStyle.pattern_uncall(::$ctor_name, _, tparams, targs, args) =
                            isempty(targs) && isempty(tparams) && isempty(args) && (
                                return $MLStyle.AbstractPatterns.literal($case)
                            ) || error($err_msg)),
                        :($MLStyle.is_enum(::$ctor_name) = true),
                        :(const $case = $reinterpret($ctor_name, nothing)),
                        :($ctor_name() = $case),
                        :($Base.show(io::$IO, ::$ctor_name) = $Base.print(io, $(string(case))))
                    )
                else
                    as_record(ctor_name, ln, mod)
                end

                push!(
                    suite,
                    expr_struct_def,
                    expr_infer_constructor,
                    expr_setup_record
                )
                continue
            @case :($case{$(_...)}) && if error("invalid enum constructor $each, use $case instead.") end
            @case _
                error("unrecognised data constructor $each.")
        end
    end
end
@specialize
end # module


text:
module Pool1

using JuMP, BARON

m = Model(with_optimizer(BARON.Optimizer))
lb = [0,3,1,2,0,0,0]
ub = [10,20,2,4,10,201,100]
@variable(m, lb[i] ≤ x[i=1:7] ≤ ub[i])

@constraints(m, begin
    x[3]^2 + x[4]^2 ≤ 12
    x[1]^2 - x[2]^2 + x[4]^2 ≥ 3
    x[1]^2 - x[2]^2 + x[4]^2 ≤ 100
    x[1] + x[2] - 5x[3] + 2x[4] ≥ 10
    x[1] + x[2] - 5x[3] + 2x[4] ≤ 20
    x[4] - x[6]^2 + 4x[7]^2 ≤ 0
end)

@objective(m, Min, 6x[1] + 16x[2] - 9x[3] - 10x[4])

optimize!(m)

end # module


text:
# ---------------------------------------------------------------------------------------------------
"""
    blendimg!(color::GMTimage, shade::GMTimage; new=false)

Blend the RGB `color` GMTimage with the `shade` intensity image (normally obtained with gdaldem)
The `new` argument determines if we return a new RGB image or update the `color` argument.

The blending method is the one explained in https://gis.stackexchange.com/questions/255537/merging-hillshade-dem-data-into-color-relief-single-geotiff-with-qgis-and-gdal/255574#255574

### Returns
A GMT RGB Image

    blendimg!(img1::GMTimage, img2::GMTimage; new=false, transparency=0.5)

Blend two 2D UInt8 or 2 RGB images using transparency.
  - **transparency** The default value, 0.5, gives equal weight to both images. 0.75 will make
    `img` weight 3/4 of the total sum, and so forth.
  - **new** If true returns a new GMTimage object, otherwise it cahnges the `img` content.

### Returns
A GMT intensity Image
"""
function blendimg!(color::GMTimage{UInt8, 3}, shade::GMTimage{UInt8, 2}; new=false)

	blend = (new) ? Array{UInt8,3}(undef, size(shade,1), size(shade,2), 3) : color.image

	n_pix = length(shade)
	if (color.layout[3] == 'B')			# Band interleaved
		for n = 1:3
			off = (n - 1) * n_pix
			@inbounds @simd for k = 1:n_pix
				t = shade.image[k] / 255
				blend[k+off] = (t < 0.5) ? round(UInt8, 2t * color.image[k+off]) : round(UInt8, (1 - 2*(1 - t) * (1 - color.image[k+off]/255)) * 255)
			end
		end
	else								# Assume Pixel interleaved
		nk = 1
		@inbounds @simd for k = 1:n_pix
			t = shade.image[k] / 255
			if (t < 0.5)
				blend[nk] = round(UInt8, 2t * color.image[nk]);		nk += 1
				blend[nk] = round(UInt8, 2t * color.image[nk]);		nk += 1
				blend[nk] = round(UInt8, 2t * color.image[nk]);		nk += 1
			else
				blend[nk] = round(UInt8, (1 - 2*(1 - t) * (1 - color.image[nk]/255)) * 255);	nk += 1
				blend[nk] = round(UInt8, (1 - 2*(1 - t) * (1 - color.image[nk]/255)) * 255);	nk += 1
				blend[nk] = round(UInt8, (1 - 2*(1 - t) * (1 - color.image[nk]/255)) * 255);	nk += 1
			end
		end
	end
	return (new) ? mat2img(blend, color) : color
end

function blendimg!(img1::GMTimage, img2::GMTimage; new=false, transparency=0.5)
	# This method blends two UInt8 images with transparency
	@assert eltype(img1) == eltype(img2)
	@assert length(img1) == length(img2)
	same_layout = (img1.layout[1:2] == img2.layout[1:2])
	blend = (new) ? Array{eltype(img1), ndims(img1)}(undef, size(img1)) : img1.image
	t, o = transparency, 1. - transparency
	if (same_layout)
		@inbounds @simd for k = 1:length(img1)
			blend[k] = round(UInt8, t * img1.image[k] + o * img2.image[k])
		end
	else
		#(size(img1,3) == 1) && error("Sorry, blending RGB images of different mem layouts is not yet implemented")
		flip, transp = img1.layout[1] != img2.layout[1], img1.layout[2] != img2.layout[2]
		if     (flip && !transp)  blend = reverse(img2.image, dims=1)
		elseif (!flip && transp)  blend = collect(img2.image')
		else                      blend = reverse(img2.image', dims=1)
		end
		@inbounds @simd for k = 1:length(img1)
			blend[k] = round(UInt8, t * img1.image[k] + o * blend[k])
		end
		(!new) && (img1.image = blend)
	end
	return (new) ? mat2img(blend, img1) : img1
end

# ---------------------------------------------------------------------------------------------------
"""
    gammacorrection(I::GMTimage, gamma; contrast=[0.0, 1.0], brightness=[0.0, 1.0])

Apply a gamma correction to a 2D (intensity) GMTimage using the exponent `gamma`.
Optionally set also `contrast` and/or `brightness`

### Returns
A GMT intensity Image
"""
function gammacorrection(I::GMTimage, gamma; contrast=[0.0, 1.0], brightness=[0.0, 1.0])

	@assert 0.0 <= contrast[1] < 1.0;	@assert 0.0 < contrast[2] <= 1.0;	@assert contrast[2] > contrast[1]
	@assert 0.0 <= brightness[1] < 1.0;	@assert 0.0 < brightness[2] <= 1.0;	@assert brightness[2] > brightness[1]
	contrast_min, contrast_max, brightness_min, brightness_max = 0.0, 1.0, 0.0, 1.0
	lut = (eltype(I) == UInt8) ? linspace(0.0, 1, 256) : linspace(0.0, 1, 65536)
	lut = max.(contrast[1], min.(contrast[2], lut))
	lut = ((lut .- contrast[1]) ./ (contrast[2] - contrast[1])) .^ gamma;
	lut = lut .* (brightness[2] - brightness[1]) .+ brightness[1];		# If brightness[1] != 0 || brightness[2] != 1
	lut = (eltype(I) == UInt8) ? round.(UInt8, lut .* 255) : round.(UInt16, lut .* 65536)
	intlut(I, lut)
end

"""
    intlut(I, lut)

Creates an array containing new values of `I` based on the lookup table, `lut`. `I` can be a GMTimage or an uint matrix.
The types of `I` and `lut` must be the same and the number of elements of `lut` is eaqual to intmax of that type.
E.g. if eltype(lut) == UInt8 then it must contain 256 elements.

### Returns
An object of the same type as I
"""
function intlut(I, lut)
	@assert eltype(I) == eltype(lut)
	mat = Array{eltype(I), ndims(I)}(undef, size(I))
	@inbounds for n = 1:length(I)
		mat[n] = lut[I[n]+1]
	end
	return (isa(I, GMTimage)) ? mat2img(mat, I) : mat
end

# ---------------------------------------------------------------------------------------------------
"""
    texture_img(G::GMTgrid; detail=1.0, contrast=2.0, uint16=false)

Compute the Texture Shading calling functions from the software from Leland Brown at
http://www.textureshading.com/Home.html

  - **detail** is the amount of texture detail. Lower values of detail retain more elevation information,
    giving more sense of the overall, large structures and elevation trends in the terrain, at the expense
	of fine texture detail. Higher detail enhances the texture but gives an overall "flatter" general appearance,
	with elevation changes and large structure less apparent.
  - **contrast** is a parameter called “vertical enhancement.” Higher numbers increase contrast in the midtones,
    but may lose detail in the lightest and darkest features. Lower numbers highlight only the sharpest ridges
	and deepest canyons but reduce contrast overall.
  - **uint16** controls if output is a UIn16 or a UInt8 image (the dafault). Note that the original code writes
    only UInt16 images bur if we want to combine this with with the hillshade computed with `gdaldem`, a UInt8
	image is more handy.

### Returns
A UInt8 (or 16) GMT Image
"""
function texture_img(G::GMTgrid; detail=1.0, contrast=2.0, uint16=false, intensity=false)
	if (G.layout[2] == 'C')  texture = reverse(G.z', dims=1)
	else                     texture = deepcopy(G.z)
	end
	terrain_filter(texture, detail, size(G,1), size(G,2), G.inc[1], G.inc[2], 0)
	(startswith(G.proj4, "+proj=merc")) && fix_mercator(texture, detail, size(G,1), size(G,2), G.range[3], G.range[4])
	(intensity) && (uint16 = true)
	terrain_image_data(texture, contrast, size(G,1), size(G,2), 0.0, (uint16) ? 65535.0 : 255.0)
	if (intensity)
		texture = texture ./ 65535 .* 2 .- 1
		Go = mat2grid(texture, G)
		Go.range[5:6] .= extrema(Go.z)
	else
		mat = (uint16) ? round.(UInt16, texture) : round.(UInt8, texture)
		Go = mat2img(mat, hdr=grid2pix(G), proj4=G.proj4, wkt=G.wkt, noconv=true, layout=G.layout*"a")
		Go.range[5:6] .= extrema(Go.image)
	end
	Go
end


text:
import TDA
import Clustering
using Test

@testset "Clustering" begin

    # test knee point selection
    x = [0.470711, 0.215764, 0.17828, 0.085663, 0.0557541, 0.0279894, 0.0288017, 0.0178112, 0.0176353]
    @test last(TDA.knee(x)) == 4

    # test elbow selection method for clustering
    clusterfn = Clustering.kmeans
    clusterselectionfn = TDA.elbow
    Y = [-0.0790049 -0.10533 -0.16009 -0.17787 -0.199888 -0.235282 -0.244151 -0.283751 -0.25926 -0.214562 -0.2007 -0.176357 -0.115027 -0.0928281 -0.0955865;
          0.503826 0.510191 0.48794 0.488375 0.457246 0.47511 0.437912 0.436693 -0.411324 -0.418725 -0.405257 -0.439537 -0.463596 -0.438397 -0.457734]
    @test all(sort!(clusterselectionfn(clusterfn, Y)) .== vcat(fill(1,8), fill(2,7)))
    @test all(sort!(clusterselectionfn(clusterfn, Y, maxk=4)) .== vcat(fill(1,8), fill(2,7)))

end


text:
"""
    _master_download(master_path::String, master_url::String)

Downloads (and unzips) a master table from HEASARC given its `url`
and a destination `path`
"""
function _master_download(master_path::String, master_url::String)
    if !isdir(dirname(master_path))
        mkpath(dirname(master_path))
    end

    @info "Downloading latest master catalog"
    Base.download(master_url, master_path)

    # Windows (used to) unzip .gz during download, unzip now if Linux
    unzip!(master_path)
end

"""
    _type_master_df!(master_df)

Slightly janky way to strongly type columns in the master table, this
needs to be done to ensure the `.feather` file is saved/read correctly

TODO: Make this less... stupid
"""
function _type_master_df!(master_df)
    nt_converter = (name=String, ra=Float64, dec=Float64, lii=Float64, bii=Float64, roll_angle=Float64,
        time=Dates.DateTime, end_time=Dates.DateTime, obsid=String, exposure=Float64, exposure_a=Float64,
        exposure_b=Float64, ontime_a=Float64, ontime_b=Float64, observation_mode=String, instrument_mode=String,
        spacecraft_mode=String, slew_mode=String, time_awarded=Float64, num_fpm=Int64,
        processing_status=String, processing_date=Dates.DateTime,public_date=Dates.DateTime,
        processing_version=String, num_processed=Int64, caldb_version=String, software_version=String,
        prnb=String, abstract=String, subject_category=String, category_code=Int64, priority=String,
        country=String, data_gap=Meta.parse, nupsdout=Meta.parse, solar_activity=String, coordinated=String,
        issue_flag=Meta.parse, comments=String, satus=String, pi_lname=String, pi_fname=String,
        cycle=Int64, obs_type=String, title=String, remarks=Union{String,Missing})

    for (name, coltype) in pairs(nt_converter)
        try
            if coltype in [String, Dates.DateTime, Union{String,Missing}]
                master_df[name] = Array{coltype,1}(master_df[name])
            else
                master_df[name] = parse.(coltype, master_df[name])
            end
        catch e
            if !in(typeof(e), [KeyError, ArgumentError])
                @warn "$name - $coltype - $e"
            end
        end
    end

    return master_df
end

"""
    _master_read_tdat(master_path::String)

Reads a raw `.tdat` table from HEASARC mastertable archives,
parses the ASCII data, finds and stores column names, cleans punctuation,
converts to `DataFrame`, strongly types the columsn, and finally returns
cleaned table as `DataFrame`
"""
function _master_read_tdat(master_path::String)
    master_ascii = readdlm(master_path, '\n')

    data_start = Int(findfirst(master_ascii .== "<DATA>")[1] + 1)
    data_end   = Int(findfirst(master_ascii .== "<END>")[1] - 1)
    keys_line  = Int(findfirst(master_ascii .== "# Data Format Specification")[1] + 2)
    field_line = Int(findfirst(master_ascii .== "# Table Parameters")[1] + 2)

    # Key names are given on the keys_line, split and make into symbols for use later
    key_names = Symbol.(split(master_ascii[keys_line][11:end])) # 11:end to remove 'line[1] = '
    no_cols   = length(key_names)
    key_obsid = findfirst(key_names .== :obsid)[1]
    key_archv = findfirst(key_names .== :processing_status)

    key_types = [line[3] for line in split.(master_ascii[field_line:field_line+no_cols-1], " ")]

    master_ascii_data = master_ascii[data_start:data_end]

    master_df = DataFrame(zeros(1, no_cols), key_names)

    deleterows!(master_df, 1) # Remove row, only made to get column names

    for (row_i, row) in enumerate(master_ascii_data)
        obs_values = split(row, "|")[1:end - 1] # Split row by | delims

        if length(obs_values) != no_cols # Some rows don't have the proper no. of columns, skip them
            @warn "Skipped row $row_i due to malformed columns, ObsID: $(obs_values[key_obsid])"
            continue
        end

        df_tmp = DataFrame()

        for (itr, key) in enumerate(key_names) # Create DataFrame of key and val for row
            cleaned = replace(obs_values[itr], "," => ".. ") # Remove some punctuation, screw with CSV
            cleaned = replace(cleaned, ";" => ".. ")

            if cleaned != ""
                if key in [:time, :end_time, :processing_date, :public_date]
                    cleaned = _mjd2datetime(Meta.parse(obs_values[itr]))
                end
            else
                cleaned = missing
            end

            df_tmp[key] = cleaned
        end

        master_df = [master_df; df_tmp] # Concat
    end

    sort!(master_df, :name)

    _type_master_df!(master_df)

    return master_df
end

"""
    master(mission::Mission; update=false)

Reads in a previously created `.feather` master table for a specific `mission_name`
using a path provided by `_mission_master_url(mission))`
"""
function master_base(mission::Mission; update=false, reload_cache=false)
    path_jaxtam         = mission_paths(mission)[:jaxtam]
    path_master_tdat    = joinpath(path_jaxtam, "master.tdat")
    path_master_feather = joinpath(path_jaxtam, "master.feather")

    if (!isfile(path_master_tdat) && !isfile(path_master_tdat)) || update
        _master_download(path_master_tdat, _mission_master_url(mission))
    end

    if isfile(path_master_feather) && !update
        @info "Loading $path_master_feather"
        master_data = Feather.read(path_master_feather)
    elseif isfile(path_master_tdat)
        @info "Loading $(path_master_tdat)"
        master_data = _master_read_tdat(path_master_tdat)
        @info "Saving $path_master_feather"
        Feather.write(path_master_feather, master_data)
    end

    # Convert Primitive's to Dates.DateTime so they don't get returned as an Arrow timestamp
    for (col_key, col_data) in DataFrames.eachcol(master_data, true)
        col_type = typeof(col_data) # 'col' is a Pair of (:key, col_data)
        if col_type <: Primitive{T} where T <: Timestamp
            master_data[col_key] = convert(Array{Dates.DateTime,1}, master_data[col_key])
        end
    end

    if isdefined(JAXTAM, Symbol(mission, "_master_df")) && update
        if reload_cache || update
            master(mission; cache=true, reload_cache=true, update=false) # Reload cache if master is updated
        end
    end

    return master_data
end

"""
    _add_append_publicity!(append_df::DataFrames.DataFrame, master_df::DataFrames.DataFrame)

Appends column of `Union{Bool,Missing}`, true if `public_date <=`now()`
"""
function _add_append_publicity!(mission::Mission, append_df::DataFrames.DataFrame, master_df::DataFrames.DataFrame)
    n  = Dates.now()
    pd = Array{DateTime,1}(master_df[:public_date])

    @info("Checking publicity")
    append_df[:publicity] = map(t->n>t, pd)

    return append_df
end

function _add_append_logged_vars!(mission::Mission, append_df::DataFrames.DataFrame, master_df::DataFrames.DataFrame)
    obs_count = size(append_df, 1)

    append_logged        = falses(obs_count)
    append_downloaded    = falses(obs_count)
    append_error         = falses(obs_count)
    append_error_stage   = Array{String,1}(undef, obs_count)
    append_report_exists = falses(obs_count)
    append_report_path   = Array{String,1}(undef, obs_count)

    full_e_range = _mission_good_e_range(mission)

    lpad_digits = floor(Int, log10(obs_count))

    @info("Looping through logs")
    for (i, obs_row) in enumerate(DataFrames.eachrow(master_df))
        print("\r\t Log $(lpad(string(i), lpad_digits, "0"))/$obs_count")
        log_path         = _log_path(mission, obs_row)
        log_exists       = isfile(log_path)
        append_logged[i] = log_exists

        append_error_stage[i] = ""
        append_report_path[i] = ""

        if log_exists
            log_contents = _log_query(mission, master_df[i, :])

            if haskey(log_contents, "meta")
                append_downloaded[i] = haskey(log_contents["meta"], :downloaded) ? log_contents["meta"][:downloaded] : false
            end

            if haskey(log_contents, "errors")
                append_error[i] = true
                append_error_stage[i] = join([string(k) for k in keys(log_contents["errors"])], ", ")
            end

            if haskey(log_contents, "web")
                if length(log_contents["web"]) > 0
                    append_report_exists[i] = true
                    if haskey(log_contents["web"], full_e_range)
                        append_report_path[i] = log_contents["web"][full_e_range]
                    else
                        append_report_path[i] = first(log_contents["web"])[2]  # First value in dict
                    end
                end
            end
        else
            continue
        end
    end

    print(" - done\n")

    append_df[:logged]        = append_logged
    append_df[:downloaded]    = append_downloaded
    append_df[:error]         = append_error
    append_df[:error_stage]   = append_error_stage
    append_df[:report_exists] = append_report_exists
    append_df[:report_path]   = append_report_path
end

"""
    _append_gen(mission, master_df)

Runs all the `_add_append` functions, returns the full `append_df`
"""
function _append_gen(mission::Mission, master_df::DataFrames.DataFrame)
    append_df = DataFrame(obsid=master_df[:obsid])

    @info("Generating append table, this may take some time")

    _add_append_publicity!(mission, append_df, master_df)
    _add_append_logged_vars!(mission, append_df, master_df)

    return append_df
end

function master_append(mission::Mission; update=false)
    path_jaxtam         = mission_paths(mission)[:jaxtam]
    path_append_feather = joinpath(path_jaxtam, "append.feather")

    if !isfile(path_append_feather) || update
        append_df = _append_gen(mission, master_base(mission))
        Feather.write(path_append_feather, append_df)

        if isdefined(JAXTAM, Symbol(mission, "_master_df"))
            master(mission; cache=true, reload_cache=true, update=false) # Reload cache if append is updated
        end
        return append_df
    else
        @info "Loading $path_append_feather"
        return Feather.read(path_append_feather)
    end
end

"""
    master(::Mission; cache::Bool, reload_cache::Bool, update::Bool)

Wrapper function which deals with downloading and updating both the `master` and `append` files, as well as
the caching functionality

Will download and set up the master and append files when first called

Subsequent calls will use the downloaded files

If you want to update the master table, simply call `master(::Mission, update=true)`
"""
function master(mission::Mission; cache::Bool=true, reload_cache::Bool=false, update::Bool=false)
    master_df_var = Symbol(_mission_name(mission), "_master_df")

    if update
        master_base(mission; update=true)
        master_append(mission; update=true)
    end

    if cache
        if reload_cache
            @info "Reloading master_append cache"
        end

        if isdefined(JAXTAM, master_df_var) && !reload_cache
            @assert typeof(getproperty(JAXTAM, master_df_var)) == DataFrames.DataFrame
            return getproperty(JAXTAM, master_df_var)
        else cache
            master_df = master_base(mission)
            append_df = master_append(mission)

            # If reloading cache after master update, then current append_df may not be the correct size
            # for some odd edge cases, this is a hacky fix
            if size(master_df, 1) != size(append_df, 1)
                append_df = master_append(mission; update=true)
            end

            master_df_a = join(master_df, append_df, on=:obsid)

            eval(:(global const $master_df_var = $master_df_a))

            return getproperty(JAXTAM, master_df_var)
        end
    else
        master_df = master_base(mission)
        append_df = master_append(mission)
        master_df_a = join(master_df, append_df, on=:obsid)
        return master_df_a
    end
end


text:
using DelimitedFiles

filename = "./data_files/simple_runtest_01.dat"

# Reads the execution data
data = readdlm(filename)

# Dimension of the function domain
n = 2

# Number of points in sample set
m = 3

# Splits the useful data
data_c = data[1:(end - 1), 14]
data_g = data[1:(end - 1), 15:(14 + n)]
data_xbase = data[1:(end - 1), (15 + n):(14 + 2 * n)]
data_xopt = data[1:(end - 1), (15 + 2 * n):(14 + 3 * n)]
data_Y = data[1:(end - 1), (15 + 3 * n + m):(14 + 3 * n + m + n * (m - 1))]
data_d = data[1:(end - 1), (end - 1):end ]

it = 64
println("xbase: ", data_xbase[it, :])
println("xopt : ", data_xopt[it, :])
println("Y_{1}: ", data_Y[it, 1:2] + data_xbase[it, :])
println("Y_{2}: ", data_Y[it, 3:4] + data_xbase[it, :])
println("d: ", data_d[it, :])


text:
abstract type CropMethod end
struct RedMax <: CropMethod end

function isredmax(p::T) where T <: ColorTypes.AbstractRGB
    pmax = max(p.r, p.g, p.b)
    # disallow masked and overexposed pixels, assume plants reflect better than dirt
    pmax == 1 && return false
    return p.r == pmax
end

gapfraction(pixs::AbstractArray, thresh::RedMax) = sum(isredmax(p) for p in pixs) / length(pixs)


text:
using DataDrivenDiffEq
using ModelingToolkit
using LinearAlgebra
using OrdinaryDiffEq
using Test
using DataDrivenDiffEq.Optimise

include("./basis.jl")
include("./koopman.jl")
include("./sindy.jl")
include("./isindy.jl")
include("./utils.jl")
include("./optimise.jl")


text:
using Base.Test

using Rifraf

@testset "BandedArrays" begin
    @testset "inband" begin
        m = BandedArray(Int, (13, 11), 5)
        @test inband(m, 1, 1)
        @test inband(m, 8, 1)
        @test !inband(m, 9, 1)
        @test inband(m, 1, 6)
        @test !inband(m, 1, 7)
    end

    @testset "data_row" begin
        m = BandedArray(Int, (3, 3), 1)
        @test data_row(m, 1, 1) == 2
        @test data_row(m, 2, 1) == 3
        @test data_row(m, 1, 2) == 1
        @test data_row(m, 2, 2) == 2
        @test data_row(m, 3, 2) == 3
        @test data_row(m, 2, 3) == 1
        @test data_row(m, 3, 3) == 2

        m = BandedArray(Int, (3, 5), 1)
        @test data_row(m, 1, 1) == 4
        @test data_row(m, 3, 5) == 2

        m = BandedArray(Int, (5, 3), 1)
        @test data_row(m, 1, 1) == 2
        @test data_row(m, 5, 3) == 4
    end

    @testset "row_range" begin
        m = BandedArray(Int, (3, 5), 1)
        @test row_range(m, 1) == (1, 2)
        @test row_range(m, 2) == (1, 3)

        m = BandedArray(Int, (5, 3), 1)
        @test row_range(m, 1) == (1, 4)
        @test row_range(m, 2) == (1, 5)
    end

    @testset "data_row_range" begin
        m = BandedArray(Int, (3, 5), 1)
        @test data_row_range(m, 1) == (4, 5)
        @test data_row_range(m, 2) == (3, 5)

        m = BandedArray(Int, (5, 3), 1)
        @test data_row_range(m, 1) == (2, 5)
        @test data_row_range(m, 2) == (1, 5)
    end

    @testset "sparsecol" begin
        m = BandedArray(Int, (5, 3), 1)
        m[1, 1] = 1
        @test sparsecol(m, 1) == [1, 0, 0, 0]
    end

    @testset "flip" begin
        m = BandedArray(Int, (5, 3), 1, row_padding=5, col_padding=5)
        m[1, 1] = 1
        flip!(m)
        @test m[5, 3] == 1
    end

    @testset "sym_band" begin
        m = BandedArray(Int, (3, 3), 1)
        m.data[:] = 1
        expected = ones(Int, (3, 3))
        expected[3, 1] = 0
        expected[1, 3] = 0
        @test full(m) == expected
    end

    @testset "test_wide" begin
        m = BandedArray(Int, (3, 4), 1)
        m.data[:] = 1
        expected = ones(Int, (3, 4))
        expected[1, end] = 0
        expected[end, 1] = 0
        @test full(m) == expected
    end

    @testset "wide col" begin
        m = BandedArray(Int, (3, 5), 1)
        m.data[:] = 1
        first = ones(Int, 2)
        middle = ones(Int, 3)
        last = first
        @test sparsecol(m, 1) == first
        @test sparsecol(m, 2) == middle
        @test sparsecol(m, 3) == middle
        @test sparsecol(m, 4) == middle
        @test sparsecol(m, 5) == last
    end

    @testset "tall" begin
        m = BandedArray(Int, (4, 3), 1)
        m.data[:] = 1
        expected = ones(Int, (4, 3))
        expected[1, end] = 0
        expected[end, 1] = 0
        @test full(m) == expected
    end

    @testset "tall band" begin
        m = BandedArray(Int, (5, 3), 1)
        m.data[:] = 1
        expected = ones(Int, (5, 3))
        expected[5, 1] = 0
        expected[1, 3] = 0
        @test full(m) == expected
    end

    @testset "individual setting" begin
        m = BandedArray(Int, (3, 3), 1)
        m[1, 2] = 3
        m[2, 1] = 5
        expected = zeros(Int, (3, 3))
        expected[1, 2] = 3
        expected[2, 1] = 5
        @test full(m) == expected
    end

    @testset "individual setting of entire band" begin
        m = BandedArray(Int, (3, 3), 1)
        m[1, 1] = 1
        m[2, 1] = 1
        m[1, 2] = 2
        m[2, 2] = 2
        m[3, 2] = 2
        m[2, 3] = 3
        m[3, 3] = 3
        expected = zeros(Int, (3, 3))
        expected[1:2, 1] = 1
        expected[1:3, 2] = 2
        expected[2:3, 3] = 3
        @test full(m) == expected
    end

    @testset "individual setting of entire band, wide array" begin
        m = BandedArray(Int, (3, 5), 1)
        m[1, 1] = 1
        m[1, 2] = 1
        m[1, 3] = 1
        m[3, 5] = 2
        expected = zeros(Int, (3, 5))
        expected[1, 1:3] = 1
        expected[3, 5] = 2
        @test full(m) == expected
    end
    @testset "individual setting of entire band, tall array" begin
        m = BandedArray(Int, (5, 3), 1)
        m[1, 1] = 1
        m[2, 1] = 1
        m[3, 1] = 1
        m[4, 1] = 1
        m[5, 3] = 2
        expected = zeros(Int, (5, 3))
        expected[1:4, 1] = 1
        expected[5, 3] = 2
        @test full(m) == expected
    end

    @testset "resize up within padding" begin
        m = BandedArray(Int, (3, 3), 1, row_padding=10, col_padding=10)
        old = m.data
        resize!(m, (5, 5))
        @test m.data == old
    end

    @testset "resize down" begin
        m = BandedArray(Int, (5, 5), 1)
        old = m.data
        resize!(m, (3, 3))
        @test m.data == old
    end

    @testset "resize with reallocate" begin
        m = BandedArray(Int, (3, 3), 1, row_padding=5, col_padding=5)
        @test Rifraf.row_range(m, 1) == (1, 2)
        old = m.data
        resize!(m, (5, 10))
        @test m.data != old
        @test Rifraf.row_range(m, 1) == (1, 2)
        @test Rifraf.row_range(m, 3) == (1, 4)
        @test Rifraf.row_range(m, 5) == (1, 5)
        @test Rifraf.row_range(m, 10) == (4, 5)
    end

    @testset "resize columns" begin
        m = BandedArray(Int, (3, 3), 1, row_padding=5, col_padding=10)
        @test Rifraf.row_range(m, 1) == (1, 2)
        old = m.data
        resize!(m, (5, 10))
        @test m.data == old
    end

    @testset "test get and set" begin
        default = Score(-Inf)
        m = BandedArray(Score, (13, 11), 5, default=default)
        m[1, 1] = Score(1.0)
        @test m[1, 1] == Score(1.0)
        @test m[end, 1] == default
    end

    @testset "equal_ranges" begin
        @test Rifraf.equal_ranges((3, 5), (4, 6)) == ((2, 3), (1, 2))
        @test Rifraf.equal_ranges((1, 5), (1, 2)) == ((1, 2), (1, 2))
        @test Rifraf.equal_ranges((1, 5), (4, 5)) == ((4, 5), (1, 2))
    end

end


text:
#=
Generates an array of draws from a discrete random variable with a
specified vector of probabilities.

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 2014-07-10

References
----------

http://quant-econ.net/jl/finite_markov.html?highlight=discrete_rv

TODO: as of 07/10/2014 it is not possible to define the property
      interface we see in the python version. Once issue 1974 from
      the main Julia repository is resolved, we can revisit this and
      implement that feature.
=#

"""
Generates an array of draws from a discrete random variable with
vector of probabilities given by q.

##### Fields

- `q::Vector{T<:Real}`: A vector of non-negative probabilities that sum to 1
- `Q::Vector{T<:Real}`: The cumulative sum of q
"""
type DiscreteRV{T<:Real}
    q::Vector{T}
    Q::Vector{T}
    DiscreteRV(x::Vector{T}) = new(x, cumsum(x))
end

# outer constructor so people don't have to put {T} themselves
DiscreteRV{T<:Real}(x::Vector{T}) = DiscreteRV{T}(x)

"""
Make a single draw from the discrete distribution

##### Arguments

- `d::DiscreteRV`: The `DiscreteRV` type represetning the distribution

##### Returns

- `out::Int`: One draw from the discrete distribution
"""
draw(d::DiscreteRV) = searchsortedfirst(d.Q, rand())

"""
Make multiple draws from the discrete distribution represented by a
`DiscreteRV` instance

##### Arguments

- `d::DiscreteRV`: The `DiscreteRV` type representing the distribution
- `k::Int`:

##### Returns

- `out::Vector{Int}`: `k` draws from `d`
"""
draw{T}(d::DiscreteRV{T}, k::Int) = Int[draw(d) for i=1:k]


text:
# The default `AbstractArr`, more-or-less `Array`
struct Arr{T,N} <: AbstractArr{T, N}
    a::Array{T, N}
end

const Vec{T} = Arr{T, 1}
Vec(v::Vector{T}) where {T} = Vec{T}(v)

const Mat{T} = Arr{T, 2}
Mat(m::Matrix{T}) where {T} = Mat{T}(m)

indices(a::Arr) = keys(a.a)

# Support getindex, etc, like currently exists
@propagate_inbounds getindex(a::Arr, i::Int) = a.a[i]
@propagate_inbounds getindex(a::Arr{<:Any, N}, i::Vararg{Int, N}) where {N} = a.a[i...]
@propagate_inbounds getindex(a::Arr{<:Any, N}, i::CartesianIndex{N}) where {N} = a.a[i]
@propagate_inbounds getindex(a::Arr, i::Union{Int,AbstractVector{Int},Colon}...) = Arr(a.a[i])

@propagate_inbounds setindex!(a::Arr{T}, v, i...) where {T} = (a.a[i...] = v; v)

similar(a::Assoc, lengths::Vararg{Int, N}) where {T, N} = similar(a, eltype(a), lengths...)
similar(a::Assoc, ::Type{T}, lengths::Vararg{Int, N}) where {T, N} = similar(a, T, CartesianRange(lengths))
similar(::Assoc, ::Type{T}, inds::Base.OneTo{Int}) where {T} = Vec{T}(Vector{T}(last(inds)))
similar(::Assoc, ::Type{T}, inds::CartesianRange{N}) where {T, N} = Arr{T, N}(Array{T, N}(map(last, inds)))

empty(::Assoc{T,N}, ::Type{CartesianIndex{1}}, ::Type{T}) where {T, N} = Vec{T}(Vector{T}(ntuple(i->0, Val(N))))
# No point making empty 0D, 2D, 3D, etc arrays... makes no sense!


text:
using SimJulia

@resumable function client(sim::Simulation, res::Resource, i::Int, priority::Int)
  println("$(now(sim)), client $i is waiting")
  @yield request(res, priority=priority)
  println("$(now(sim)), client $i is being served")
  @yield timeout(sim, rand())
  println("$(now(sim)), client $i has been served")
  @yield release(res)
end

@resumable function generate(sim::Simulation, res::Resource)
  for i in 1:10
    @process client(sim, res, i, 10-i)
    @yield timeout(sim, 0.5*rand())
  end
end

sim = Simulation()
res = Resource(sim, 2; level=1)
println(res)
@process generate(sim, res)
run(sim)

@resumable function my_consumer(sim::Simulation, con::Container)
  for i in 1:10
    amount = 3*rand()
    println("$(now(sim)), consumer is demanding $amount")
    @yield timeout(sim, 1.0*rand())
    get_ev = get(con, amount)
    val = @yield get_ev | timeout(sim, rand())
    if val[get_ev].state == SimJulia.processed
      println("$(now(sim)), consumer is being served, level is ", con.level)
      delay = 5.0*rand()
      @yield timeout(sim, delay)
    else
      println("$(now(sim)), consumer has timed out")
      cancel(con, get_ev)
    end
  end
end

@resumable function my_producer(sim::Simulation, con::Container)
  for i in 1:10
    amount = 2*rand()
    println("$(now(sim)), producer is offering $amount")
    @yield timeout(sim, 1.0*rand())
    @yield put(con, amount)
    level = con.level
    println("$(now(sim)), producer is being served, level is ", level)
    delay = 5.0*rand()
    @yield timeout(sim, delay)
  end
end

sim = Simulation()
con = Container(sim, 10.0; level=5.0)
@process my_consumer(sim, con)
@process my_producer(sim, con)
run(sim)


text:
"""
    SpeGuid(tolerance::Float64)

SpeGuid, Specification-Guided Safety Verification, performs over-approximated reachability analysis to compute the over-approximated output reachable set for a network.

# Problem requirement
1. Network: any depth, any activation that is monotone
2. Input: `HPolytope`
3. Output: `HPolytope`

# Return
`BasicResult`

# Method
Search and Reachability

# Property
Sound but not complete.

# Reference
[W. Xiang, H.-D. Tran, and T. T. Johnson,
"Specification-Guided Safety Verification for Feedforward Neural Networks,"
*ArXiv Preprint ArXiv:1812.06161*, 2018.](https://arxiv.org/abs/1812.06161)
"""
@with_kw struct SpeGuid
    tolerance::Float64 = 1.0
end

# This is the main function

function solve(solver::SpeGuid, problem::Problem)
    input = problem.input
    stack = Vector{Hyperrectangle}(undef, 0)
    push!(stack, input)
    #count = 1
    while !isempty(stack)
        interval = popfirst!(stack)
        reach = forward_network(solver, problem.network, interval)
        if issubset(reach, problem.output)
            continue
        else
            if get_largest_width(interval) > solver.tolerance
                sections = bisect(interval)
                for i in 1:2
                    push!(stack, sections[i])
                    #count += 1
                end
            else
                return BasicResult(:violated)
            end
        end
    end
    #println(count)
    return BasicResult(:holds)
end

function get_largest_width(input::Hyperrectangle)
    width = high(input) - low(input)
    return max(width...)
end

#interval refinement
function bisect(input::Hyperrectangle)
    lower, upper = low(input), high(input)
    width = upper - lower
    index = 0
    temp = 0
    for i in 1:length(width)
        if width[i] >= temp
            temp = width[i]
            index = i
        end
    end
    upper[index] = input.center[index]
    lower_part = Hyperrectangle(low = lower, high = upper)
    lower[index] = input.center[index]
    upper[index] = input.center[index] + input.radius[index]
    upper_part = Hyperrectangle(low = lower, high = upper)
    return (lower_part, upper_part)
end

function forward_layer(solver::SpeGuid, L::Layer, input::Hyperrectangle)
    (W, b, act) = (L.weights, L.bias, L.activation)
    center = zeros(size(W, 1))
    gamma  = zeros(size(W, 1))
    for j in 1:size(W, 1)
        node = Node(W[j,:], b[j], act)
        center[j], gamma[j] = forward_node(solver, node, input)
    end
    return Hyperrectangle(center, gamma)
end

function forward_node(solver::SpeGuid, node::Node, input::Hyperrectangle)
    output    = node.w' * input.center + node.b
    deviation = sum(abs.(node.w) .* input.radius)
    βmax = node.act(output + deviation)
    βmin = node.act(output - deviation)
    return ((βmax + βmin)/2, (βmax - βmin)/2)
end


text:
"""
    SurfaceFluxes

## Interface
  - [`surface_conditions`](@ref) computes
    - Monin-Obukhov length
    - Potential temperature flux (if not given) using Monin-Obukhov theory
    - transport fluxes using Monin-Obukhov theory
    - friction velocity/temperature scale/tracer scales
    - exchange coefficients

## References
 - [Nishizawa2018](@cite)

"""
module SurfaceFluxes

using NonlinearSolvers
using KernelAbstractions: @print

using ..Thermodynamics
using DocStringExtensions
using CLIMAParameters: AbstractEarthParameterSet
using CLIMAParameters.Planet: molmass_ratio, grav
using CLIMAParameters.SubgridScale: von_karman_const
using StaticArrays

include("UniversalFunctions.jl")
using .UniversalFunctions
const UF = UniversalFunctions

abstract type SurfaceFluxesModel end

struct FVScheme end
struct DGScheme end

export surface_conditions,
    exchange_coefficients, recover_profile, monin_obukhov_length

"""
    SurfaceFluxConditions{FT}

Surface flux conditions, returned from `surface_conditions`.

# Fields

$(DocStringExtensions.FIELDS)
"""
struct SurfaceFluxConditions{FT, VFT}
    L_MO::FT
    wθ_flux_star::FT
    flux::VFT
    x_star::VFT
    C_exchange::VFT
end

function Base.show(io::IO, sfc::SurfaceFluxConditions)
    println(io, "----------------------- SurfaceFluxConditions")
    println(io, "L_MO           = ", sfc.L_MO)
    println(io, "wθ_flux_star = ", sfc.wθ_flux_star)
    println(io, "flux           = ", sfc.flux)
    println(io, "x_star         = ", sfc.x_star)
    println(io, "C_exchange     = ", sfc.C_exchange)
    println(io, "-----------------------")
end

function surface_fluxes_f!(F, x, nt)
    param_set = nt.param_set
    wθ_flux_star = nt.wθ_flux_star
    z_in = nt.z_in
    z_0 = length(nt.z_0) > 1 ? Tuple(nt.z_0) : nt.z_0
    x_in = Tuple(nt.x_in)
    x_s = Tuple(nt.x_s)
    n_vars = nt.n_vars
    scheme = nt.scheme
    universal_func = nt.universal_func
    θ_scale = nt.θ_scale

    x_tup = Tuple(x)

    u_star, θ_star = x_tup[2], x_tup[3]
    if wθ_flux_star == nothing
        wθ_surf_flux = -u_star * θ_star
    else
        wθ_surf_flux = wθ_flux_star
    end
    L_MO = monin_obukhov_length(param_set, u_star, θ_scale, wθ_surf_flux)
    uf = universal_func(param_set, L_MO)
    F_nt = ntuple(Val(n_vars + 1)) do i
        if i == 1
            F_i =
                x_tup[1] - monin_obukhov_length(
                    param_set,
                    u_star,
                    θ_scale,
                    wθ_surf_flux,
                )
        else
            ϕ = x_tup[i]
            transport = i - 1 == 1 ? MomentumTransport() : HeatTransport()
            F_i =
                ϕ - compute_physical_scale(
                    uf,
                    z_in,
                    length(nt.z_0) > 1 ? z_0[i - 1] : z_0,
                    x_in[i - 1],
                    x_s[i - 1],
                    transport,
                    scheme,
                )
        end
        F_i
    end
    F .= F_nt
end

"""
    surface_conditions
Surface conditions given
 - `MO_param_guess` initial guess for solution (`L_MO, u_star, θ_star, ϕ_star, ...`)
 - `x_in` Inner values for state variable array `x`. They correspond to volume averages in FV and
    nodal point values in DG
 - `x_s` surface values for state variable array `x`
 - `z_0` roughness lengths for state variable array `x`
 - `θ_scale` virtual dry static energy (i.e., basic potential temperature)
 - `z_in` Input height for the similarity functions. It is Δz for FV and the height
    of the first nodal point for DG
 - `wθ_flux_star` potential temperature flux (optional)
  - `universal_func` family of flux-gradient universal functions used (optional)

If `wθ_flux_star` is not given, then it is computed by iteration
of equations 3, 17, and 18 in Nishizawa2018.
"""
function surface_conditions(
    param_set::AbstractEarthParameterSet,
    MO_param_guess::AbstractVector,
    x_in::AbstractVector,
    x_s::AbstractVector,
    z_0::Union{AbstractVector, FT},
    θ_scale::FT,
    z_in::FT,
    scheme,
    wθ_flux_star::Union{Nothing, FT} = nothing,
    universal_func::Union{Nothing, F} = Businger,
) where {FT <: AbstractFloat, AbstractEarthParameterSet, F}

    n_vars = length(MO_param_guess) - 1
    @assert length(MO_param_guess) == n_vars + 1
    @assert length(x_in) == n_vars
    @assert length(x_s) == n_vars
    local sol

    args = (;
        param_set,
        wθ_flux_star,
        z_in,
        z_0,
        n_vars,
        x_in,
        x_s,
        scheme,
        θ_scale,
        MO_param_guess,
        universal_func,
    )

    # Define closure over args
    f!(F, x_all) = surface_fluxes_f!(F, x_all, args)

    nls = NewtonsMethodAD(f!, MO_param_guess)
    # sol = solve!(nls, CompactSolution(), ResidualTolerance(FT(10)), 1)
    sol = solve!(nls, CompactSolution())

    root_tup = Tuple(sol.root)
    if sol.converged
        L_MO, x_star = root_tup[1], SVector(root_tup[2:end])
        u_star, θ_star = x_star[1], x_star[2]
    else
        @print("Warning: Unconverged Surface Fluxes\n")
        L_MO, x_star = root_tup[1], SVector(root_tup[2:end])
        u_star, θ_star = x_star[1], x_star[2]
    end

    _grav::FT = grav(param_set)
    _von_karman_const::FT = von_karman_const(param_set)
    wθ_flux_star = -u_star^3 * θ_scale / (_von_karman_const * _grav * L_MO)
    flux = -u_star * x_star

    C_exchange = get_flux_coefficients(
        param_set,
        z_in,
        x_star,
        x_s,
        L_MO,
        z_0,
        scheme,
        universal_func,
    )

    VFT = typeof(flux)
    return SurfaceFluxConditions{FT, VFT}(
        L_MO,
        wθ_flux_star,
        flux,
        x_star,
        C_exchange,
    )
end

"""
    get_energy_flux(surf_conds::SurfaceFluxConditions{FT, VFT})

Returns the potential temperature and q_tot fluxes needed to
compute the energy flux.
"""
function get_energy_flux(surf_conds::SurfaceFluxConditions)
    θ_flux = surf_conds.flux[2]
    q_tot_flux = surf_conds.flux[3]
    return θ_flux, q_tot_flux
end

"""
    compute_physical_scale(uf, z_in, z_0, x_in, x_s, transport, ::FVScheme)

Returns u_star, θ_star, ... given equations (17), (18) for FV and (8), (9) for DG
from Nishizawa & Kitamura (2018).

## Arguments
 - `uf` universal function family
 - `z_in` Input height for the similarity functions. It is Δz for FV and the height
    of the first nodal point for DG
 - `z_0` roughness lengths for state variable array `x`
 - `x_in` Inner values for state variable array `x`. They correspond to volume averages in FV and
    nodal point values in DG
 - `x_s` surface values for state variable array `x`
 - `transport` Type of transport (momentum or heat)

"""
function compute_physical_scale(
    uf::AbstractUniversalFunction{FT},
    z_in,
    z_0,
    x_in,
    x_s,
    transport,
    ::FVScheme,
) where {FT}
    _von_karman_const::FT = von_karman_const(uf.param_set)
    _π_group = FT(UF.π_group(uf, transport))
    _π_group⁻¹ = (1 / _π_group)
    R_z0 = 1 - z_0 / z_in
    temp1 = log(z_in / z_0)
    temp2 = -Psi(uf, z_in / uf.L, transport)
    temp3 = z_0 / z_in * Psi(uf, z_0 / uf.L, transport)
    temp4 = R_z0 * (psi(uf, z_0 / uf.L, transport) - 1)
    Σterms = temp1 + temp2 + temp3 + temp4
    return _π_group⁻¹ * _von_karman_const / Σterms * (x_in - x_s)
end

function compute_physical_scale(
    uf::AbstractUniversalFunction{FT},
    z_in,
    z_0,
    x_in,
    x_s,
    transport,
    ::DGScheme,
) where {FT}
    _von_karman_const::FT = von_karman_const(uf.param_set)
    _π_group = FT(UF.π_group(uf, transport))
    _π_group⁻¹ = (1 / _π_group)
    temp1 = log(z_in / z_0)
    temp2 = -psi(uf, z_in / uf.L, transport)
    temp3 = psi(uf, z_0 / uf.L, transport)
    Σterms = temp1 + temp2 + temp3
    return _π_group⁻¹ * _von_karman_const / Σterms * (x_in - x_s)
end

"""
    recover_profile(z, x_star, x_s, z_0, L_MO, transport, ::DGScheme, uf)

Recover vertical profiles u(z), θ(z), ... using equations (4) and (5)
for DG and (12), (13) for FV from Nishizawa & Kitamura (2018).

## Arguments
 - `z` Input height for evaluation
 - `x_star` Surface layer scales for state variable array `x`.
 - `x_s` surface values for state variable array `x`
 - `z_0` roughness lengths for state variable array `x`
 - `L_MO` Obukhov length
 - `transport` Type of transport (momentum or heat)
 - `uf` universal function family

"""
function recover_profile(
    param_set::AbstractEarthParameterSet,
    z,
    x_star,
    x_s,
    z_0::Union{AbstractVector, FT},
    L_MO,
    transport,
    ::DGScheme,
    universal_func = Businger,
) where {FT}
    uf = universal_func(param_set, L_MO)
    _von_karman_const::FT = von_karman_const(param_set)
    _π_group = FT(UF.π_group(uf, transport))
    temp1 = log(z / z_0)
    temp2 = -psi(uf, z / uf.L, transport)
    temp3 = psi(uf, z_0 / uf.L, transport)
    Σterms = temp1 + temp2 + temp3
    return _π_group * x_star * Σterms / _von_karman_const + x_s
end

function recover_profile(
    param_set::AbstractEarthParameterSet,
    z,
    x_star,
    x_s,
    z_0::Union{AbstractVector, FT},
    L_MO,
    transport,
    ::FVScheme,
    universal_func = Businger,
) where {FT}
    uf = universal_func(param_set, L_MO)
    _von_karman_const::FT = von_karman_const(param_set)
    _π_group = FT(UF.π_group(uf, transport))
    R_z0 = 1 - z_0 / z
    temp1 = log(z / z_0)
    temp2 = -Psi(uf, z / uf.L, transport)
    temp3 = z_0 / z * Psi(uf, z_0 / uf.L, transport)
    temp4 = R_z0 * (psi(uf, z_0 / uf.L, transport) - 1)
    Σterms = temp1 + temp2 + temp3 + temp4
    return _π_group * x_star * Σterms / _von_karman_const + x_s
end

### Generic terms

"""
    monin_obukhov_length(param_set, u_star, θ_scale, wθ_surf_flux)

Returns the Monin-Obukhov length.
"""
function monin_obukhov_length(
    param_set::AbstractEarthParameterSet,
    u_star,
    θ_scale,
    wθ_surf_flux,
)
    FT = typeof(u_star)
    _grav::FT = grav(param_set)
    _von_karman_const::FT = von_karman_const(param_set)
    return -u_star^3 * θ_scale /
           (_von_karman_const * _grav * wθ_surf_flux + eps(FT))
end

monin_obukhov_length(sfc::SurfaceFluxConditions) = sfc.L_MO

"""
    exchange_coefficients(
        param_set,
        z,
        F_exchange,
        x_star::VFT,
        θ_scale,
        L_MO,
        universal_func,
    )

Computes exchange transfer coefficients
  - `K_exchange` exchange coefficients
"""
function exchange_coefficients(
    param_set,
    z,
    F_exchange,
    x_star::VFT,
    θ_scale,
    L_MO,
    universal_func,
) where {VFT}
    N = length(F_exchange)
    FT = typeof(z)
    _von_karman_const::FT = von_karman_const(param_set)
    uf = universal_func(param_set, L_MO)
    x_star_tup = Tuple(x_star)
    K_exchange = similar(x_star)
    F_exchange_tup = Tuple(F_exchange)
    K_exchange .= ntuple(Val(length(x_star))) do i
        transport = i == 1 ? MomentumTransport() : HeatTransport()
        phi_t = phi(uf, z / L_MO, transport)
        _π_group = FT(UF.π_group(uf, transport))
        num = -F_exchange_tup[i] * _von_karman_const * z
        den = _π_group * (x_star_tup[i] * phi_t)
        K_exch = num / den # Eq. 19 in
    end
    return K_exchange
end

"""
    get_flux_coefficients(
        param_set,
        z,
        x_star::VFT,
        θ_scale,
        L_MO,
        z0,
        universal_func,
    )

Returns the exchange coefficients for bulk transfer formulas associated
with the frictional parameters `x_star` and the surface layer similarity
profiles. Taken from equations (36) and (37) of Byun (1990).
"""
function get_flux_coefficients(
    param_set,
    z_in,
    x_star::VFT,
    x_s,
    L_MO,
    z0::Union{AbstractVector, FT},
    scheme,
    universal_func,
) where {VFT, FT}
    N = length(x_star)
    z0_tup = Tuple(z0)
    x_star_tup = Tuple(x_star)
    x_s_tup = Tuple(x_s)
    u_in = recover_profile(
        param_set,
        z_in,
        x_star_tup[1],
        x_s_tup[1],
        (length(z0) > 1 ? z0_tup[1] : z0),
        L_MO,
        MomentumTransport(),
        scheme,
        universal_func,
    )
    C = similar(x_star)
    C .= ntuple(Val(length(x_star))) do i
        if i == 1
            C_i = x_star[i]^2 / (u_in - x_s_tup[i])^2
        else
            ϕ_in = recover_profile(
                param_set,
                z_in,
                x_star_tup[i],
                x_s_tup[i],
                (length(z0) > 1 ? z0_tup[i] : z0),
                L_MO,
                HeatTransport(),
                scheme,
                universal_func,
            )
            C_i =
                x_star[1] * x_star[i] / (u_in - x_s_tup[1]) /
                (ϕ_in - x_s_tup[i])
        end
        C_i
    end
    return C
end

end # SurfaceFluxes module


text:
# Very simple Scoring systems created similar to the select-regress-round technique in:
#  Jung et al, "Simple Rules for Complex Decisions" https://arxiv.org/pdf/1702.04690.pdf
# However, we do LassoPath-round to get multiple scoring systems from a single run of
# Lasso and we use the Lasso path on all features rather than forward regression to
# select features.
using Lasso

abstract type LearntModel end
abstract type Classifier <: LearntModel end

struct LinearIntegerModel <: Classifier
    L::Int                  # Max beta value used when creating the scoring system
    β::AbstractVector{Int}
    treshold::Int # If sum of betas is higher that treshold we predict the class
end
marginconvert(v::Number) = (v < 0) ? -1 : 1
zerooneconvert(v::Number) = (v < 0) ? 0 : 1
extremascores(β::AbstractVector{N}) where {N<:Number} =
    (sum(filter(v -> v < 0, β)), sum(filter(v -> v > 0, β)))
function classify(ss::LinearIntegerModel, X::AbstractMatrix{N}) where {N<:Number}
    zerooneconvert.(ss.treshold .+ X * ss.β)
end

using StatsBase
StatsBase.predict(m::LinearIntegerModel, X::AbstractMatrix) = classify(m, X)

# Calculate integer score in range from -L to L, given beta value and absolute max beta.
rescale(beta, absmaxbeta, L) = round(Int, L*beta/absmaxbeta)

function rescale_lasso_coefficients(betas::AbstractVector{Float64}, beta0 = 0; L::Int = 10)
    amb = maximum(abs.(betas))
    rsbetas = rescale.(betas, amb, L)
    return rsbetas, rescale(beta0, amb, L)
end

# Returns a set of zero or more models that have been optimized given
# a dataset.
abstract type ModelOptimizer end

# A Lasso-based Linear Integer Model optimizer:
struct LassoLIMOptimizer <: ModelOptimizer
    withintercept::Bool
    maxnummodels::Int
    L::Int
    LassoLIMOptimizer(; L::Int = 10, withintercept = true, maxnummodels = 30) =
        new(withintercept, maxnummodels, L)
end

function optimize(::Type{LinearIntegerModel}, o::LassoLIMOptimizer,
    X::AbstractMatrix, y::AbstractVector)
    lp = fit(LassoPath, X, y, Binomial(), LogitLink(); intercept = o.withintercept)
    nmodels = min(o.maxnummodels, size(lp.coefs, 2))
    wi = o.withintercept
    coefs = unique(map(i ->
        rescale_lasso_coefficients(lp.coefs[:, i], (wi ? lp.b0[i] : 0); L = o.L),
        2:nmodels))
    LinearIntegerModel[LinearIntegerModel(o.L, c...) for c in coefs]
end

using ROC

function rankbyauc(models::AbstractVector{M}, X, y) where {M<:LearntModel}
    preds = map(m -> predict(m, X), models)
    rocs = map(pr -> roc(pr, y), preds)
    aucs = map(AUC, rocs)
    p = sortperm(aucs, rev=true)
    return models[p], aucs[p], preds[p]
end

using SparseArrays
N, P = 10_000, 5_000
realcoefs = zeros(Int, P)
realcoefs[1] = 4
realcoefs[42] = -2
realcoefs[244] = -8
@time sparsecoefs = sparsevec(realcoefs)
X = rand(N, P) .< 0.5

@time yraw = X * realcoefs
#@time yrawsparse = X * sparsecoefs
using MLDataUtils
@time y = zerooneconvert.(yraw)
models = optimize(LinearIntegerModel, LassoLIMOptimizer(), X, y)

ms, aucs, preds = rankbyauc(models, X, y)

# Let's try SLOPE and see if similar result. Assumes RCall available and will install
# SLOPE R package which requires a functioning compiler build chain.
using RCall
function rlibrary(libname)
    @rput libname
    R"require(libname, character.only = TRUE)"
end
rlibrary("SLOPE")

struct SlopeLIMOptimizer <: ModelOptimizer
    withintercept::Bool
    maxnummodels::Int
    L::Int
    SlopeLIMOptimizer(; L::Int = 5, withintercept = true, maxnummodels = 30) =
        new(withintercept, maxnummodels, L)
end

function optimize(::Type{LinearIntegerModel}, o::SlopeLIMOptimizer,
    X::AbstractMatrix, y::AbstractVector)

    @rput X
    @rput y
    wintercept = o.withintercept
    @rput wintercept
    R"slopefit <- SLOPE(X, y, family='binomial', intercept=wintercept)"
    @rget slopefit

    nmodels = min(o.maxnummodels, size(lp.coefs, 2))
    wi = o.withintercept
    coefs = unique(map(i ->
        rescale_lasso_coefficients(lp.coefs[:, i], (wi ? lp.b0[i] : 0); L = o.L),
        2:nmodels))
    LinearIntegerModel[LinearIntegerModel(o.L, c...) for c in coefs]
end