text:
__precompile__()

module ProjectedGPE

using Reexport

@reexport using DifferentialEquations
@reexport using FastGaussQuadrature
@reexport using Parameters

@with_kw type Params @deftype Float64
#== paste this template into timeevolution and modify ==#
#fundamental constants/units
  ħ = 1.0545718e-34
  kB = 1.38064852e-23
  amu = 1.660339040e-27
  a₀ = 5.29e-11
#Rb87 mass and scattering length
  m = 86.909180527*amu
  as = 100*a₀
#trap frequencies
  ωx = 2π
  ωy = 0.
  ωz = 0.
#choice of time and length units
  t0 = 1.0/ωx
  x0 = sqrt(ħ*t0/m)
  E0 = ħ/t0
#interactions
  g  = (4*pi*ħ^2*as/m)*x0^3/E0 #dimensionless 3D
#damping parameters (dimensionless)
  γ  = 1e-4; @assert γ<=1.0
  ℳ  = 0.0;  @assert ℳ<=1.0
#chemical potential (dimensionless)
  μ  = 12.0
#time evolution parameters
  ti = 0.0
  tf = 100t0
  Nt::Int64 = 50
  t::Vector = collect(linspace(ti,tf,Nt))
  dt = 0.01π/μ #integrate step size [ - should have dt ≪ 2π/μ]
#== end template ==#
end

@with_kw type Cinfo
  γ::Float64=1e-4
  ℳ::Float64=0.0
  g::Float64=0.1
  μ::Float64=12.0
  basis::String="Hermite"
  Ω::Vector=[1]
  ecut::Number=10
  e0::Number=0.5
  Mmax::Number=10
  Mult::Number=10
  M::Vector=[10]
  en::Array=collect(0:9)+0.5
  P::BitArray=[true]
end

@with_kw type Tinfo
    Tx::Array{Float64,2}=zeros(2,2)
    Ty::Array{Float64,2}=zeros(2,2)
    Tz::Array{Float64,2}=zeros(2,2)
    W::Array=zeros(2,2,2)
end

include("eigmat.jl")
include("makespecops.jl")
include("nenergy.jl")
include("ladderops.jl")
include("makeangularmomentum.jl")
include("nfieldtrans.jl")
include("x2c.jl")
include("x2c!.jl")
include("xgrid2c.jl")
include("c2x!.jl")
include("c2x.jl")
include("growthrate.jl")
include("Skernel.jl")
include("makescatteringtrans.jl")
include("makescatteringnoisetrans.jl")
include("makecinfo.jl")
include("maketinfo.jl")
include("maketinfoplot.jl")
include("makealltrans.jl")
#include("pgpe.jl")
include("evolve.jl")

export eigmat, nfieldtrans, c2x, c2x!, x2c, x2c!, xgrid2c, growthrate,
Skernel, makescatteringtrans, makescatteringnoisetrans,
makecinfo, maketinfo, maketinfoplot, makealltrans, makespecops, ladderops,
makeangularmomentum, nenergy, gausshermite, evolve, Cinfo, Params, Tinfo

end # module


text:
using LightGraphs
using LightGraphs.SimpleGraphs
using LinearAlgebra
using SparseArrays
using DifferentialEquations
using ProgressMeter
#using DiffEqOperators
using Statistics

include("utils\\randutils.jl")
include("utils\\utils.jl")

struct Metapopulation{T<:AbstractSimpleGraph}
    g::T
    V::Vector{<:Integer} # node capacities
    β::Float64 # infection rate
    μ::Float64 # migration rate
end

"""
    metapopulation_gillespie(mp, s_0, i_0; kwargs)

Perform a single simulation of SI spreading on the metapopulation described by `mp`.

`s0` and `i0` contain respectively the initial numbers of susceptible and infected individuals per node in the metapopulation.

Keyword arguments
  * `nmax=1000`: maximum number of iterations of the algorithm
  * `tmax=100.0`: maximum time of the simulation
"""
function metapopulation_gillespie(mp::Metapopulation{SimpleGraph{T}}, s0, i0; nmax=1000, tmax=100.0) where T
    t = 0.0
    n = 1
    β = mp.β
    μ = mp.μ
    V = mp.V
    N = nv(mp.g)
    ts = Vector{typeof(t)}(undef, nmax)
    ts[1] = t
    ss = Array{typeof(s0[1]),2}(undef, nmax, length(s0))
    is = Array{typeof(i0[1]),2}(undef, nmax, length(i0))
    ss[1,:] = s0
    is[1,:] = i0
    s = copy(s0)
    i = copy(i0)
    d = degree(mp.g) .> 0 # which nodes have neighbors
    a = Vector{Float64}(undef, 3*N)
    a[1:N] = β*s.*i./(V.^2)
    a[N+1:2*N] = μ*s./V.*d # NB. aggregated probabilities of migration FROM each node
    a[2*N+1:3*N] = μ*i./V.*d # isolated nodes have migration probability of 0

    a0 = sum(a)
    while t < tmax && n < nmax && a0 > 0
        #println(a)
        τ = log(1/rand())/a0
        k = rand_categorical(a)
        if k <= N # infection in node k
            s[k] -= 1
            i[k] += 1
            a[k] = β*s[k]*i[k]/(V[k]^2)
            if d[k]
                a[N+k] = μ*s[k]/V[k]
                a[2*N+k] = μ*i[k]/V[k]
            end
        elseif k <= 2*N # migration of susceptible individual
            u = k-N
            v = rand(neighbors(mp.g,u)) # NB. will fail if there are no neighbors
            s[u] -= 1
            s[v] += 1
            a[u] = β*s[u]*i[u]/(V[u]^2)
            a[N+u] = μ*s[u]/V[u]
            a[v] = β*s[v]*i[v]/(V[v]^2)
            a[N+v] = μ*s[v]/V[v]
        else # migration of infected individual
            u = k-2*N
            v = rand(neighbors(mp.g,u)) # NB. will fail if there are no neighbors
            i[u] -= 1
            i[v] += 1
            a[u] = β*s[u]*i[u]/(V[u]^2)
            a[2*N+u] = μ*i[u]/V[u]
            a[v] = β*s[v]*i[v]/(V[v]^2)
            a[2*N+v] = μ*i[v]/V[v]
        end
        a0 = sum(a)
        t += τ
        ts[n+1] = t
        ss[n+1,:] = s
        is[n+1,:] = i
        n += 1
    end
    return ts[1:n], ss[1:n,:], is[1:n,:]
end

"""
    metapopulation_ode(mp, s0, i0; kwargs)

Integrate the Mean Field equations for SI spreading on the metapopulation described by `mp`.

`s0` and `i0` contain respectively the initial numbers of susceptible and infected individuals per node.

Keyword arguments:
  * `tmax=100.0`: the maximum time of the simulation
"""
function metapopulation_ode(mp::Metapopulation{<:AbstractSimpleGraph{<:Integer}}, s0, i0; tmax = 100.0, saveat=[])
    L = normalized_laplacian(mp.g)
    # define a special operator for solving the linear part of the equation exactly
    N = nv(mp.g)
    f! = function(dx,x,p,t)
        n = p
        tmp = mp.β*x[1:n].*x[n+1:2*n]./(mp.V).^2
        dx[1:n] = -tmp - mp.μ*L'*x[1:n]./mp.V
        dx[n+1:2*n] = tmp - mp.μ*L'*x[n+1:2*n]./mp.V
    end
    prob = ODEProblem(f!, vcat(float(s0),float(i0)), (0.0,tmax), N)
    sol = solve(prob, Tsit5(), saveat=saveat) # order 5 solver
    u = hcat(sol.u...)'
    s = u[:,1:N]
    i = u[:,N+1:2*N]
    return sol.t, s, i, sol
end

"""
    metapopulation_montecarlo(mp, s0, i0; kwargs)

Perform multiple simulations of SI spreading on the metapopulation described by `mp`, returning the average curve.

`s0` and `i0` contain respectively the initial numbers of susceptible and infected individuals in each node.

Keyword arguments:
  * `nmax=1000`: maximum number of iterations of each simulation
  * `tmax=100.0`: maximum time of the simulations
  * `nsims=100`: number of simulations to perform
  * `nbins`: the number of time steps at which the average is computed
"""
function metapopulation_montecarlo(mp::Metapopulation{SimpleGraph{T}}, s0, i0; nmax = 1000, tmax = 100.0, nsims = 100, nbins = 100) where T <: Integer
    N = nv(mp.g)
    ts = LinRange(0.0, tmax, nbins)
    Xs_sum = zeros(nbins, N)
    Ms_sum = zeros(nbins, N)
    Xi_sum = zeros(nbins, N)
    Mi_sum = zeros(nbins, N)
    Xn = Vector(undef, N)
    Δn = Vector(undef, N)
    p = Progress(nsims, dt=1.0)
    for n in 1:nsims
        t, s, i = metapopulation_gillespie(mp, s0, i0, tmax=tmax, nmax=nmax)
        l = 1
        for k in 1:nbins
            while l < length(t) && t[l+1] < ts[k]
                l += 1
            end
            Xn .= s[l,:]
            Δn .= Xn .- Xs_sum[k,:]/(n>1 ? n-1 : 1)
            Xs_sum[k,:] .+= Xn
            Ms_sum[k,:] .+= Δn.*(Xn.-Xs_sum[k]/n)
            Xn .= i[l,:]
            Δn .= Xn .- Xi_sum[k,:]/(n>1 ? n-1 : 1)
            Xi_sum[k,:] .+= Xn
            Mi_sum[k,:] .+= Δn.*(Xn.-Xi_sum[k]/n)
        end
        ProgressMeter.next!(p; showvalues= [(:n,n)])
    end
    return ts, Xs_sum/nsims, Xi_sum/nsims, sqrt.(abs.(Ms_sum/(nsims-1))), sqrt.(abs.(Mi_sum/(nsims-1)))
end

ext: jl

max_stars_repo_name: csimal/metaplex-networks

max_issues_repo_name: csimal/metaplex-networks

max_forks_repo_name: csimal/metaplex-networks


text:
function init_constraint_struct(com, ::GeqSetInternal, internals)
    GeqSetConstraint(internals, internals.indices[1], internals.indices[2:end], [])
end

"""
    init_constraint!(com::CS.CoM, constraint::GeqSetConstraint, fct::MOI.VectorOfVariables, set::GeqSetInternal)

Initialize the GeqSetConstraint by filling `sub_constraint_idxs` with all constraints
which variables are fully included in this constraint.
"""
function init_constraint!(
    com::CS.CoM,
    constraint::GeqSetConstraint,
    fct::MOI.VectorOfVariables,
    set::GeqSetInternal;
)
    for other_constraint in com.constraints
        if other_constraint.indices ⊆ constraint.indices[2:end]
            push!(constraint.sub_constraint_idxs, other_constraint.idx)
        end
    end
    return true
end

"""
    update_init_constraint!(com::CS.CoM, constraint::GeqSetConstraint, fct::MOI.VectorOfVariables,
        set::GeqSetInternal, constraints::Vector{<:Constraint})

Updates the `sub_constraint_idxs` when new constraints were added.
"""
function update_init_constraint!(
    com::CS.CoM,
    constraint::GeqSetConstraint,
    fct::MOI.VectorOfVariables,
    set::GeqSetInternal,
    constraints::Vector{<:Constraint},
)
    for other_constraint in constraints
        if other_constraint.indices ⊆ constraint.indices[2:end]
            push!(constraint.sub_constraint_idxs, other_constraint.idx)
        end
    end
    return true
end

"""
    prune_constraint!(
        com::CS.CoM,
        constraint::GeqSetConstraint,
        fct::MOI.VectorOfVariables,
        set::GeqSetInternal;
        logs = true
    )

Prune the constraint with:
- Finding the maximum of the minima as a lower bound for the `a` in `a >= X`
- Remove values bigger than the allowed maximum from `X`
- Set the lower bound of `a` even higher if there fully included `AllDifferentConstraints` allow it
"""
function prune_constraint!(
    com::CS.CoM,
    constraint::GeqSetConstraint,
    fct::MOI.VectorOfVariables,
    set::GeqSetInternal;
    logs = true,
)
    # find the maximum of the minima
    max_val = -typemax(Int)
    variables = com.search_space
    for vidx in constraint.greater_than
        max_val = max(variables[vidx].min, max_val)
    end
    feasible = remove_below!(com, variables[constraint.vidx], max_val)
    !feasible && return false

    # remove values bigger than the maximum of constraint.vidx
    # as this violates constraint.vidx >= constraint.greater_than
    max_val_variable = variables[constraint.vidx].max
    for vidx in constraint.greater_than
        feasible = remove_above!(com, variables[vidx], max_val_variable)
        !feasible && return false
    end

    # remove more values from constraint.vidx by checking the maximal
    # minimum value in all different constraints that are completely inside
    max_min = -typemax(Int)
    for cidx in constraint.sub_constraint_idxs
        sub_constraint = com.constraints[cidx]
        !(sub_constraint isa AllDifferentConstraint) && continue
        min_vals, max_vals = get_sorted_extrema(com, sub_constraint, 0, 0, 0)
        cmax_min = min_vals[1]
        for i in 2:length(sub_constraint.indices)
            if min_vals[i] <= cmax_min
                cmax_min += 1
            else
                cmax_min = min_vals[i]
            end
        end
        max_min = max(max_min, cmax_min)
    end
    feasible = remove_below!(com, variables[constraint.vidx], max_min)
    return feasible
end

"""
    still_feasible(
        com::CoM,
        constraint::GeqSetConstraint,
        fct::MOI.VectorOfVariables,
        set::GeqSetInternal,
        vidx::Int,
        value::Int,
    )

Check if the constraint is still feasible when setting vidx to value
"""
function still_feasible(
    com::CoM,
    constraint::GeqSetConstraint,
    fct::MOI.VectorOfVariables,
    set::GeqSetInternal,
    vidx::Int,
    value::Int,
)
    if vidx == constraint.vidx
        return true
    end
    variables = com.search_space
    return value <= variables[constraint.vidx].max
end

"""
    is_constraint_solved(
        constraint::GeqSetConstraint,
        fct::MOI.VectorOfVariables,
        set::GeqSetInternal,
        values::Vector{Int}
    )

Return true if `values` fulfills the constraint
"""
function is_constraint_solved(
    constraint::GeqSetConstraint,
    fct::MOI.VectorOfVariables,
    set::GeqSetInternal,
    values::Vector{Int},
)
    max_val = values[1]
    for i in 2:length(values)
        values[i] > max_val && return false
    end
    return true
end

"""
    is_constraint_violated(
        com::CoM,
        constraint::GeqSetConstraint,
        fct::MOI.VectorOfVariables,
        set::GeqSetInternal,
    )

Checks if the constraint is violated as it is currently set. This can happen inside an
inactive reified or indicator constraint.
"""
function is_constraint_violated(
    com::CoM,
    constraint::GeqSetConstraint,
    fct::MOI.VectorOfVariables,
    set::GeqSetInternal,
)
    for var in com.search_space[constraint.indices]
        if var.min > com.search_space[constraint.indices[1]].max
            return true
        end
    end
    return false
end


text:
export FullSumSampler

struct FullSumSampler <: FullSpaceSampler end

## Cache
mutable struct FullSumSamplerCache{T} <: SamplerCache{FullSumSampler}
    last_position::Int
    interval::T
end

_sampler_cache(s::FullSumSampler, v::FiniteBasisState, net, part) =
    FullSumSamplerCache(0, 1:spacedimension(v))

function init_sampler!(sampler::FullSumSampler, net, σ::FiniteBasisState, c::FullSumSamplerCache)
    c.last_position = 1
    # only initialize if it's really bigger.
    length(c.interval) >0 && set_index!(σ, c.interval[1])
    return c
end

done(s::FullSumSampler, σ, c) = c.last_position >= length(c.interval)

function samplenext!(σ, s::FullSumSampler, net, c)
    done(s, σ, c) && return false
    c.last_position += 1
    set_index!(σ,c.interval[c.last_position])
    return true
end

######## Multithreaded
_mt_recompute_sampler_params!(samplers, s::FullSumSampler) = nothing

function _divide_in_blocks(interval, rank, n_par)
    rank         = rank -1
    n_min, extra = divrem(length(interval), n_par)
    iter_length  = n_min + (rank < extra ? 1 : 0)
    iter_start   = n_min * rank + min(rank, extra) + 1
    iter_end     = iter_start + iter_length - 1

    return iter_start:iter_end
end

_sampler_cache(s::FullSumSampler, v::FiniteBasisState, net, ::ParallelThreaded, thread_i) =
    FullSumSamplerCache(0, _divide_in_blocks(1:spacedimension(v), thread_i, Threads.nthreads()))


text:
# %%
include("plotrecipe.jl")
using .MyPlot
# %%
MyPlot.cp
plot(rand(10))
# %%
MyPlot.set()
plot(rand(10))


text:
module TestInterface
using Plots
plotly()
using GaussianDistributions
using SARA
############################ Loading Data  #####################################

path = "/Users/sebastianament/Documents/SEA/XRD Analysis/SARA/Bi2O3_19F44_01/"
log10_τ, T, center_real, pos_wafer, coeffs = SARA.load_data(path)
num_coeffs, num_pos, num_stripes = size(coeffs)

###########################
position = pos_wafer[:, 1] .- center_real[1]
optical = [coeffs[c,:,1] for c in 1:size(coeffs, 1)]

# subsample
using Random: randperm
n = 8
rand_ind = randperm(num_pos)[1:n]
sub_pos = position[rand_ind]
sub_opt = [optical[i][rand_ind] for i in 1:num_coeffs]

# pre-process
domain = (minimum(position), maximum(position))
############################## run inner loop ##################################
c = 1/2
μ = 0.
σ = .25
β = 4.
parameters = [(c, μ, σ, β)]
rescale = SARA.get_rescale(parameters)
@time next, posteriors, max_uncertainty =
                        innerloop(sub_pos, sub_opt, domain;
                            grid = position, rescale_parameters = parameters)
# ################################ plotting ######################################
xs = position
for i in 1:1 #eachindex(posteriors)
    mu, std = evaluate(posteriors[i], xs)
    plot(position, optical[i], label = "complete data", legend = :bottomright)
    scatter!(sub_pos, sub_opt[i], label = "measured data")
    plot!(xs, mu, ribbon = 2std, label = "GP")
    scatter!(next, zero(next), label = "next")
    gui()
end

using GaussianDistributions: sample, mean, cov
using LinearAlgebra
# P = posteriors[1]
# M = marginal(P, xs)
# M = Gaussian(mean(M), Symmetric(cov(M)))
# tol = 1e-12
# factorize?
# K = cholesky(Symmetric(cov(M)), Val(true), check = false)
# P = cholesky(P)
# plot!(xs, sample(M, 3))
# gui()

######################## convert to T ###############################
ind = sortperm(T[:, 1])
temperature = T[ind, 1]/maximum(T[:,1])
position = position[ind]
domain = (minimum(temperature), maximum(temperature))
DT, UT = SARA.get_global_data(position, temperature, temperature, posteriors, domain)

########################## gradient map ##############################
n = 128
domain_scale = 1
temp = domain_scale * rand(n)
dwell = domain_scale * rand(n)
# create synthetic gradient map (two stripes)
∇(x) = exp(-((x[1]-.25)/.05)^2/2) + exp(-((x[1]-.75)/.05)^2/2)

# create synthetic gradient map (donut)
# m = 0.5
# s = .1
# ∇(x) = exp(-(norm(x.-fill(.5, length(x)))-m)^2/2s^2)

DT = ∇.(tuple.(temp, dwell))
UT = 1e-2exp.(randn(size(DT))) # non-uniform noise

grid = 0:1/50:1
grid = tuple.(grid, grid')
temperature_grid = [g[1] for g in grid]
dwelltime_grid = [g[2] for g in grid] # TODO: lazy grid?

domain = ((0., domain_scale * 1.), (0., domain_scale * 1.))
gradient_map = SARA.get_gradient_map(temp, dwell, DT, UT, domain)

M, S = SARA.evaluate_gradient_map(gradient_map, temperature_grid,
                                    dwelltime_grid, domain)

########################## outer loop ####################################
next_temp, next_dwell = SARA.outerloop(temp, dwell, gradient_map, temperature_grid[:], dwelltime_grid[:],
                                        domain, num_next = 1)
temp, dwell = 0:1/50:1, 0:1/50:1
next = (next_dwell, next_temp)
heatmap(dwell, temp, M)
scatter!([next])
gui()

end # TestSARA


text:
using PGEN
using Test, BGEN
const data = PGEN.datadir("bgen_example.16bits.pgen")
@testset "PGEN.jl" begin

@testset "Header" begin
    p = PGEN.Pgen(data)
    h = p.header
    @test h.bits_per_variant_type == 8
    @test h.bytes_per_record_length == 2
    @test h.n_blocks == 1
    @test h.n_samples == 0x01f4
    @test h.n_variants == 0xc7
    @test h.provisional_reference == 0x01
    @test h.provisional_reference_flags === nothing
    @test h.allele_counts === nothing
    @test h.storage_mode == 0x10
    @test h.variant_block_offsets[1] == 0x0269
    @test h.variant_lengths[1] == 0x04a2
    @test h.variant_lengths[end] == 0x04a0
    @test h.variant_types[1] == 0x60
    @test h.variant_types[end-1] == 0x41
    @test h.bytes_per_sample_id == 0x02
end

@testset "Difflist" begin
    # Write an example difflist given in PGEN spec.
    # Note that an offset of 1 is added to sample indexes to make it 1-based.
    io = open("dummy", "w")

    write(io, 0x4f)

    write(io, 0x88)
    write(io, 0x13)
    write(io, 0x00)

    write(io, 0x88)
    write(io, 0xf5)
    write(io, 0x04)

    write(io, 0x3f)

    for i in 1:20
        write(io, 0x00)
    end

    for i in 1:77
        write(io, 0x88)
        write(io, 0x27)
    end
    close(io)

    d = read("dummy")
    dl, offset = PGEN.parse_difflist(d, UInt(0), 3, true)
    @test dl.len == 79
    @test all(dl.genotypes .== 0)
    @test dl.has_genotypes
    @test unsafe_load(dl.last_component_sizes, 1) == 0x3f
    @test dl.sample_id_bases[1] == 5000
    @test dl.sample_id_bases[2] == 325000
    #@test length(dl.sample_id_increments[]) == 154
    idx = Vector{UInt32}(undef, 64)
    idx_incr = Vector{UInt32}(undef, 64)
    PGEN.parse_difflist_sampleids!(idx, idx_incr, dl, 1)
    @test all(idx .== [5000 * i for i in 1:64] .+ 1)
    PGEN.parse_difflist_sampleids!(idx, idx_incr, dl, 2)
    @test all(idx[1:15] .== [5000 * (64 + i) for i in 1:15] .+ 1) # for idx 65..79
    @test all(idx[16:end] .== 0)

    rm("dummy", force=true)
end

@testset "dosage" begin
    # NOTE: First alleles in the BGEN file are encoded as alternate allele in the transformation.
    # Some record types are not covered by this test, LD-compressions in particular.
    # They have been tested on a private UK Biobank data file.
    b = Bgen(PGEN.datadir("example.16bits.bgen"))
    p = PGEN.Pgen(data)
    g_pgen = Array{UInt8}(undef, p.header.n_samples)
    g_pgen_ld = similar(g_pgen)
    d_pgen = Array{Float64}(undef, p.header.n_samples)
    for (v_bgen, v_pgen) in zip(BGEN.iterator(b), PGEN.iterator(p)) #
        d_bgen = BGEN.ref_allele_dosage!(b, v_bgen)
        PGEN.alt_allele_dosage!(d_pgen, g_pgen, p, v_pgen)
        @test all(isapprox.(d_bgen, d_pgen; atol=5e-5, nans=true))
        PGEN.alt_allele_dosage!(d_pgen, g_pgen, p, v_pgen; genoldbuf=g_pgen_ld)
        @test all(isapprox.(d_bgen, d_pgen; atol=5e-5, nans=true))
        v_rt = v_pgen.record_type & 0x07
        if v_rt != 0x02 && v_rt != 0x03 # non-LD-compressed. See Format description.
            g_pgen_ld .= g_pgen
        end
    end
end
end


text:
## Formulations for Devices

# Dispatch formulations
abstract type DispatchForm <: DeviceFormulation end

struct LinearDispatchForm <: DispatchForm end
struct UnitCommittmentForm <: DispatchForm end

# Expansion Formulations
abstract type ExpansionForm{T<:DispatchForm} <: DeviceFormulation end
struct LinearExpansionForm{T<:DispatchForm} <: ExpansionForm{T} end
struct LinearExpansionDispatchForm{T<:DispatchForm} <: ExpansionForm{T} end
struct LinearExpansionInvestmentForm{T<:DispatchForm} <: ExpansionForm{T} end

struct BlockExpansionForm{T<:DispatchForm} <: ExpansionForm{T} end
struct BlockExpansionDispatchForm{T<:DispatchForm} <: ExpansionForm{T} end
struct BlockExpansionInvestmentForm{T<:DispatchForm} <: ExpansionForm{T} end

struct BinaryExpansionForm{T<:DispatchForm} <: ExpansionForm{T} end
struct BinaryExpansionDispatchForm{T<:DispatchForm} <: ExpansionForm{T} end
struct BinaryExpansionInvestmentForm{T<:DispatchForm} <: ExpansionForm{T} end

struct FormulationDescription
    name::Symbol
    form::DataType
end

FormulationDescription(tup) = FormulationDescription(tup...)

const formulations = FormulationDescription.([
    (:lindisp, LinearDispatchForm),
    (:ucdisp, UnitCommittmentForm),
    (:linexp_lindisp, LinearExpansionForm{LinearDispatchForm}),
    (:linexp_ucdisp, LinearExpansionForm{UnitCommittmentForm})
])

function resolve(::Type{DeviceFormulation}, name::Symbol)
    j = findfirst(fd -> fd.name == name, formulations)
    !isnothing(j) ? formulations[j].form : error("Formulation $name has not been registered in `formulations`")
end

function naming(::Type{T}) where T <: DeviceFormulation
    j = findfirst(fd -> fd.form == T, formulations)
    !isnothing(j) ? formulations[j].name : error("Formulation $T has not been registered in `formulations`")
end

formulation(::Device{DF}) where DF <: DeviceFormulation = DF

# TODO multi-year formulations

## ModelTypes

abstract type DispatchModel <: ModelType end
struct LinearDispatchModel <: DispatchModel end
struct UnitCommittmentModel <: DispatchModel end

abstract type ExpansionModel <: ModelType end

# If nothing else is said, let's say it's fine
demote_formulation(::Type{ExpansionModel}, ::Type{DF}) where DF <: ExpansionForm = DF
demote_formulation(::Type{ExpansionModel}, ::Type{DF}) where DF <: DispatchForm = DF
demote_formulation(::Type{<:DispatchModel},  ::Type{DF}) where DF <: DispatchForm = DF

# An expansion form should demote to its respective DispatchForm
demote_formulation(::Type{MT}, ::Type{DF}) where {MT <: DispatchModel, DDF, DF <: ExpansionForm{DDF}} = demote_formulation(MT, DDF)

# We demote UnitCommittment to LinearDispatch in a LinearDispatchModel
demote_formulation(::Type{LinearDispatchModel}, ::Type{UnitCommittmentForm}) = LinearDispatchForm

function demote_formulation(m::EnergyModel{MT}, d::Device{DF}) where {MT <: ModelType, DF <: DeviceFormulation}
    # TODO test whether conditionally changing the type would be faster!
    with_formulation(d, demote_formulation(MT, DF))
end

with_formulation(d::T, ::Type{DF}) where {T <: Device, DF <: DeviceFormulation} = T(d, DF)

max_stars_repo_name: coroa/EnergyModels

max_issues_repo_name: coroa/EnergyModels

max_forks_repo_name: coroa/EnergyModels


text:
using NCDatasets
using DataStructures

project_root_dir = @__DIR__
data_dir = joinpath(project_root_dir, "data")
domain_dir = joinpath(project_root_dir, "CESM_domains")
casename = "Sandbox"
Nz = 33
config = OrderedDict{Any, Any}(

    "DRIVER" => Dict(
        "casename"           => casename,
        "caseroot"           => joinpath(project_root_dir, casename, "caseroot"),
        "caserun"            => joinpath(project_root_dir, casename, "caserun"),
        "archive_root"       => joinpath(project_root_dir, casename, "archive"),
        "compute_QFLX_direct_method" => true,
    ),

    "MODEL_MISC" => Dict(
        "timetype"               => "DateTimeNoLeap",
        "init_file"              => joinpath(data_dir, "init_ocn.jld2"),
        "rpointer_file"          => "rpointer.iom",
        "daily_record"           => [],
        "monthly_record"         => ["{ESSENTIAL}", "QFLXT", "QFLXS"],
        "enable_archive"         => true,
    ),

    "MODEL_CORE" => Dict(

        #"domain_file"                  => joinpath(domain_dir, "domain.ocn.gx3v7.120323.nc"),
        "domain_file"                  => joinpath(domain_dir, "domain.ocn.gx1v6.090206.nc"),
        "topo_file"                    => joinpath(data_dir, "Nz_bot.nc"),
        "cdata_file"                   => joinpath(data_dir, "forcing.g37.nc"),

        "cdata_beg_time"               => "0001-01-01 00:00:00",
        "cdata_end_time"               => "0002-01-01 00:00:00",
        "cdata_align_time"             => "0001-01-01 00:00:00",

        "ϵ"                            => 0.7/86400,
        "γ"                            => 0.06,
        "z_w"                          => nothing,

        "substeps"                     => 8,
        "MLD_scheme"                   => "static",
        "Qflx"                         => "on",
        "Qflx_finding"                 => "off",
        "convective_adjustment"        => "on",
        "advection_scheme"             => "ekman_AGA2020",

        "weak_restoring"               => "off",
        "τwk_TEMP"                     => 86400.0 * 365 * 1000,
        "τwk_SALT"                     => 86400.0 * 365 * 1000,

        "τ_frz"                        => 3600.0,
        "Ekman_layers"                 => 5,
        "Returnflow_layers"            => 28,

        "transform_vector_field"       => true,
    ),

)

Dataset("data/coord.nc", "r") do ds

    println("Using $(Nz) vertical layers...")

    z_w_top = - nomissing(ds["z_w_top"][1:Nz], NaN) / 100.0
    z_w_bot = - nomissing(ds["z_w_bot"][1:Nz], NaN) / 100.0

    z_w      = zeros(Float64, length(z_w_top)+1)
    z_w[1:end-1]   = z_w_top
    z_w[end] = z_w_bot[end]

    global config["MODEL_CORE"]["z_w"]  = z_w
end

using TOML

open("data/config.toml", "w") do io
    TOML.print(io, config)
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
for i in 1:5
    println("Hiii number $i")
    sleep(1.0)
end
println("DONE")

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
#===============================================================================

    supsmu.jl

    wrapper for Fortran subroutine `supsmu` (Friedman 1984)

===============================================================================#

import Memento: debug, error

## location of executable
const libsupsmu = "$LOAD_FROM_DIR/_supsmu.so"

"Wrapper function for Fortran subroutine `supsmu` (Friedman 1984)."
function supsmu(
    X           ::AbstractVector,
    Y           ::AbstractVector,
    span        ::Real,
    wt          ::AbstractVector = ones(length(X)), ## weight
    periodic    ::Bool = false,
    alpha_bass  ::Real = 0, ## controls high frequency (small span) penality used with automatic span selection (bass tone control). (alpha.le.0.0 or alpha.gt.10.0 => no effect.)
    x_sorted    ::Bool = true
    # IntT_Fortran::DataType=Int32, FloatT_Fortran::DataType=Float32 ## Julia DataTypes of integer and float compatible with the Fortran subroutine `supsmu`
)
    debug(logger, "at supsmu()")
    #
    n = length(X)
    (ndims(X) > 1 || ndims(Y) > 1) && throw(ArgumentError, "X and Y must be 1-dimensional")
    (length(Y) != n)               && throw(ArgumentError, "X and Y must be equal in length")
    (n < 3)                        && warn(logger, "lengths of X and Y are less than 3, " *
                                            "no smoothing will be performed")
    (span < 0 || span > 1)         && throw(ArgumentError, "`span` must be between 0 and 1")
    if periodic && (minimum(x) < 0 || maximum(x) > 1) ## iper = 2
        throw(ArgumentError, "if X is assumed to be periodic, its values must be between 0 and 1")
    end
    iper = 1
    #
    ## check done
    FloatT_in = typeof(Y[1])
    if !x_sorted
        indice_sorted = sortperm(X)
        indice_back = sortperm(indice_sorted)
        X = X[indice_sorted]
        Y = Y[indice_sorted]
    end
    #
    ## convert types for compatibility with Fortran `supsmu`
    X, Y, span, wt, alpha_bass =
        (Array(X), Array(Y), span, wt, alpha_bass) |> cast(cast(Float32))
    n, iper = map(Int32, (n, iper))
    #
    ## variables to be modified
    smo = zeros(Float32, n)
    sc = zeros(Float32, n, 7)
    #
    ## call Fortran subroutine `supsmu`
    ccall((:supsmu_, libsupsmu), Void, (
            Ref{Int32}, Ref{Float32}, Ref{Float32}, ## n, X, Y
            Ref{Float32}, Ref{Int32}, ## wt, iper
            Ref{Float32}, Ref{Float32}, ## span, alpha_bass
            Ref{Float32}, Ref{Float32} ## smo, sc
        ), ## must be literal, cannot be variable(s) to be evaluated
        n, X, Y, wt, iper, span, alpha_bass, smo, sc
    ) ## smo and sc change values but not
    #
    ## convert back to the DataType of input
    smo = map(FloatT_in, smo)
    #
    return x_sorted ? smo : smo[indice_back]
end

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -19.57782063716216 <= q <= 17.16108472444359)

                     add_NL_constraint(m, :(log(1 + exp(-0.6503895050792288 + 0.5000448231799557*log(1 + exp(0.48086890691524875 + 0.8604652607112211*$(x[1]) + -0.9524910094017911*$(x[2]) + -0.8411501540974782*$(x[3]))) + 0.8223864208279696*log(1 + exp(0.4084009324685236 + 0.9106618206937234*$(x[1]) + -0.8345217492048942*$(x[2]) + 0.6009855404051798*$(x[3]))) + -0.39198040344556784*log(1 + exp(0.29882092390692216 + -0.9088484335369627*$(x[1]) + -0.948042223821898*$(x[2]) + 0.5765779510962643*$(x[3]))) + 0.4922116220556738*log(1 + exp(0.9223388106039003 + 0.14028925335959963*$(x[1]) + -0.026802684789215103*$(x[2]) + 0.9753923898070198*$(x[3]))))) + log(1 + exp(0.1539410450835783 + -0.07568393712605781*log(1 + exp(0.48086890691524875 + 0.8604652607112211*$(x[1]) + -0.9524910094017911*$(x[2]) + -0.8411501540974782*$(x[3]))) + -0.902375002810281*log(1 + exp(0.4084009324685236 + 0.9106618206937234*$(x[1]) + -0.8345217492048942*$(x[2]) + 0.6009855404051798*$(x[3]))) + 0.7665443117707627*log(1 + exp(0.29882092390692216 + -0.9088484335369627*$(x[1]) + -0.948042223821898*$(x[2]) + 0.5765779510962643*$(x[3]))) + -0.49047645888895497*log(1 + exp(0.9223388106039003 + 0.14028925335959963*$(x[1]) + -0.026802684789215103*$(x[2]) + 0.9753923898070198*$(x[3]))))) + log(1 + exp(0.6062617117178619 + -0.6503507748090844*log(1 + exp(0.48086890691524875 + 0.8604652607112211*$(x[1]) + -0.9524910094017911*$(x[2]) + -0.8411501540974782*$(x[3]))) + -0.9378359170676176*log(1 + exp(0.4084009324685236 + 0.9106618206937234*$(x[1]) + -0.8345217492048942*$(x[2]) + 0.6009855404051798*$(x[3]))) + -0.23896524724444168*log(1 + exp(0.29882092390692216 + -0.9088484335369627*$(x[1]) + -0.948042223821898*$(x[2]) + 0.5765779510962643*$(x[3]))) + -0.3347498056547815*log(1 + exp(0.9223388106039003 + 0.14028925335959963*$(x[1]) + -0.026802684789215103*$(x[2]) + 0.9753923898070198*$(x[3]))))) + log(1 + exp(-0.5413645906370745 + 0.4055343512605676*log(1 + exp(0.48086890691524875 + 0.8604652607112211*$(x[1]) + -0.9524910094017911*$(x[2]) + -0.8411501540974782*$(x[3]))) + -0.13388786735965885*log(1 + exp(0.4084009324685236 + 0.9106618206937234*$(x[1]) + -0.8345217492048942*$(x[2]) + 0.6009855404051798*$(x[3]))) + -0.9731450409901217*log(1 + exp(0.29882092390692216 + -0.9088484335369627*$(x[1]) + -0.948042223821898*$(x[2]) + 0.5765779510962643*$(x[3]))) + 0.17849747613123057*log(1 + exp(0.9223388106039003 + 0.14028925335959963*$(x[1]) + -0.026802684789215103*$(x[2]) + 0.9753923898070198*$(x[3]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
struct CykModel
    nrule::Int
    base2::Int
    labelid::Td{String, Int}
    idlabel::Vector{String}
    fa2lr::Vector{Vector{Tuple{Int, Int, Tf}}}
    ss2fa::Vector{Vector{Tuple{Int, Tf}}}
    lr2fa::Td{Int, Vector{Tuple{Int, Tf}}}
    l2far::Vector{Vector{Tuple{Int, Int, Tf}}}
end

function CykModel(ctb::ChTreebank)
    labelid = Td{String, Int}()
    falr = Td{Tuple{Int, Int, Int}, Tf}()
    fass = Td{Tuple{Int, Int}, Tf}()
    function visitor(cur::ChTree)::Nothing
        nchd = length(cur.adj)
        if isleaf(cur) || isposn(cur)
            ;
        elseif nchd == 1
            key = (cur.label, cur.adj[1].label)
            for pos in key if !haskey(labelid, pos) labelid[pos] = length(labelid) + 1 end end
            key = map(x -> labelid[x], key)
            fass[key] = Base.get(fass, key, 0) + 1
        else
            key = (cur.label, cur.adj[1].label, cur.adj[2].label)
            for pos in key if !haskey(labelid, pos) labelid[pos] = length(labelid) + 1 end end
            key = map(x -> labelid[x], key)
            falr[key] = Base.get(falr, key, 0) + 1
        end
        return nothing
    end
    for vt in ctb
        for tree in vt
            tree = cnf(tree)
            dfstraverse(tree, visitor)
        end
    end
    npos = length(labelid)
    base2 = 1; while (1 << base2) <= npos base2 += 1 end;
    # smoothing?
    fatot = Td{Int, Tf}()
    for (key, cnt) in falr
        fa = key[1]
        fatot[fa] = Base.get(fatot, fa, 0) + 1
    end
    for (key, cnt) in fass
        fa = key[1]
        fatot[fa] = Base.get(fatot, fa, 0) + 1
    end
    for (key, cnt) in falr
        fa = key[1]
        prob = -log(cnt / fatot[fa])
        falr[key] = prob
    end
    for (key, cnt) in fass
        fa = key[1]
        prob = -log(cnt / fatot[fa])
        fass[key] = prob
    end
    return CykModel(base2, labelid, falr, fass)
end
import Base.length
length(model::CykModel) = model.nrule

function CykModel(base2::Int, labelid::Td{String, Int}, falr::Td{Tuple{Int, Int, Int}, Tf}, fass::Td{Tuple{Int, Int}, Tf})
    idlabel = Vector{String}(undef, length(labelid))
    for (pos, id) in labelid idlabel[id] = pos end
    lr2fa = Td{Int, Vector{Tuple{Int, Tf}}}()
    for ((fa, l, r), prob) in falr
        lr = (l << base2) + r
        if !haskey(lr2fa, lr) lr2fa[lr] = Vector{Tuple{Int, Tf}}() end
        fas = lr2fa[lr]
        push!(fas, (fa, prob))
    end
    for ((fa, ss), prob) in fass
        if !haskey(lr2fa, ss) lr2fa[ss] = Vector{Tuple{Int, Tf}}() end
        fas = lr2fa[ss]
        push!(fas, (fa, prob))
    end
    nlabel = length(labelid)
    l2far = Vector{Vector{Tuple{Int, Int, Tf}}}(undef, nlabel)
    for i = 1:nlabel l2far[i] = Vector{Tuple{Int, Int, Tf}}() end
    for ((fa, l, r), prob) in falr
        push!(l2far[l], (fa, r, prob))
    end
    fa2lr = Vector{Vector{Tuple{Int, Int, Tf}}}(undef, nlabel)
    for i = 1:nlabel fa2lr[i] = Vector{Tuple{Int, Int, Tf}}() end
    for ((fa, l, r), prob) in falr
        push!(fa2lr[fa], (l, r, prob))
    end
    ss2fa = Vector{Vector{Tuple{Int, Tf}}}(undef, nlabel)
    for i = 1:nlabel ss2fa[i] = Vector{Tuple{Int, Tf}}() end
    for ((fa, ss), prob) in fass
        push!(ss2fa[ss], (fa, prob))
    end
    return CykModel(length(falr), base2, labelid, idlabel, fa2lr, ss2fa, lr2fa, l2far)
end

import Base: read, write, ==
==(x::CykModel, y::CykModel) = length(x) == length(y) && x.idlabel == y.idlabel && x.base2 == y.base2 && x.fa2lr == y.fa2lr && x.ss2fa == y.ss2fa

function read(io::IO, ::Type{CykModel})
    base2 = 0; labelid = Td{String, Int}(); lr2fa = Td{Int, Td{Int, Tf}}()
    fa2lr = Td{Tuple{Int, Int, Int}, Tf}()
    fa2ss = Td{Tuple{Int, Int}, Tf}()

    for line in eachline(io)
        cells = split(line); ncell = length(cells)
        if ncell == 1 base2 = parse(Int, cells[1]) end
        if ncell == 2 labelid[cells[1]] = parse(Int, cells[2]) end
        if ncell == 3
            fa = parse(Int, cells[1])
            ss = parse(Int, cells[2])
            prob = parse(Tf, cells[2])
            fa2ss[(fa, ss)] = prob
        end
        if ncell == 4
            fa = parse(Int, cells[1])
            l = parse(Int, cells[2])
            r = parse(Int, cells[3])
            prob = parse(Tf, cells[4])
            fa2lr[(fa, l, r)] = prob
        end
    end
    return CykModel(base2, labelid, fa2lr, fa2ss)
end

function write(io::IO, obj::CykModel)
    println(io, obj.base2)
    for (pos, id) in obj.labelid
        println(io, pos, " ", id)
    end
    for (falr, prob) in obj.fa2lr
        fa = falr[1]; l = falr[2]; r = falr[3];
        println(io, fa, " ", l, " ", r, " ", prob)
    end
    for (fass, prob) in obj.fa2ss
        fa = fass[1]; ss = fass[2];
        println(io, fa, " ", ss, " ", prob)
    end
end

import Base.get

const single = Td{String, Tf}()
get(model::CykModel, l::Int, r::Int) = get(model.lr2fa, (l << model.base2) + r, single)
get(model::CykModel, ss::Int) = get(model.lr2fa, ss, single)

function cyk(poswords::Vector{Tuple{String, String}}, model::CykModel)
    nrule = model.nrule
    labelid = model.labelid
    fa2lr = model.fa2lr
    ss2fa = model.ss2fa
    l2far = model.l2far
    nlabel = length(labelid)
    nword = length(poswords)
    Tf = Float64
    Td = Dict{Int, Tf}
    dp = Matrix{Td}(undef, nword, nword)
    Td2 = Dict{Int, Union{Tuple{Int, Int, Int}, Tuple{Int, Int}, Tuple{String}}}
    nx = Matrix{Td2}(undef, nword, nword)
    @inbounds begin
        for len = 1:nword
            @show len
            for i = 1:nword - len + 1
                res = Td()
                res2 = Td2()
                if len == 1
                    pos = labelid[poswords[i][1]]
                    res[pos] = 0.0
                    res2[pos] = (poswords[i][2],)
                else
                    for k = 1:len - 1
                        dpl = dp[i, k]; dpr = dp[i + k, len - k]
                        # @show length(dpl) * length(dpr), length(model.fa2lr)
                        if length(dpl) * length(dpr) << 1 >= nrule
                            # if false
                            for fa = 1:nlabel
                                for (l, r, prob) in fa2lr[fa]
                                    if haskey(dpl, l) && haskey(dpr, r)
                                        old = get(res, fa, Inf)
                                        lp = dpl[l]
                                        rp = dpr[r]
                                        new = lp + rp + prob
                                        if new < old
                                            res[fa] = new
                                            res2[fa] = (l, r, k)
                                        end
                                    end
                                end
                            end
                        else
                            for (l, pl) in dpl
                                for (fa, r, prob) in l2far[l]
                                    if haskey(dpr, r)
                                        old = get(res, fa, Inf)
                                        new = pl + dpr[r] + prob
                                        if new < old
                                            res[fa] = new
                                            res2[fa] = (l, r, k)
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
                while true
                    upd = false
                    for (ss, pss) in res
                        for (fa, prule) in ss2fa[ss]
                            old = get(res, fa, Inf)
                            new = prule + pss
                            if new < old res[fa] = new; res2[fa] = (ss, len); upd = true end
                        end
                    end
                    if !upd break end
                end
                dp[i, len] = res
                nx[i, len] = res2
            end
        end
    end
    start = findmin(dp[1, nword])
    function dfs(s, o, ilabel)::ChTree
        cur = ChTree(model.idlabel[ilabel], ChTree[])
        nxt = nx[s, o][ilabel]; nnxt = length(nxt)
        if nnxt == 1
            push!(cur.adj, ChTree(nxt[1], ChTree[]))
        elseif nnxt == 2
            push!(cur.adj, dfs(s, o, nxt[1]))
        else
            no = nxt[3]
            push!(cur.adj, dfs(s, no, nxt[1]))
            push!(cur.adj, dfs(s + no, o - no, nxt[2]))
        end
        return cur
    end
    return (start[1], dfs(1, nword, start[2]))
end


text:
using Test, Dates

using Mongoc: BSON, BSONObjectId
using BSONSerializer

include("TestModule.jl")

@testset "Resolve" begin
    @test BSONSerializer.typepathref(TestModule.ChildType) == ["TestModule", "ChildType"]
    @test BSONSerializer.resolve_typepath(BSONSerializer.typepathref(TestModule.ChildType)) == TestModule.ChildType
    @test_throws MethodError BSONSerializer.resolve_typepath("heyyou")
end

@BSONSerializable(TestModule.ChildType)
@BSONSerializable(TestModule.FatherType)
@BSONSerializable(TestModule.ManyDicts)
@BSONSerializable(TestModule.Periods)
@BSONSerializable(TestModule.SingletonStruct)
@BSONSerializable(TestModule.OptionNothing)
@BSONSerializable(TestModule.OptionMissing)
@BSONSerializable(TestModule.DateEncodedAsString)
@BSONSerializable(TestModule.Submodule.SubStruct)
@BSONSerializable(TestModule.StructFloat)
@BSONSerializable(TestModule.StructUInt64)
@BSONSerializable(TestModule.StructAbsTypes)
@BSONSerializable(TestModule.UnionPeriods)
@BSONSerializable(TestModule.Lift2)

function encode_roundtrip(v::T) where {T}
    BSONSerializer.decode(BSONSerializer.encode(v, T), T)
end

@testset "encode" begin
    let
        d = Dict{String, Int}("a" => 1, "b" => 2)
        new_d = encode_roundtrip(d)
        @test new_d == d
    end
end

@testset "serialize structs" begin
    @testset "ChildType" begin
        @test BSONSerializer.encode_type(TestModule.ChildType) == BSON

        instance = TestModule.ChildType(
            "Hello from ChildType",
            101,
            UInt8(255),
            Int8(120),
            UInt16(545),
            Int16(-12),
            DateTime(Dates.today()),
            Dates.today(),
            [Date(2019, 1, 1), Date(2019, 1, 2)],
            UInt32(12))

        bson = BSONSerializer.serialize(instance)
        #println(bson)
        @test bson["type"] == ["TestModule", "ChildType"]
        @test isa(bson["args"], Dict)
        args = bson["args"]
        @test args["c1"] == "Hello from ChildType"
        @test args["c2"] == 101
        @test args["c3"] == 255
        @test args["c4"] == 120
        @test args["c5"] == 545
        @test args["c6"] == -12
        @test args["c7"] == DateTime(Dates.today())
        @test args["c8"] == args["c7"]
        @test isa(args["c8"], DateTime)
        @test args["c9"] == [ DateTime(Date(2019, 1, 1)), DateTime(Date(2019, 1, 2)) ]
        @test args["c10"] == 12

        new_instance = BSONSerializer.deserialize(bson)
        @test new_instance == instance
    end

    @testset "FatherType" begin
        child_instance = TestModule.ChildType(
            "hey",
            101,
            UInt8(1),
            Int8(-1),
            UInt16(2),
            Int16(-2),
            DateTime(Dates.today()),
            Dates.today(),
            [Date(2019, 1, 1), Date(2019, 1, 2)],
            UInt32(12))

        oid = BSONObjectId()

        father_instance = TestModule.FatherType(
            "Hello from father",
            child_instance,
            [1, 2, 3],
            oid,
            :sym,
            '∈')

        bson = BSONSerializer.serialize(father_instance)
        #println(bson)
        @test bson["type"] == ["TestModule", "FatherType"]
        args = bson["args"]
        @test args["f1"] == "Hello from father"
        @test isa(args["f2"], Dict)
        @test args["f3"] == [ 1, 2, 3 ]
        child_dict = args["f2"]
        @test child_dict["type"] == ["TestModule", "ChildType"]
        @test isa(child_dict["args"], Dict)
        child_args = child_dict["args"]
        @test child_args["c1"] == "hey"
        @test child_args["c2"] == 101
        @test child_args["c3"] == 1
        @test child_args["c4"] == -1
        @test child_args["c5"] == 2
        @test child_args["c6"] == -2
        @test child_args["c7"] == DateTime(Dates.today())
        @test child_args["c8"] == child_args["c7"]
        @test child_args["c9"] == [ DateTime(Date(2019, 1, 1)), DateTime(Date(2019, 1, 2)) ]
        @test child_args["c10"] == 12
        @test args["f4"] == oid
        @test args["f5"] == "sym"
        @test args["f6"] == "∈"

        new_father_instance = BSONSerializer.deserialize(bson)
        @test new_father_instance == father_instance
    end

    @testset "IO" begin
        child_instance = TestModule.ChildType(
            "hey",
            101,
            UInt8(1),
            Int8(-1),
            UInt16(2),
            Int16(-2),
            DateTime(Dates.today()),
            Dates.today(),
            [Date(2019, 1, 1), Date(2019, 1, 2)],
            UInt32(12))

        oid = BSONObjectId()

        father_instance = TestModule.FatherType(
            "Hello from father",
            child_instance,
            [1, 2, 3],
            oid,
            :sym,
            '∈')

        io = IOBuffer()
        BSONSerializer.serialize(io, father_instance)
        seekstart(io)
        @test father_instance == BSONSerializer.deserialize(io)
    end

    @testset "Periods" begin
        instance = TestModule.Periods(Year(2000), Month(12), Day(20), Hour(23), Minute(59))
        new_instance = BSONSerializer.roundtrip(instance)
        @test new_instance == instance
    end

    @testset "ManyDicts" begin
        instance = TestModule.ManyDicts(
            Dict(:today => Dates.today(), :tomorrow => (Dates.today() + Dates.Day(1))),
            Dict(1 => Year(2000), 2 => Year(2001)),
            Dict(1 => TestModule.Periods(Year(2000), Month(10), Day(1), Hour(20), Minute(20))))
        new_instance = BSONSerializer.roundtrip(instance)
        @test new_instance == instance
    end

    @testset "SingletonStruct" begin
        instance = TestModule.SingletonStruct()
        new_instance = BSONSerializer.roundtrip(instance)
        @test new_instance == instance
    end

    @testset "Option" begin
        let
            instance = TestModule.OptionNothing(1)
            @test instance == BSONSerializer.roundtrip(instance)
        end

        let
            instance = TestModule.OptionNothing(nothing)
            @test instance == BSONSerializer.roundtrip(instance)
        end

        let
            instance = TestModule.OptionMissing("hey")
            @test instance == BSONSerializer.roundtrip(instance)
        end

        let
            instance = TestModule.OptionMissing(missing)
            @test instance == BSONSerializer.roundtrip(instance)
        end
    end

    @testset "DateEncodedAsString" begin
        instance = TestModule.DateEncodedAsString(Date(2019,12,25))
        bson_with_str = BSON("""
{
    "type" : ["TestModule", "DateEncodedAsString"],
    "args" : {
        "date" : "2019-12-25"
    }
}
""")
        new_instance = BSONSerializer.deserialize(bson_with_str)
        @test new_instance == instance
    end

    @testset "Submodule" begin
        instance = TestModule.Submodule.SubStruct(1)
        new_instance = BSONSerializer.roundtrip(instance)
        @test new_instance == instance
    end

    @testset "decode Int as Float" begin
        bson = BSON("""
{
    "type" : ["TestModule","StructFloat"],
    "args" : {
        "val" : 10
    }
}
""")
        new_instance = BSONSerializer.deserialize(bson)
        @test new_instance.val == 10
        @test isa(new_instance.val, Float64)
    end

    @testset "UInt64" begin
        let
            instance = TestModule.StructUInt64(typemax(UInt64))
            new_instance = BSONSerializer.roundtrip(instance)
            @test new_instance == instance
        end

        let
            instance = TestModule.StructUInt64(typemin(UInt64))
            new_instance = BSONSerializer.roundtrip(instance)
            @test new_instance == instance
        end
    end

    @testset "UnionPeriods" begin
        instance = TestModule.UnionPeriods(Dates.Day(1))
        bson = BSONSerializer.serialize(instance)
        @test bson["type"] == ["TestModule","UnionPeriods"]
        val = bson["args"]["val"]
        @test isa(val, Dict)
        @test haskey(val, "type")
        @test val["type"] == ["Dates", "Day"]
        @test haskey(val, "value")
        @test val["value"] == 1
        new_instance = BSONSerializer.deserialize(bson)
        @test new_instance == instance
    end
end

@BSONSerializable struct Tagged
    a::Int
    b::String
end

@testset "tagged struct" begin

    instance = Tagged(1, "hey")
    @test instance == BSONSerializer.roundtrip(instance)
end

@testset "abstract types" begin
    let
        instance = TestModule.StructAbsTypes(1)
        bson = BSONSerializer.serialize(instance)
        new_instance = BSONSerializer.deserialize(bson)
        @test new_instance == instance
    end

    let
        instance = TestModule.StructAbsTypes(Tagged(1, "hey"))
        new_instance = BSONSerializer.roundtrip(instance)
        @test new_instance == instance
    end

    let
        instance = TestModule.StructAbsTypes(TestModule.StructAbsTypes("omg"))
        new_instance = BSONSerializer.roundtrip(instance)
        @test new_instance == instance
    end

    let
        instance = TestModule.StructAbsTypes([TestModule.StructAbsTypes("omg"), 1])
        new_instance = BSONSerializer.roundtrip(instance)
        @test new_instance == instance
    end
end

@testset "encode functions" begin
    instance = TestModule.Lift2(+, 2, 3)
    new_instance = BSONSerializer.roundtrip(instance)
    @test TestModule.exec_op(new_instance) == 5
end

@testset "Usage" begin
    include("usage.jl")
end


text:
using Fst
using Base.Test

include("compose_tests.jl")
include("trim_tests.jl")

#=
println(Fst.states_with_no_in_arcs(a))
println(Fst.states_with_no_in_arcs(b))
println(Fst.states_with_no_in_arcs(c))

#println(topological_sort(a))
println(topological_sort(b))
println(topological_sort(c))
=#


text:
## Copyright (c) 2015, Júlio Hoffimann Mendes <juliohm@stanford.edu>
##
## Permission to use, copy, modify, and/or distribute this software for any
## purpose with or without fee is hereby granted, provided that the above
## copyright notice and this permission notice appear in all copies.
##
## THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
## WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
## MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
## ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
## WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
## ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
## OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

__precompile__()

module ImageQuilting

using Base: @nexprs, @nloops, @nref
using Images: imfilter_fft, padarray, dilate
using StatsBase: sample, weights
using LightGraphs

if VERSION > v"0.5-"
  using Combinatorics: nthperm!
end

try # optional dependencies
  using OpenCL
  using CLFFT
  global cl = OpenCL
  global clfft = CLFFT
catch
  global cl = nothing
  global clfft = nothing
end

include("utils.jl")
include("utils_gpu.jl")
include("datatypes.jl")
include("imfilter_gpu.jl")
include("relaxation.jl")
include("tau_model.jl")
include("dijkstra_cut.jl")
include("boykov_kolmogorov_cut.jl")
include("simplex_transform.jl")
include("mean_voxel_reuse.jl")
include("iqsim.jl")

export
  # functions
  iqsim,
  meanvoxreuse,

  # types
  SoftData,
  HardData

end


text:
"""
    NNCache{N}

The base abstract type holding a cache for the network type `N<:NeuralNetwork`.
"""
abstract type NNBatchedCache{N} <: NNCache{N}
    #valid::Bool
end


text:
#===============================================================================

    macros.jl

    container for data associated with a field in an Input struct

    Author: Tom Price
    Date:   July 2019

===============================================================================#

"Make a struct from predefined collection of Fields"
macro make_struct_from_SCHEMA(structname, super = Any, mutable = false)
    fields = SCHEMA |> mold(var_arg)
    esc(
        Expr(:type,
            mutable,
            Expr(:(<:), structname, super),
            Expr(:block, fields...)))
end ## macro

"Make a constructor method with default values from predefined collection of Fields"
macro make_constructor_from_SCHEMA(structname)
    no_constructor = SCHEMA |> mold(no_default) |> all
    if !no_constructor
        kw_args  = SCHEMA |> sift(!no_default) |> mold(kw_arg)
        var_args = SCHEMA |> sift( no_default) |> mold(var_arg)
        varnames = SCHEMA |> their(:name)
        esc(
            Expr(:(=),
                Expr(:call,
                    structname,
                    Expr(:parameters, kw_args...),
                    var_args...),
                Expr(:block,
                    Expr(:call, structname, varnames...))))
    end
end ## macro

"Get calibration data from request."
macro get_calibration_data_from_req(action)
    esc( :(
        begin
            req_key = curry(haskey)(req)
            req_key(CALIBRATION_INFO_KEY) ||
                return ArgumentError(
                    "no calibration data for " * string($action) * " analysis in request")
            calibration_data = try
                CalibrationData(req[CALIBRATION_INFO_KEY])
            catch err
                warn(logger, "calibration data exception: " * sprint(showerror, err))
                return ArgumentError(
                    "cannot parse calibration data for " * string($action) * " analysis")
            end ## try
        end ) )
end ## macro

"Parse raw data data from request."
macro parse_raw_data_from_req(action)
    esc( :(
        begin
            req_key(RAW_DATA_KEY) ||
                return ArgumentError(
                    "no raw data for " * string($action) * " analysis in request")
            parsed_raw_data = try
                parse_raw_data(Val{$action}, req[RAW_DATA_KEY])
            catch err
                warn(logger, "parse_raw_data exception: " * sprint(showerror, err))
                return ArgumentError(
                    "cannot parse raw data for " * string($action) * " analysis")
            end ## try
        end ) )
end ## macro

## usage:
# println(@macroexpand @make_struct MyStruct (a, Int_T) (b, Float64) )
# @makestruct MyStruct (a, Int_T) (b, Float64)
# macro make_struct(name, schema...)
#     fields = [  :( $(entry.args[1])::$(entry.args[2]) )
#                 for entry in schema ]
#     esc(:(
#         struct $name
#             $(fields...)
#         end ## end of struct
#         ) )
# end ## end of macro


text:
function SABRVolatility(F, X, T, ATMVol, Beta, VolVol, rho)

        alphaSABR(F, X, T, FindAlpha(F, X, T, ATMVol, Beta, VolVol, rho), Beta, VolVol, rho)

end

function CRoot(cubic, quadratic, linear, constant)

        roots = zeros(Float64, 3)
        a, b, C, r, Q, capA, capB, theta = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

        a = quadratic / cubic
        b = linear / cubic
        C = constant / cubic
        Q = (a^2 -3*b) / 9
        r = (2*a^3 - 9*a*b + 27*C) / 54

        if r^2 - Q^3 >= 0
                capA = -sign(r) * (abs(r) + sqrt(r^2 - Q^3))^(1/3)
                if capA == 0
                        capB = 0
                else
                        capB = Q/capA
                end

                croot = capA + capB - a/3

        else
                theta = acos(r/Q^1.5)

                roots[1] = -2 * sqrt(Q) * cos(theta / 3) - a / 3
                roots[2] = -2 * sqrt(Q) * cos(theta / 3 + 2.0943951023932) - a / 3
                roots[3] = -2 * sqrt(Q) * cos(theta / 3 - 2.0943951023932) - a / 3

                if roots[1] > 0
                        croot = roots[1]
                elseif roots[2] > 0
                        croot = roots[2]
                elseif roots[3] > 0
                        croot = roots[3]
                end

                if roots[2] > 0 && roots[2] < croot
                        croot = roots[2]
                end

                if roots[3] > 0 && roots[3] < croot
                        croot = roots[3]
                end

        end

return croot

end

function FindAlpha(F, X, T, ATMVol, Beta, VolVol, rho)

        cubic = (1-Beta)^2 * T / (24 * F^(2 - 2 * Beta))
        quadratic = 0.25 * rho * VolVol * Beta * T / F^(1 - Beta)
        linear = 1 + (2 - 3 * rho^2) / 24 * VolVol^2 * T
        constant = -ATMVol * F^(1 - Beta)

        CRoot(cubic, quadratic, linear, constant)

end

function alphaSABR(F, X, T, Alpha, Beta, VolVol, rho)

        dSABR = zeros(Float64, 3)
        #sabrz, y

        dSABR[1] = Alpha / ((F * X)^((1 - Beta) / 2) * (1 + (((1 - Beta)^2) / 24) * (log(F/X)^2) + ((1 - Beta)^4 / 1920) * (log(F / X)^4)))

        if abs(F - X)^2 > 10.0^(-8)
                sabrz = (VolVol / Alpha) * (F * X)^((1 - Beta) / 2) * log(F / X)
                y = (sqrt(1 - 2 * rho * sabrz + sabrz^2) + sabrz - rho) / (1 - rho)

                if abs(y - 1) < 10.0^(-8)
                        dSABR[2] = 1
                elseif y > 0
                        dSABR[2] = sabrz / log(y)
                else
                        dSABR[2] = 1
                end

        else
                dSABR[2] = 1

        end

        dSABR[3] = 1 + ((((1 - Beta)^2 / 24) * Alpha^2 / ((F * X)^(1 - Beta))) +
        0.25 * rho * Beta * VolVol * Alpha / ((F * X)^((1 - Beta) / 2)) +
        (2 - 3 * rho^2) * VolVol^2 / 24) * T

        dSABR[1] * dSABR[2] * dSABR[3]

end


text:
using Test
using BenchmarkTools
using LinearAlgebra
using BQCESubroutine
using YaoLocations

N = 20
M = 4
st = rand(Float64, 1<<N)
#st_mat = rand(Float64, 5, 1<<N)
U = rand(Float64, 1<<M, 1<<M)
locs = Locations(1:M)

subspace = BQCESubroutine.bsubspace(N, locs)
comspace = BQCESubroutine.bcomspace(N, locs)

nthreads = Threads.nthreads()
print("nthreads=$(nthreads), N=$N, M=$M")

if nthreads == 1
    BQCESubroutine.disable_threads()
    @btime st1 = BQCESubroutine.subspace_mul!(st, comspace, U, subspace)
else
    BQCESubroutine.enable_threads()
    @btime st1 = BQCESubroutine.threaded_subspace_mul_generic!(st, comspace, U, subspace)
end


text:
using MittagLefflerFunctions
using PyPlot

α = 1.85
β = 1.0
x = range(0, 10, length=201)
E1 = MLQuad3(α, β, 6)
E2 = MLQuad3(α, β, 8)

y = mlf.(α, β, -x, 10)
y1 = E1.(x)
y2 = E2.(x)

figure(1)
#plot(x, y1-y, x, y2-y)
semilogy(x, abs.(y1-y), x, abs.(y2-y))
grid(true)


text:
################################################################################
#                                ___  ____ _  _                                #
#                                |__] [__  |_/                                 #
#                                |    ___] | \_                                #
################################################################################

type PSK <: Modulation
    M::Int                    # Number of symbols in the constellation
    m::Int                    # Bits per symbol
    constellation::Vector     # Symbol constellation
    grayConstellation::Vector # Symbol constellation with gray code
    bitsMap::Vector           # bitsMap[i] maps to constellation[i]
end

function PSK( M::Integer )
    ispow2( M ) || error( "M must be a power of 2" )
    Δϕ                = 2π/M
    bitsMap           = [ encode( Gray, i ) for i in 0:M-1 ]
    # constellation     = Complex128[ exp(Δϕ*im*i) for i in 0:M-1 ]
    # grayConstellation = Complex128[ exp(Δϕ*im*decode( Gray, i)) for i in 0:M-1 ]
    constellation     = Array( Complex128, M )
    grayConstellation = Array( Complex128, M )
    for i in 0:M-1
       constellation[i+1]     = exp(Δϕ*im*i)
       grayConstellation[i+1] = exp(Δϕ*im*decode( Gray, i))
    end
    m                 = log2( M )
    PSK( M, m, constellation, grayConstellation, bitsMap )
end

################################################################################
#                                _  _ ____ ___                                 #
#                                |\/| |  | |  \                                #
#                                |  | |__| |__/                                #
################################################################################

function modulate( psk::PSK, bits::Integer )
    psk.grayConstellation[ bits+1 ]
end

function modulate( psk::PSK, data::AbstractVector )
    T = eltype(psk.grayConstellation)
    T[ modulate( psk, datum ) for datum in data ]
end

################################################################################
#                           ___  ____ _  _ ____ ___                            #
#                           |  \ |___ |\/| |  | |  \                           #
#                           |__/ |___ |  | |__| |__/                           #
################################################################################

function demodulate( psk::PSK, symbol::Complex )
    ϕ = angle( symbol )
    ϕ = ϕ < 0 ? ϕ += 2π : ϕ

    index = int( ϕ*psk.M / 2π ) + 1
    index = mod1( index, psk.M )
    psk.bitsMap[index]
end

function demodulate( psk::PSK, symbols::AbstractVector )
    Int[ demodulate( psk, symbol ) for symbol in symbols ]
end

################################################################################
#                      ____ _  _ ____ _  _ ___  _    ____                      #
#                      |___  \/  |__| |\/| |__] |    |___                      #
#                      |___ _/\_ |  | |  | |    |___ |___                      #
################################################################################

#=

using PyPlot

function PyPlot.plot( psk::PSK )
    fig = figure()
    scatter( real(psk.constellation), imag(psk.constellation) )
    plot( real(psk.grayConstellation), imag(psk.grayConstellation) )

    for i in 0:psk.M-1
        symbol = psk.constellation[ i + 1 ]
        x      = real( symbol )
        y      = imag( symbol )
        bits   = psk.bitsMap[i+1]

        annotate( bin(bits, psk.m ), (x, y+0.05), ha="center" )
    end

    setp( axes(), aspect = 1  )
    xlabel( "In Phase" )
    ylabel( "Quadrature" )

    return fig
end

psk   = PSK( 4 )

plot( psk )

modData   = 0:psk.M-1 #rand(0:8-1, 10)
symbols   = modulate( psk, modData)
demodData = demodulate( psk, symbols )

Test.@test_approx_eq modData demodData

=#


text:
using D91PVInversion
using IterativeSolvers
using LinearAlgebra
using SparseArrays
using Random
import PyPlot; const plt = PyPlot

# Define functional form of PV anomaly
function q′fun(x, y, z, A, L, H, p::Params)
    πc = p.π0 - 0.5
    Π = p.Π
    return A*exp(-(x^2 + y^2)/L^2)*exp(-(z-πc)^2/H^2)
end

# Set up inversion domain
params = Params(Float64)
domain = Domain(params, size = (3, 3, 3), x = (-1, 1)./2, y = (-1, 1)./2)
println(params)
println(domain)

# Set up background fields
ψ0, ϕ0, q0 = allocate_fields(domain)
set_background_ψ!(ψ0, domain, params)
set_background_ϕ!(ϕ0, domain, params)
fill_ψ_halos!(ψ0, domain, params)
fill_ϕ_halos!(ϕ0, domain, params)

# Allocate storage for fields and vectors
ψ′, ϕ′, q′ = allocate_fields(domain)
x, b = allocate_linearized_rhs(domain)

# Set initial values
@. ψ′ = 0
@. ϕ′ = 0

# Set PV field
A = 0.1
L = 0.1
H = 0.1
set_q′!(q′, domain, params, q′fun, A, L, H, params)

# Plot fields
plt.rc("font", size = 8)
fig, axes = plt.subplots(
    figsize = (9.5, 2), nrows = 1, ncols = 5,
    sharey = true, sharex = true, constrained_layout = true
)
for i = 1:5
    im = axes[i].imshow(ψ′[:,i-1,:]')
    plt.colorbar(im, ax = axes[i], location = "bottom", aspect = 5, shrink = 0.7)
    axes[i].set_title(string("ψ′, j = ", i))
end
axes[1].invert_yaxis()
plt.show()

# Construct linear operators
L = generate_linearized_L(ψ0, ϕ0, domain)
S = generate_sparse_linearized_L(ψ0, ϕ0, domain)

# Plot operators
ML = Matrix(L)
plt.figure()
plt.imshow(ML)
plt.title(string("L, rank = ", rank(ML), ", κ = ", cond(ML)))
plt.show()

MS = Matrix(S)
plt.figure()
plt.imshow(ML)
plt.title(string("S, rank = ", rank(MS), ", κ = ", cond(MS)))
plt.show()

plt.figure()
plt.imshow(ML - MS)
plt.show()

# Fill halos
fill_ψ′_halos!(ψ′, domain)
fill_ϕ′_halos!(ϕ′, domain)

# Re-plot fields
fig, axes = plt.subplots(
    figsize = (9.5, 2), nrows = 1, ncols = 5,
    sharey = true, sharex = true, constrained_layout = true
)
for i = 1:5
    im = axes[i].imshow(ψ′[:,i-1,:]')
    plt.colorbar(im, ax = axes[i], location = "bottom", aspect = 5, shrink = 0.7)
    axes[i].set_title(string("ψ′, j = ", i))
end
axes[1].invert_yaxis()
plt.show()

# Re-plot operators
M = Matrix(L)
Minv = inv(M)
plt.figure()
plt.imshow(M)
plt.show()
plt.figure()
plt.imshow(Minv)
plt.show()
plt.figure()
plt.imshow(Minv*M)
plt.show()

# Set RHS
set_linearized_b!(b, q′, domain)

# Plot q and RHS
fig, axes = plt.subplots(
    figsize = (9.5, 2), nrows = 1, ncols = 5,
    sharey = true, sharex = true, constrained_layout = true
)
for i = 1:5
    im = axes[i].imshow(q′[:,i-1,:]')
    plt.colorbar(im, ax = axes[i], location = "bottom", aspect = 5, shrink = 0.7)
    axes[i].set_title(string("q′, j = ", i))
end

# Solve
Random.seed!(1234)
idrs!(x, L, b; log = false, verbose = true)
fields_from_linearized_rhs!(ϕ′, ψ′, x, domain)

# Plot solution
fig, axes = plt.subplots(
    figsize = (9.5, 2), nrows = 1, ncols = 5,
    sharey = true, sharex = true, constrained_layout = true
)
for i = 1:5
    im = axes[i].imshow(ψ′[:,i-1,:]')
    plt.colorbar(im, ax = axes[i], location = "bottom", aspect = 5, shrink = 0.7)
    axes[i].set_title(string("ψ′, j = ", i))
end
axes[1].invert_yaxis()
plt.show()

# Solve directly with sparse LU decomposition
LU = lu(S)
println("S sparsity = ", nnz(S)/length(S))
println("L sparsity = ", nnz(LU.L)/length(LU.L))
println("U sparsity = ", nnz(LU.U)/length(LU.U))
ldiv!(x, LU, b)
fields_from_linearized_rhs!(ϕ′, ψ′, x, domain)

# Plot solution
fig, axes = plt.subplots(
    figsize = (9.5, 2), nrows = 1, ncols = 5,
    sharey = true, sharex = true, constrained_layout = true
)
for i = 1:5
    im = axes[i].imshow(ψ′[:,i-1,:]')
    plt.colorbar(im, ax = axes[i], location = "bottom", aspect = 5, shrink = 0.7)
    axes[i].set_title(string("ψ′, j = ", i))
end
axes[1].invert_yaxis()
plt.show()


text:
"""
abstract type AbstractReward end

Used to customize the reward function. If you want to use your own reward, you have to create a struct
(called `CustomReward` for example) and define the following methods:
- set_reward!(::StepPhase, env::RLEnv{CustomReward}, model::CPModel, current_status::Union{Nothing, Symbol})
- set_reward!(::DecisionPhase, env::RLEnv{CustomReward}, model::CPModel)
- set_reward!(::EndingPhase, env::RLEnv{CustomReward}, symbol::Symbol)

Then, when creating the `LearnedHeuristic`, you define it using `LearnedHeuristic{CustomReward}(agent::RL.Agent)`
and your functions will be called instead of the default ones.
"""

abstract type AbstractReward end

"""
    LearnedHeuristic{SR<:AbstractStateRepresentation, R<:AbstractReward, A<:ActionOutput}

The LearnedHeuristic is a value selection heuristic which is learned thanks to a training made by solving a certain amount
of problem instances from files or from an `AbstractModelGenerator`. From the RL point of view, this is an agent which is
learning how to choose an appropriate action (value to assign) when observing an `AbstractStateRepresentation`. The agent
learns thanks to rewards that are given regularly during the search. He wil try to maximize the total reward.

From the RL point of view, this LearnedHeuristic also plays part of the role normally played by the environment. Indeed, the
LearnedHeuristic stores the action space, the current state and reward. The other role that a classic RL environment has is describing
the consequences of an action: in SeaPearl, this is done by the CP part - branching, running fixPoint!, backtracking, etc...
"""
mutable struct LearnedHeuristic{SR<:AbstractStateRepresentation, R<:AbstractReward, A<:ActionOutput} <: ValueSelection
    agent::RL.Agent
    fitted_problem::Union{Nothing, Type{G}} where G
    fitted_strategy::Union{Nothing, Type{S}} where S <: SearchStrategy
    action_space::Union{Nothing, Array{Int64,1}}
    current_state::Union{Nothing, SR}
    reward::Union{Nothing, R}
    search_metrics::Union{Nothing, SearchMetrics}
    firstActionTaken::Bool
    LearnedHeuristic{SR, R, A}(agent::RL.Agent) where {SR, R, A}= new{SR, R, A}(agent, nothing, nothing, nothing, nothing, nothing, nothing, false)
end

"""
    (valueSelection::LearnedHeuristic)(::InitializingPhase, model::CPModel, x::Union{Nothing, AbstractIntVar}, current_status::Union{Nothing, Symbol})

Update the part of the LearnedHeuristic which act like an RL environment, by initializing it with informations caracteritic of the new CPModel. Creates an
initial observation with a false variable. A fixPoint! will be called before the LearnedHeuristic takes its first decision.
Finally, makes the agent call the process of the RL pre_episode_stage (basically making sure that the buffer is empty).
"""
function (valueSelection::LearnedHeuristic)(::Type{InitializingPhase}, model::CPModel)
    # create the environment
    valueSelection.firstActionTaken = false
    update_with_cpmodel!(valueSelection, model)
    # FIXME get rid of this => prone to bugs
    false_x = first(values(branchable_variables(model)))
    env = get_observation!(valueSelection, model, false_x)

    # Reset the agent, useful for things like recurrent networks
    Flux.reset!(valueSelection.agent)

    valueSelection.agent(RL.PRE_EPISODE_STAGE, env)
end

"""
    (valueSelection::LearnedHeuristic)(::StepPhase, model::CPModel, x::Union{Nothing, AbstractIntVar}, current_status::Union{Nothing, Symbol})

The step phase is the phase between every procedure call. It is the best place to get stats about the search tree, that's why
we put a set_metrics! function here. As those metrics could be used to design a reward, we also put a set_reward! function.
ATM, the metrics are updated after the reward assignment as the current_status given totally decribes the changes that are to be made.
Another possibility would be to have old and new metrics in memory.
"""
function (valueSelection::LearnedHeuristic)(PHASE::Type{StepPhase}, model::CPModel, current_status::Union{Nothing, Symbol})
    set_reward!(PHASE, valueSelection, model, current_status)
    # incremental metrics, set after reward is updated
    set_metrics!(PHASE, valueSelection.search_metrics, model, current_status)
    nothing
end

"""
    (valueSelection::LearnedHeuristic)(::DecisionPhase, model::CPModel, x::Union{Nothing, AbstractIntVar}, current_status::Union{Nothing, Symbol})

Observe, store useful informations in the buffer with agent(POST_ACT_STAGE, ...) and take a decision with the call to agent(PRE_ACT_STAGE).
The metrics that aren't updated in the StepPhase, which are more related to variables domains, are updated here. Once updated, they can
be used in the other set_reward! function as the reward of the last action will only be collected in the RL.POST_ACT_STAGE.
"""
function (valueSelection::LearnedHeuristic)(PHASE::Type{DecisionPhase}, model::CPModel, x::Union{Nothing, AbstractIntVar})
    # domain change metrics, set before reward is updated
    set_metrics!(PHASE, valueSelection.search_metrics, model, x)
    set_reward!(PHASE, valueSelection, model)
    model.statistics.lastVar = x
    env = get_observation!(valueSelection, model, x)

    #println("Decision  ", env.reward, " ", env.terminal, " ", env.legal_actions, " ", env.legal_actions_mask)
    if valueSelection.firstActionTaken
        valueSelection.agent(RL.POST_ACT_STAGE, env) # get terminal and reward
    else
        valueSelection.firstActionTaken = true
    end

    action = valueSelection.agent(env) # Choose action
    # TODO: swap to async computation once in deployment
    #@async valueSelection.agent(RL.PRE_ACT_STAGE, env, action) # Store state and action
    valueSelection.agent(RL.PRE_ACT_STAGE, env, action)

    return action_to_value(valueSelection, action, state(env), model)
end

"""
    (valueSelection::LearnedHeuristic)(::EndingPhase, model::CPModel, x::Union{Nothing, AbstractIntVar}, current_status::Union{Nothing, Symbol})

Set the final reward, do last observation.
"""
function (valueSelection::LearnedHeuristic)(PHASE::Type{EndingPhase}, model::CPModel, current_status::Union{Nothing, Symbol})
    # the RL EPISODE stops
    if valueSelection.firstActionTaken
        set_reward!(DecisionPhase,valueSelection, model)  #last decision reward for the previous action taken just before the ending Phase
    end
    set_reward!(PHASE, valueSelection, model, current_status)
    false_x = first(values(branchable_variables(model)))
    env = get_observation!(valueSelection, model, false_x, true)
    #println("EndingPhase  ", env.reward, " ", env.terminal, " ", env.legal_actions, " ", env.legal_actions_mask)

    valueSelection.agent(RL.POST_ACT_STAGE, env) # get terminal and reward

    valueSelection.agent(RL.POST_EPISODE_STAGE, env)  # let the agent see the last observation

    if CUDA.has_cuda()
        CUDA.reclaim()
    end
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: amp
using AWS.Compat
using AWS.UUIDs

"""
    create_alert_manager_definition(data, workspace_id)
    create_alert_manager_definition(data, workspace_id, params::Dict{String,<:Any})

Create an alert manager definition.

# Arguments
- `data`: The alert manager definition data.
- `workspace_id`: The ID of the workspace in which to create the alert manager definition.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: Optional, unique, case-sensitive, user-provided identifier to ensure the
  idempotency of the request.
"""
function create_alert_manager_definition(
    data, workspaceId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return amp(
        "POST",
        "/workspaces/$(workspaceId)/alertmanager/definition",
        Dict{String,Any}("data" => data, "clientToken" => string(uuid4()));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_alert_manager_definition(
    data,
    workspaceId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return amp(
        "POST",
        "/workspaces/$(workspaceId)/alertmanager/definition",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("data" => data, "clientToken" => string(uuid4())),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_rule_groups_namespace(data, name, workspace_id)
    create_rule_groups_namespace(data, name, workspace_id, params::Dict{String,<:Any})

Create a rule group namespace.

# Arguments
- `data`: The namespace data that define the rule groups.
- `name`: The rule groups namespace name.
- `workspace_id`: The ID of the workspace in which to create the rule group namespace.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: Optional, unique, case-sensitive, user-provided identifier to ensure the
  idempotency of the request.
- `"tags"`: Optional, user-provided tags for this rule groups namespace.
"""
function create_rule_groups_namespace(
    data, name, workspaceId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return amp(
        "POST",
        "/workspaces/$(workspaceId)/rulegroupsnamespaces",
        Dict{String,Any}("data" => data, "name" => name, "clientToken" => string(uuid4()));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_rule_groups_namespace(
    data,
    name,
    workspaceId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return amp(
        "POST",
        "/workspaces/$(workspaceId)/rulegroupsnamespaces",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "data" => data, "name" => name, "clientToken" => string(uuid4())
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_workspace()
    create_workspace(params::Dict{String,<:Any})

Creates a new AMP workspace.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"alias"`: An optional user-assigned alias for this workspace. This alias is for user
  reference and does not need to be unique.
- `"clientToken"`: Optional, unique, case-sensitive, user-provided identifier to ensure the
  idempotency of the request.
- `"tags"`: Optional, user-provided tags for this workspace.
"""
function create_workspace(; aws_config::AbstractAWSConfig=global_aws_config())
    return amp(
        "POST",
        "/workspaces",
        Dict{String,Any}("clientToken" => string(uuid4()));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_workspace(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return amp(
        "POST",
        "/workspaces",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("clientToken" => string(uuid4())), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_alert_manager_definition(workspace_id)
    delete_alert_manager_definition(workspace_id, params::Dict{String,<:Any})

Deletes an alert manager definition.

# Arguments
- `workspace_id`: The ID of the workspace in which to delete the alert manager definition.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: Optional, unique, case-sensitive, user-provided identifier to ensure the
  idempotency of the request.
"""
function delete_alert_manager_definition(
    workspaceId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return amp(
        "DELETE",
        "/workspaces/$(workspaceId)/alertmanager/definition",
        Dict{String,Any}("clientToken" => string(uuid4()));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_alert_manager_definition(
    workspaceId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return amp(
        "DELETE",
        "/workspaces/$(workspaceId)/alertmanager/definition",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("clientToken" => string(uuid4())), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_rule_groups_namespace(name, workspace_id)
    delete_rule_groups_namespace(name, workspace_id, params::Dict{String,<:Any})

Delete a rule groups namespace.

# Arguments
- `name`: The rule groups namespace name.
- `workspace_id`: The ID of the workspace to delete rule group definition.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: Optional, unique, case-sensitive, user-provided identifier to ensure the
  idempotency of the request.
"""
function delete_rule_groups_namespace(
    name, workspaceId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return amp(
        "DELETE",
        "/workspaces/$(workspaceId)/rulegroupsnamespaces/$(name)",
        Dict{String,Any}("clientToken" => string(uuid4()));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_rule_groups_namespace(
    name,
    workspaceId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return amp(
        "DELETE",
        "/workspaces/$(workspaceId)/rulegroupsnamespaces/$(name)",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("clientToken" => string(uuid4())), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_workspace(workspace_id)
    delete_workspace(workspace_id, params::Dict{String,<:Any})

Deletes an AMP workspace.

# Arguments
- `workspace_id`: The ID of the workspace to delete.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: Optional, unique, case-sensitive, user-provided identifier to ensure the
  idempotency of the request.
"""
function delete_workspace(workspaceId; aws_config::AbstractAWSConfig=global_aws_config())
    return amp(
        "DELETE",
        "/workspaces/$(workspaceId)",
        Dict{String,Any}("clientToken" => string(uuid4()));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_workspace(
    workspaceId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return amp(
        "DELETE",
        "/workspaces/$(workspaceId)",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("clientToken" => string(uuid4())), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_alert_manager_definition(workspace_id)
    describe_alert_manager_definition(workspace_id, params::Dict{String,<:Any})

Describes an alert manager definition.

# Arguments
- `workspace_id`: The ID of the workspace to describe.

"""
function describe_alert_manager_definition(
    workspaceId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return amp(
        "GET",
        "/workspaces/$(workspaceId)/alertmanager/definition";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_alert_manager_definition(
    workspaceId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return amp(
        "GET",
        "/workspaces/$(workspaceId)/alertmanager/definition",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_rule_groups_namespace(name, workspace_id)
    describe_rule_groups_namespace(name, workspace_id, params::Dict{String,<:Any})

Describe a rule groups namespace.

# Arguments
- `name`: The rule groups namespace.
- `workspace_id`: The ID of the workspace to describe.

"""
function describe_rule_groups_namespace(
    name, workspaceId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return amp(
        "GET",
        "/workspaces/$(workspaceId)/rulegroupsnamespaces/$(name)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_rule_groups_namespace(
    name,
    workspaceId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return amp(
        "GET",
        "/workspaces/$(workspaceId)/rulegroupsnamespaces/$(name)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_workspace(workspace_id)
    describe_workspace(workspace_id, params::Dict{String,<:Any})

Describes an existing AMP workspace.

# Arguments
- `workspace_id`: The ID of the workspace to describe.

"""
function describe_workspace(workspaceId; aws_config::AbstractAWSConfig=global_aws_config())
    return amp(
        "GET",
        "/workspaces/$(workspaceId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_workspace(
    workspaceId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return amp(
        "GET",
        "/workspaces/$(workspaceId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_rule_groups_namespaces(workspace_id)
    list_rule_groups_namespaces(workspace_id, params::Dict{String,<:Any})

Lists rule groups namespaces.

# Arguments
- `workspace_id`: The ID of the workspace.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: Maximum results to return in response (default=100, maximum=1000).
- `"name"`: Optional filter for rule groups namespace name. Only the rule groups namespace
  that begin with this value will be returned.
- `"nextToken"`: Pagination token to request the next page in a paginated list. This token
  is obtained from the output of the previous ListRuleGroupsNamespaces request.
"""
function list_rule_groups_namespaces(
    workspaceId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return amp(
        "GET",
        "/workspaces/$(workspaceId)/rulegroupsnamespaces";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_rule_groups_namespaces(
    workspaceId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return amp(
        "GET",
        "/workspaces/$(workspaceId)/rulegroupsnamespaces",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_tags_for_resource(resource_arn)
    list_tags_for_resource(resource_arn, params::Dict{String,<:Any})

Lists the tags you have assigned to the resource.

# Arguments
- `resource_arn`: The ARN of the resource.

"""
function list_tags_for_resource(
    resourceArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return amp(
        "GET",
        "/tags/$(resourceArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_tags_for_resource(
    resourceArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return amp(
        "GET",
        "/tags/$(resourceArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_workspaces()
    list_workspaces(params::Dict{String,<:Any})

Lists all AMP workspaces, including workspaces being created or deleted.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"alias"`: Optional filter for workspace alias. Only the workspaces with aliases that
  begin with this value will be returned.
- `"maxResults"`: Maximum results to return in response (default=100, maximum=1000).
- `"nextToken"`: Pagination token to request the next page in a paginated list. This token
  is obtained from the output of the previous ListWorkspaces request.
"""
function list_workspaces(; aws_config::AbstractAWSConfig=global_aws_config())
    return amp("GET", "/workspaces"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
end
function list_workspaces(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return amp(
        "GET", "/workspaces", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end

"""
    put_alert_manager_definition(data, workspace_id)
    put_alert_manager_definition(data, workspace_id, params::Dict{String,<:Any})

Update an alert manager definition.

# Arguments
- `data`: The alert manager definition data.
- `workspace_id`: The ID of the workspace in which to update the alert manager definition.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: Optional, unique, case-sensitive, user-provided identifier to ensure the
  idempotency of the request.
"""
function put_alert_manager_definition(
    data, workspaceId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return amp(
        "PUT",
        "/workspaces/$(workspaceId)/alertmanager/definition",
        Dict{String,Any}("data" => data, "clientToken" => string(uuid4()));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function put_alert_manager_definition(
    data,
    workspaceId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return amp(
        "PUT",
        "/workspaces/$(workspaceId)/alertmanager/definition",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("data" => data, "clientToken" => string(uuid4())),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    put_rule_groups_namespace(data, name, workspace_id)
    put_rule_groups_namespace(data, name, workspace_id, params::Dict{String,<:Any})

Update a rule groups namespace.

# Arguments
- `data`: The namespace data that define the rule groups.
- `name`: The rule groups namespace name.
- `workspace_id`: The ID of the workspace in which to update the rule group namespace.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: Optional, unique, case-sensitive, user-provided identifier to ensure the
  idempotency of the request.
"""
function put_rule_groups_namespace(
    data, name, workspaceId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return amp(
        "PUT",
        "/workspaces/$(workspaceId)/rulegroupsnamespaces/$(name)",
        Dict{String,Any}("data" => data, "clientToken" => string(uuid4()));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function put_rule_groups_namespace(
    data,
    name,
    workspaceId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return amp(
        "PUT",
        "/workspaces/$(workspaceId)/rulegroupsnamespaces/$(name)",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("data" => data, "clientToken" => string(uuid4())),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    tag_resource(resource_arn, tags)
    tag_resource(resource_arn, tags, params::Dict{String,<:Any})

Creates tags for the specified resource.

# Arguments
- `resource_arn`: The ARN of the resource.
- `tags`:

"""
function tag_resource(resourceArn, tags; aws_config::AbstractAWSConfig=global_aws_config())
    return amp(
        "POST",
        "/tags/$(resourceArn)",
        Dict{String,Any}("tags" => tags);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function tag_resource(
    resourceArn,
    tags,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return amp(
        "POST",
        "/tags/$(resourceArn)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("tags" => tags), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Deletes tags from the specified resource.

# Arguments
- `resource_arn`: The ARN of the resource.
- `tag_keys`: One or more tag keys

"""
function untag_resource(
    resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()
)
    return amp(
        "DELETE",
        "/tags/$(resourceArn)",
        Dict{String,Any}("tagKeys" => tagKeys);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function untag_resource(
    resourceArn,
    tagKeys,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return amp(
        "DELETE",
        "/tags/$(resourceArn)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("tagKeys" => tagKeys), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    update_workspace_alias(workspace_id)
    update_workspace_alias(workspace_id, params::Dict{String,<:Any})

Updates an AMP workspace alias.

# Arguments
- `workspace_id`: The ID of the workspace being updated.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"alias"`: The new alias of the workspace.
- `"clientToken"`: Optional, unique, case-sensitive, user-provided identifier to ensure the
  idempotency of the request.
"""
function update_workspace_alias(
    workspaceId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return amp(
        "POST",
        "/workspaces/$(workspaceId)/alias",
        Dict{String,Any}("clientToken" => string(uuid4()));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_workspace_alias(
    workspaceId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return amp(
        "POST",
        "/workspaces/$(workspaceId)/alias",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("clientToken" => string(uuid4())), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
    @eval Base.Experimental.@optlevel 1
end

using Core: SimpleVector, CodeInfo, NewvarNode, GotoNode
using Base.Meta: isexpr

const SSAValues = Union{Core.Compiler.SSAValue, JuliaInterpreter.SSAValue}

const trackedheads = (:method,)
const structdecls = (:_structtype, :_abstracttype, :_primitivetype)

export signature, rename_framemethods!, methoddef!, methoddefs!, bodymethod
export CodeEdges, lines_required, lines_required!, selective_eval!, selective_eval_fromstart!

include("utils.jl")
include("signatures.jl")
include("codeedges.jl")

# precompilation

if ccall(:jl_generating_output, Cint, ()) == 1
    kwdefine = NamedTuple{(:define,),Tuple{Bool}}
    for ct in (Vector{Any}, Set{Any})
        f = methoddef!
        m = which(f, Tuple{Function, ct, Frame, Expr, Int})
        @assert precompile(Tuple{typeof(f), Function, ct, Frame, Expr, Int})
        mbody = bodymethod(m)
        # @assert precompile(Tuple{mbody.sig.parameters[1], Bool, typeof(f), Function, ct, Frame, Expr, Int})
        @assert precompile(Tuple{Core.kwftype(typeof(f)), kwdefine, typeof(f), Function, ct, Frame, Expr, Int})
        f = methoddefs!
        @assert precompile(Tuple{typeof(f), Any, ct, Frame})
        @assert precompile(Tuple{Core.kwftype(typeof(f)), kwdefine, typeof(f), Function, ct, Frame})
    end
    @assert precompile(Tuple{typeof(rename_framemethods!), Any, Frame, Dict{Symbol,MethodInfo},
                             Vector{SelfCall}, Dict{Symbol,Union{Nothing, Bool, Symbol}}})
    @assert precompile(Tuple{typeof(rename_framemethods!), Any, Frame, Dict{Symbol,MethodInfo},
                             Vector{NamedTuple{(:linetop, :linebody, :callee, :caller),Tuple{Int64,Int64,Symbol,Union{Bool, Symbol}}}},
                             Dict{Symbol,Union{Bool, Symbol}}})
    @assert precompile(Tuple{typeof(identify_framemethod_calls), Frame})
    @assert precompile(Tuple{typeof(callchain), Vector{SelfCall}})
    @assert precompile(Tuple{typeof(callchain), Vector{NamedTuple{(:linetop, :linebody, :callee, :caller),Tuple{Int64,Int64,Symbol,Union{Bool, Symbol}}}}})

    @assert precompile(CodeLinks, (Int, Int))
    @assert precompile(CodeEdges, (Int,))
    @assert precompile(CodeEdges, (CodeInfo,))
    @assert precompile(add_links!, (Pair{Union{SSAValue,SlotNumber,NamedVar},Links}, Any, CodeLinks))
    @assert precompile(lines_required!, (Vector{Bool}, Set{NamedVar}, CodeInfo, CodeEdges))

    precompile(Tuple{typeof(setindex!),Dict{Union{GlobalRef, Symbol},Links},Links,Symbol})
    precompile(Tuple{typeof(setindex!),Dict{Union{GlobalRef, Symbol},Variable},Variable,Symbol})
end


text:
using ClassicalOrthogonalPolynomials, Plots, Test
import ArrayLayouts: diagonal

###
# We can solve ODEs like the Airy equation
#
# u(-1) = airyai(-1)
# u(1) = airyai(1)
# u'' = x * u
#
# using the ultraspherical spectral method.

T = Chebyshev()
C = Ultraspherical(2)
x = axes(T,1)
D = Derivative(x)

c = [T[[begin,end],:]; C \ ((D^2 - diagonal(x))*T)] \ [airyai(-1); airyai(1); zeros(∞)]
u = T*c

@test u[0.0] ≈ airyai(0.0)
plot(u)

##
# Lee & Greengard
# ε*u'' - x*u' + u = 0, u(-1) = 1, u(1) = 2
##

T = ChebyshevT()
C = Ultraspherical(2)
x = axes(T,1)
D = Derivative(x)

ε = 1/100
A = [T[[begin,end],:]; C \ ((ε*D^2 - x .* D + I) * T)]
c = A \ [1; 2; zeros(∞)]
u = T*c
plot(u)

C \ (ε*D^2 - x .* D + I)


text:
using MD5: md5

d5_input = "iwrupvqb"

function find_with_prefix(prefix::String)::Int
    curr = 1
    while true
        h = bytes2hex(md5("$d5_input$curr"))
        if startswith(h, prefix)
            break
        else
            curr += 1
        end
    end
    curr
end

part_1 = find_with_prefix("00000")
println("Part 1: $part_1")
part_2 = find_with_prefix("000000")
println("Part 2: $part_2")


text:
include("match_low_volfrac.jl")
include("match_low_freq.jl")
include("match_average_wave.jl")