text:
# linear operators on fields
# ============================================================

"""
Interface requirements for AbstractLinearOp
```
Base.:inv(O::AbstractLinearOp)
LinearAlgebra.adjoint(O::AbstractLinearOp)
Base.:*(O::AbstractLinearOp, f::Field)
Base.:\\(O::AbstractLinearOp, f::Field)
```
Optionally
```
_lmult(O::AbstractLinearOp, f::Field) = O * f
```
can be specialized to avoid a final conversion back to
the argument field type as `*` would typically do.
Useful for a sequence of operations that avoid un-necessary conversions.
"""
abstract type AbstractLinearOp end

_lmult(O::AbstractLinearOp, f) = O * f

# Chains of AbstractLinearOp
# ------------------------------------------

# `*` concatinates chains
function Base.:*(O1::AbstractLinearOp, O2::AbstractLinearOp)
	tuple(O1, O2)
end

function Base.:*(O1::NTuple{N,AbstractLinearOp}, O2::AbstractLinearOp) where N
	Base.front(O1) * (Base.last(O1) * O2)
end

function Base.:*(O1::AbstractLinearOp, O2::NTuple{N,AbstractLinearOp}) where N
	(O1 * Base.first(O2)) * Base.tail(O2)
end

function Base.:*(O1::NTuple{N,AbstractLinearOp}, O2::NTuple{M,AbstractLinearOp}) where {N,M}
	tuple(O1..., O2...)
end

# `inv` broadcasts and reverses order
function Base.:inv(O1::NTuple{N,AbstractLinearOp}) where N
	tuple((inv(op) for op in reverse(O1))...)
end

# `adjoint` broadcasts and reverses order
function Base.:adjoint(O1::NTuple{N,AbstractLinearOp}) where N
	tuple((adjoint(op) for op in reverse(O1))...)
end

# `\`
function Base.:\(O1::AbstractLinearOp, O2::AbstractLinearOp)
	inv(O1) * O2
end

function Base.:\(O1::NTuple{N,AbstractLinearOp}, O2::AbstractLinearOp) where N
	inv(O1) * O2
end

function Base.:\(O1::AbstractLinearOp, O2::NTuple{N,AbstractLinearOp}) where N
	inv(O1) * O2
end

function Base.:\(O1::NTuple{N,AbstractLinearOp}, O2::NTuple{M,AbstractLinearOp}) where {N,M}
	inv(O1) * O2
end

# `/`
function Base.:/(O1::AbstractLinearOp, O2::AbstractLinearOp)
	O1 * inv(O2)
end

function Base.:/(O1::NTuple{N,AbstractLinearOp}, O2::AbstractLinearOp) where N
	O1 * inv(O2)
end

function Base.:/(O1::AbstractLinearOp, O2::NTuple{N,AbstractLinearOp}) where N
	O1 * inv(O2)
end

function Base.:/(O1::NTuple{N,AbstractLinearOp}, O2::NTuple{M,AbstractLinearOp}) where {N,M}
	O1 * inv(O2)
end

# activate the lazy tuple when operating
function Base.:*(O1::NTuple{N,AbstractLinearOp}, f::Y) where {N,Y<:Field}
	Y(foldr(_lmult, (O1..., f)))
end

function Base.:\(O1::NTuple{N,AbstractLinearOp}, f::Y) where {N,Y<:Field}
	(inv(O1) * f)::Y
end


text:
# # Inverse Lighting Tutorial
#
# In this tutorial we shall explore the inverse lighting problem.
# Here, we shall try to reconstruct a target image by optimizing
# the parameters of the light source (using gradients).

using RayTracer, Images, Zygote, Flux, Statistics

# ## Configuring the Scene
#
# Reduce the screen_size if the optimization is taking a bit long
screen_size = (w = 300, h = 300)

# Now we shall load the scene using [`load_obj`](@ref) function. For
# this we need the [`obj`](https://en.wikipedia.org/wiki/Wavefront_.obj_file)
# and [`mtl`](https://en.wikipedia.org/wiki/Wavefront_.obj_file#Material_template_library)
# files. They can be downloaded using the following commands:
#
# ```
# wget https://raw.githubusercontent.com/tejank10/Duckietown.jl/master/src/meshes/tree.obj
# wget https://raw.githubusercontent.com/tejank10/Duckietown.jl/master/src/meshes/tree.mtl
# ```
scene = load_obj("./tree.obj")

# Let us set up the [`Camera`](@ref). For a more detailed understanding of
# the rendering process look into [Introduction to rendering using RayTracer.jl](@ref).
cam = Camera(
    lookfrom = Vec3(0.0f0, 6.0f0, -10.0f0),
    lookat   = Vec3(0.0f0, 2.0f0,  0.0f0),
    vup      = Vec3(0.0f0, 1.0f0,  0.0f0),
    vfov     = 45.0f0,
    focus    = 0.5f0,
    width    = screen_size.w,
    height   = screen_size.h
)

origin, direction = get_primary_rays(cam)

# We should define a few convenience functions. Since we are going to calculate
# the gradients only wrt to `light` we have it as an argument to the function. Having
# `scene` as an additional parameters simply allows us to test our method for other
# meshes without having to run `Zygote.refresh()` repeatedly.
function render(light, scene)
    packed_image = raytrace(origin, direction, scene, light, origin, 2)
    array_image = reshape(hcat(packed_image.x, packed_image.y, packed_image.z),
                          (screen_size.w, screen_size.h, 3, 1))
    return array_image
end

showimg(img) = colorview(RGB, permutedims(img[:,:,:,1], (3,2,1)))

# ## [Ground Truth Image](@id inv_light)
#
# For this tutorial we shall use the [`PointLight`](@ref) source.
# We define the ground truth lighting source and the rendered image. We
# will later assume that we have no information about this lighting
# condition and try to reconstruct the image.
light_gt = PointLight(
    color     = Vec3(1.0f0, 1.0f0, 1.0f0),
    intensity = 20000.0f0,
    position  = Vec3(1.0f0, 10.0f0, -50.0f0)
)

target_img = render(light_gt, scene)

# The presence of [`zeroonenorm`](@ref) is very important here. It rescales the
# values in the image to 0 to 1. If we don't perform this step `Images` will
# clamp the values while generating the image in RGB format.
showimg(zeroonenorm(target_img))

# ```@raw html
# <p align="center">
#     <img width=300 height=300 src="../../assets/inv_light_original.png">
# </p>
# ```

# ## Initial Guess of Lighting Parameters
#
# We shall make some arbitrary guess of the lighting parameters (intensity and
# position) and try to get back the image in [Ground Truth Image](@ref inv_light)
light_guess = PointLight(
    color     = Vec3(1.0f0, 1.0f0, 1.0f0),
    intensity = 1.0f0,
    position  = Vec3(-1.0f0, -10.0f0, -50.0f0)
)

showimg(zeroonenorm(render(light_guess, scene)))

# ```@raw html
# <p align="center">
#     <img width=300 height=300 src="../../assets/inv_light_initial.png">
# </p>
# ```

# We shall store the images in `results_inv_lighting` directory
mkpath("results_inv_lighting")

save("./results_inv_lighting/inv_light_original.png",
     showimg(zeroonenorm(render(light_gt, scene))))
save("./results_inv_lighting/inv_light_initial.png",
     showimg(zeroonenorm(render(light_guess, scene))))

# ## Optimization Loop
#
# We will use the ADAM optimizer from Flux. (Try experimenting with other
# optimizers as well). We can also use frameworks like Optim.jl for optimization.
# We will show how to do it in a future tutorial

for i in 1:401
    loss, back_fn = Zygote.forward(light_guess) do L
        sum((render(L, scene) .- target_img) .^ 2)
    end
    @show loss
    gs = back_fn(1.0f0)
    update!(opt, light_guess.intensity, gs[1].intensity)
    update!(opt, light_guess.position, gs[1].position)
    if i % 5 == 1
        save("./results_inv_lighting/iteration_$i.png",
             showimg(zeroonenorm(render(light_guess, scene))))
    end
end

# If we generate a `gif` for the optimization process it will look similar to this
# ```@raw html
# <p align="center">
#      <img width=300 height=300 src="../../assets/inv_lighting.gif">
# </p>
# ```

ext: jl


text:
#=
Created on Saturday 26 December 2020
Last update: Tuesday 28 September 2021

@author: Bram De Jaegher
bram.de.jaegher@gmail.com
=#

"""
The @safe macro is a hidden try-catch statement to avoid the Markdown admonitions to crash when
the user introduces an error in the exercise functions.

Usage:
-----
@safe test_func(input) == correct_output
@safe[test_func1(input) == correct_output1, test_func2(input) == correct_output2]

"""
macro safe(ex)
	safe_ex = quote
		try $(esc(ex))
		catch e
			false
		end
	end
	return safe_ex
end

"""
The @terminal prints the output in the Pluto notebooks instead of the REPL.
This is a convenience wrapper for `PlutoUI.with_terminal()`.

"""
macro terminal(ex)
	return quote
			PlutoUI.with_terminal() do
					$(esc(ex))
			end
	end
end

max_stars_repo_name: Beramos/PlutoTeaching

max_issues_repo_name: Beramos/PlutoTeaching

max_forks_repo_name: Beramos/PlutoTeaching


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
struct Param{T} _inner::Dict{Symbol, T} end
Param() = Param(Dict{Symbol, Any}())
Base.getproperty(param::Param, name::Symbol) = getfield(param, :_inner)[name]
Base.setproperty!(param::Param, name::Symbol, val) = setindex!(getfield(param, :_inner), val, name)
Base.propertynames(param::Param) = (keys(getfield(param, :_inner))...,)

# %%
p = Param()

# %%
propertynames(p)

# %%
p.a = 123
p.b = 4.56
p.c = "hoge"
p

# %%
propertynames(p)

# %%
p.

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
## Let's solve a random problem
using SigmoidalProgramming
using Random
import DataStructures: dequeue!

# problem data
Random.seed!(4)
nvar = 200
nconstr = 40

l = -rand(nvar)
u = rand(nvar)
A = rand(nconstr, nvar)
b = rand(nconstr)
z = zeros(nvar)
fs = fill(logistic, nvar)
dfs = fill(logistic_prime, nvar)
problem = LinearSP(fs, dfs, z, A, b)

# branch and bound
pq, bestnodes, lbs, ubs = solve_sp(l, u, problem)

node = dequeue!(pq)
# println("best node has node.ub = $(node.ub) and solution $(node.x)")
println("lbs: ",lbs)
println("ubs: ",ubs)


text:
# Basic scitype

ST.scitype(::Integer,        ::MLJ) = Count
ST.scitype(::AbstractFloat,  ::MLJ) = Continuous
ST.scitype(::AbstractString, ::MLJ) = Textual
ST.scitype(::TimeType,       ::MLJ) = ScientificTimeType
ST.scitype(::Time    ,       ::MLJ) = ScientificTime
ST.scitype(::Date,           ::MLJ) = ScientificDate
ST.scitype(::DateTime,       ::MLJ) = ScientificDateTime

ST.scitype(img::Arr{<:Gray,2}, ::MLJ) = GrayImage{size(img)...}
ST.scitype(img::Arr{<:AbstractRGB,2}, ::MLJ) =
ColorImage{size(img)...}

# CategoricalArray scitype

function ST.scitype(c::Cat, ::MLJ)
    nc = length(levels(c.pool))
    return ifelse(c.pool.ordered, OrderedFactor{nc}, Multiclass{nc})
end

function ST.scitype(A::CArr{T,N}, ::MLJ) where {T,N}
    nlevels = length(levels(A))
    S = ifelse(isordered(A), OrderedFactor{nlevels}, Multiclass{nlevels})
    T >: Missing && (S = Union{S,Missing})
    return AbstractArray{S,N}
end

# Table scitype

function ST.scitype(X, ::MLJ, ::Val{:table}; kw...)
    Xcol = Tables.columns(X)
    col_names = propertynames(Xcol)
    types = map(col_names) do name
        scitype(getproperty(Xcol, name); kw...)
    end
    return Table{Union{types...}}
end

# Scitype for fast array broadcasting

ST.Scitype(::Type{<:Integer},        ::MLJ) = Count
ST.Scitype(::Type{<:AbstractFloat},  ::MLJ) = Continuous
ST.Scitype(::Type{<:AbstractString}, ::MLJ) = Textual
ST.Scitype(::Type{<:TimeType},         ::MLJ) = ScientificTimeType
ST.Scitype(::Type{<:Date},             ::MLJ) = ScientificDate
ST.Scitype(::Type{<:Time},             ::MLJ) = ScientificTime
ST.Scitype(::Type{<:DateTime},         ::MLJ) = ScientificDateTime


text:
__precompile__()
module Mosek
using SparseArrays

if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("Mosek not properly installed. Please run Pkg.build(\"Mosek\")")
end

# The MathOptInterface wrapper is in MosekTools.jl but it can be created with
# `using MosekTools; Mosek.Optimizer()`.
function Optimizer(args...; kwargs...)
    error("To use Mosek with JuMP (or MathOptInterface), you need to use ",
          "the package `MosekTools` (via `using MosekTools`). You may need ",
          "to first install it via `import Pkg; Pkg.add(\"MosekTools\")`.")
end

export
    makeenv, maketask, maketask_ptr,
    MosekError

# A macro to make calling C API a little cleaner
macro msk_ccall(func, args...)
    f = Base.Meta.quot(Symbol("MSK_$(func)"))
    args = [esc(a) for a in args]
    quote
        ccall(($f,libmosek), $(args...))
    end
end

# -----
# Types
# -----
struct MosekError <: Exception
    rcode :: Int32
    msg   :: String
end

# Environment: typedef void * Env_t;
mutable struct Env
    env::Ptr{Nothing}
    streamcallbackfunc::Any
end

function msk_info_callback_wrapper end

# Task: typedef void * Task_t;
mutable struct Task
    env::Env
    task::Ptr{Nothing}
    borrowed::Bool
    # need to keep a reference to callback funcs for GC
    streamcallbackfunc:: Any
    userstreamcallbackfunc:: Any
    callbackfunc:: Any
    usercallbackfunc:: Any
    nlinfo:: Any

    function Task(env::Env)
        temp = Array{Ptr{Nothing}}(undef,1)
        res = @msk_ccall(maketask, Int32, (Ptr{Nothing}, Int32, Int32, Ptr{Nothing}), env.env, 0, 0, temp)

        if res != MSK_RES_OK.value
            throw(MosekError(res,""))
        end

        task = new(env,temp[1],false,nothing,nothing,nothing,nothing,nothing)
        task.callbackfunc = @cfunction(msk_info_callback_wrapper, Cint, (Ptr{Nothing}, Ptr{Nothing}, Int32, Ptr{Float64}, Ptr{Int32}, Ptr{Int64}))
        task.usercallbackfunc = nothing
        finalizer(deletetask,task)

        r = @msk_ccall(putcallbackfunc, Cint, (Ptr{Nothing}, Ptr{Nothing}, Any), task.task, task.callbackfunc, task)
        if r != MSK_RES_OK.value
            throw(MosekError(r,getlasterrorx(t)))
        end

        task
    end

    function Task(t::Task)
        temp = Array{Ptr{Nothing}}(undef,1)
        res = @msk_ccall(clonetask, Int32, (Ptr{Nothing}, Ptr{Nothing}), t.task, temp)

        if res != MSK_RES_OK.value
            throw(MosekError(res,""))
        end

        task = new(t.env,temp[1],false,nothing,nothing,nothing,nothing,nothing)
        task.callbackfunc = @cfunction(msk_info_callback_wrapper, Cint, (Ptr{Nothing}, Ptr{Nothing}, Int32, Ptr{Float64}, Ptr{Int32}, Ptr{Int64}))
        task.usercallbackfunc = nothing

        finalizer(deletetask,task)

        r = @msk_ccall(putcallbackfunc, Cint, (Ptr{Nothing}, Ptr{Nothing}, Any), task.task, task.callbackfunc, task)
        if r != MSK_RES_OK.value
            throw(MosekError(r,getlasterrorx(t)))
        end

        task
    end

    function Task(t::Ptr{Nothing},borrowed::Bool)
        task = new(msk_global_env,t,borrowed,nothing,nothing,nothing,nothing,nothing)
        task.callbackfunc = @cfunction(msk_info_callback_wrapper, Cint, (Ptr{Nothing}, Ptr{Nothing}, Int32, Ptr{Float64}, Ptr{Int32}, Ptr{Int64}))
        task.usercallbackfunc = nothing

        finalizer(deletetask,task)

        r = @msk_ccall(putcallbackfunc, Cint, (Ptr{Nothing}, Ptr{Nothing}, Any), task.task, task.callbackfunc, task)
        if r != MSK_RES_OK.value
            throw(MosekError(r,getlasterrorx(t)))
        end

        task
    end
end
const MSKtask = Task
const MSKenv  = Env

# ------------
# API wrappers
# ------------
# TODO: Support other argument
"""
    makeenv()

Create a MOSEK environment.
"""
function makeenv()
    temp = Array{Ptr{Nothing}}(undef, 1)
    res = @msk_ccall(makeenv, Int32, (Ptr{Ptr{Nothing}}, Ptr{UInt8}), temp, C_NULL)
    if res != 0
        # TODO: Actually use result code
        error("MOSEK: Error creating environment")
    end
    Env(temp[1],nothing)
end

"""
    makeenv(func::Function)

Create a MOSEK environment for use with `do`-syntax.
"""
function makeenv(func::Function)
    temp = Array{Ptr{Nothing}}(undef,1)
    res = @msk_ccall(makeenv, Int32, (Ptr{Ptr{Nothing}}, Ptr{UInt8}), temp, C_NULL)
    if res != 0
        # TODO: Actually use result code
        error("MOSEK: Error creating environment")
    end
    env = Env(temp[1],nothing)

    try
        func(env)
    finally
        deleteenv(env)
    end
end

# Note on initialization of msk_global_env:
#
#  When loading Mosek from source this works fine, but when loading
#  precompiled module, makeenv() is not called (and some garbage
#  value is put in msk_global_env). It appears that static
#  initializers must be called from __init__(). That is called a bad solution here:
#    https://github.com/JuliaLang/julia/issues/12010
msk_global_env = makeenv() :: Env
__init__() = (global msk_global_env = makeenv())

"""
    maketask(;env::Env = msk_global_env, filename::String = "")

Create a task. If `filename` is not 0-length, initialize the task from this file.

    maketask(func :: Function;env::Env = msk_global_env, filename::String = "")

Create a task. If `filename` is not 0-length, initialize the task from
this file. The `func` parameter is a function

```julia
func(task :: Task) :: Any
```

This can be used with the `do`-syntax:

```julia
maketask(env) do task
    readdata(task,"MyFile.task")
end
```
    maketask(task::Task)

Create a clone of `task`.
"""
function maketask end
function maketask(;env::Env = msk_global_env, filename::String = "")
    t = Task(env)
    if length(filename) > 0
        try
            readdata(t,filename)
        catch e
            deletetask(t)
            rethrow()
        end
    end
    t
end

function maketask(func :: Function;env::Env = msk_global_env, filename::String = "")
    t = Task(env)
    try
        if length(filename) > 0 readdata(t,filename) end
        func(t)
    finally
        deletetask(t)
    end
end

maketask(task::Task) = Task(task)

function maketask_ptr(t::Ptr{Nothing},borrowed::Bool)
    Task(t,borrowed)
end

"""
    deletetask(t::Task)

Destroy the task object.
"""
function deletetask(t::Task)
    if t.task != C_NULL
        if ! t.borrowed
            temp = Array{Ptr{Nothing}}(undef,1)
            temp[1] = t.task
            @msk_ccall(deletetask,Int32,(Ptr{Ptr{Nothing}},), temp)
        end
        t.task = C_NULL
    end
end

"""
    deleteenv(t::Env)

Destroy the task object.
"""
function deleteenv(e::Env)
    if e.env != C_NULL
        temp = Array{Ptr{Nothing}}(undef,1)
        temp[1] = t.env
        @msk_ccall(deleteenv,Int32,(Ptr{Ptr{Nothing}},), temp)
        e.env = C_NULL
    end
end

# function getlasterror(t::Task)
#     lasterrcode = Array{Cint}(undef,1)
#     lastmsglen = Array{Cint}(undef,1)

#     @msk_ccall(getlasterror,Cint,(Ptr{Nothing},Ptr{Cint},Cint,Ptr{Cint},Ptr{UInt8}),
#                t.task, lasterrcode, 0, lastmsglen, C_NULL)
#     lastmsg = Array{UInt8}(undef,lastmsglen[1])
#     @msk_ccall(getlasterror,Cint,(Ptr{Nothing},Ptr{Cint},Cint,Ptr{Cint},Ptr{UInt8}),
#                t.task, lasterrcode, lastmsglen[1], lastmsglen, lastmsg)
#     String(lastmsg[1:lastmsglen[1]-1])
# end

function getlasterrorx(t::Task)
    let _,_,lastmsg = getlasterror(t)
        lastmsg
    end
end

using SparseArrays

include("msk_enums.jl")
include("msk_functions.jl")
include("msk_function_ext.jl")
include("msk_callback.jl")

## General Convex optimizer has been discontinued
#include("msk_geco.jl")

include("show.jl")
include("ext_functions.jl")

#import MathProgBase
#struct MosekSolver <: MathProgBase.AbstractMathProgSolver
#    options
#end
#MosekSolver(;kwargs...) = MosekSolver(kwargs)
#export MosekSolver

#include("MosekSolverInterface.jl")

end


text:
using StructOpt
using DataStructures, LinearAlgebra, PGFPlotsX, Random

function get_problems()
    regularizers = Set([regularizer_l1, regularizer_TV])

    problems = Dict()

    nseeds = 5

    ## Lasso l1
    l1_pbs = []

    m, n, sparsity = 60, 128, 8
    Random.seed!(1234)
    xstart = rand(n) .* 10
    for seed in 1:nseeds
        pb = get_randomlasso(n, m, sparsity, reg=regularizer_l1, seed = seed)
        xorig = pb.x0

        push!(l1_pbs, (name = "pblasso_l1", pb = pb, xstart = xstart))
    end

    problems["l1_randinit"] = l1_pbs

    ## Lasso l1 zero init
    l1_zeroinit_pbs = []

    for pb in l1_pbs
        push!(l1_zeroinit_pbs, (name = "pblasso_l1", pb = pb.pb, xstart = zeros(n)))
    end

    problems["l1_zeroinit"] = l1_zeroinit_pbs

    return problems
end

function get_algorithms()
    algorithms = []

    # αuser = 0.001
    itmax = 4e4
    printstep = 1e3
    push!(algorithms, (
        name="FISTA",
        updatefunc=extra_FISTA,
        params=Dict(
            :saveiter => false,
            # :αuser => αuser,
            :printstep => printstep,
            :itmax => itmax,
            :saveiter => false,
        ),
    ))

    push!(algorithms, (
        name="ISTA",
        updatefunc=extra_ISTA,
        params=Dict(
            :saveiter => false,
            # :αuser => αuser,
            :printstep => printstep,
            :itmax => itmax,
            :saveiter => false,
        ),
    ))

    push!(algorithms, (
        name="T1",
        updatefunc=extra_CondInertia,
        params=Dict(
            :saveiter => false,
            :id_testkind=>:default,
            # :αuser => αuser,
            :printstep => printstep,
            :itmax => itmax,
            :saveiter => false,
        ),
    ))

    push!(algorithms, (
        name="T2",
        updatefunc=extra_CondPredInertia,
        params=Dict(
            :saveiter => false,
            # :αuser => αuser,
            :printstep => printstep,
            :itmax => itmax,
            :saveiter => false,
        ),
    ))

    return algorithms
end

FIGS_FOLDER = "./figs"
basename(pwd()) == "src" && (FIGS_FOLDER = joinpath("..", FIGS_FOLDER))
!ispath(FIGS_FOLDER) && mkpath(FIGS_FOLDER)

runexpnums(get_problems(), get_algorithms(), FIGS_FOLDER=FIGS_FOLDER)


text:
module FMMTrees

#using VectorBackedLists
#using AbstractTrees

#include("lists.jl")
include("trees.jl")
include("simpletrees.jl")
include("pointerbasedtrees.jl")
#include("blocktrees.jl")
#include("hmatrix.jl")
#include("cluster_tree.jl")
#include("aca.jl")

include("octree.jl")
include("leveltree.jl")
include("blocktree.jl")

end # module


text:
@mix @columns struct LocalContribution{LC}
    # Field                   | Default  | Flat  | Bounds     | Description
    localContribution::LC     | 0.5     | true  | (0.0, 1.0)  | "Contribution of local cell to the output"
end

@mix @columns struct HoodContribution{HC}
    # Field                  | Default  | Flat  | Bounds     | Description
    hoodContribution::HC     | 0.5     | true  | (0.0, 1.0)  | "Contribution of neighbors to the output"
end

@Kernel @LocalContribution @HoodContribution struct MatingPopulation{R,W} <: NeighborhoodRule{R,W} end

# PHENOTYPE and GENOTYPE are the same it's just to update frequency of element

@inline applyrule(data, rule::MatingPopulation, variable, index) = begin
    @fastmath rule.localContribution*variable + rule.hoodContribution*mean(neighborhood(rule)) / (rule.localContribution + rule.hoodContribution)
end

"""
Mating dispersal rules.
"""

@mix @columns struct Jump{JX,JY}
    # Field       | Default  | Flat  | Bounds    | Description
    jumpX::JX     | 1        | true  | (1, 100)  | "Dispersal in X"
    jumpY::JY     | 1        | true  | (1, 100)  | "Dispersal in Y"
end

@Jump struct MatingDispersal{R,W} <: ManualRule{R,W} end

@inline applyrule!(data, rule::MatingDispersal{R,W}, state, index) where {R,W} = begin
    # Ignore empty cells
    state > zero(state) || return state
    #
    jump = (rule.jumpX,rule.jumpY)
    jumpdest, is_inbounds = inbounds(jump .+ index, gridsize(data), RemoveOverflow())
    # Update spotted cell if it's on the grid
    if is_inbounds
        @inbounds data[W][jumpdest...] = state
    end
end


text:
using CSV , Tables, DataFrames
using Dates , Query , Turing

covid_data = CSV.read("data/covid-long17Mar.csv")

south_america_l = ["Brazil","Uruguay", "Argentina","Chile","Paraguay","Peru","Bolivia","Ecuador","Colombia","Venezuela","Suriname"]
covid_plot = covid_data |> @filter(_.Country ∈ south_america_l) |> DataFrame

@model exponent_mult(time_x,n_cases,country_ids) = begin
	red_ids = collect(1:length(unique(country_ids)))
	σ ~ Truncated(Cauchy(0, 1), 0, Inf)
	μ ~ Normal(0, 1)
	beta ~ Normal(0,1)
	s ~ Uniform(0,50)
 	N_countr = length(red_ids)
	a_countr = Vector{Real}(undef, N_countr)
 	a_countr ~ [Normal(μ, σ)]
	sc_parameter = [a_countr[red_ids[i]] for i = 1:N_countr]

	for j in 1:length(n_cases)
			rdeff = sc_parameter[country_ids[j]]
			n_cases[j] ~ Normal(beta*exp(rdeff*time_x[j]), s)
		end
	end

num_chains = 4
chain = mapreduce(
    c -> sample(exponent_mult(x_times, case_values, countr_id), NUTS(200, 0.65), 1000, discard_adapt=false),
    chainscat,
    1:num_chains);
describe(chain)
plot(chain)
gelmandiag(chain)

chains_new = chain[201:2500,:,:]
plot(chains_new)
corner(chains_new)

max_stars_repo_name: fargolo/covid-dash

max_issues_repo_name: fargolo/covid-dash

max_forks_repo_name: fargolo/covid-dash


text:
using Base.Test
using ODLCommonTools
import ODLCommonTools.sview
using SummationByParts  # SBP operators

push!(LOAD_PATH, abspath(joinpath(pwd(), "..")))  # get TestSystem
push!(LOAD_PATH, abspath(joinpath(pwd(), "../common")))  # get TestCommon

using PDESolver
#using Base.Test
using TestCommon
using PdePumiInterface  # common mesh interface - pumi
using EulerEquationMod
using LinearSolvers
using NonlinearSolvers   # non-linear solvers
using OptimizationInterface
using ArrayViews
import ArrayViews.view
using EulerEquationMod
using Utils
using MPI
using Input
using PETSc2
using Jacobian

#------------------------------------------------------------------------------
# define tests and tags

using TestSystem
# define tags that will be used
include("../tags.jl")

# test list
global const EulerTests = TestList()

include("test_eqn_deepcopy.jl")     # note: eqn gets written random values to it, so anything that
include("test_empty.jl")
include("test_input.jl")
include("test_complexify.jl")
include("test_lowlevel.jl")

include("test_dg.jl")
#include("test_simplemesh.jl")
include("test_GLS3.jl")
include("test_modes.jl")
include("test_3d.jl")
include("test_bc.jl")
include("test_interp.jl")
include("test_jac.jl")
include("test_functionals.jl")
include("test_adjoint.jl")
include("test_flux.jl")
include("test_ESS.jl")
include("test_curvilinear.jl")
include("test_rk4.jl")
include(joinpath("./convergence/runtests.jl"))
include("Utils.jl")
include("test_parallel_serialpart.jl")
include("test_homotopy.jl")
include("test_staggered.jl")
include("test_checkpoint.jl")
include("test_coordsfd.jl")
include("test_shock_capturing.jl")
include("test_lps.jl")

#------------------------------------------------------------------------------
# run tests
@testset "----- Running Euler tests -----" begin

  runTestSystem(EulerTests, solvePDE, ARGS)
end

#------------------------------------------------------------------------------
# cleanup


text:
first_3_of_6(v::SVector{6,T}) where {T} = SVector{3,T}(v[1], v[2], v[3])
last_3_of_6(v::SVector{6,T}) where {T} = SVector{3,T}(v[4], v[5], v[6])

function getTop(a::SMatrix{4,4,T,16}) where {T}
    return @SMatrix [
    a[1]  a[5]  a[9]  a[13];
    a[2]  a[6]  a[10] a[14];
    a[3]  a[7]  a[11] a[15]
    ]
end

unPad(a::SVector{4,T}) where {T} = SVector{3,T}(a[1], a[2], a[3])
unPad(a::SMatrix{1,4,T,4}) where {T} = SVector{3,T}(a[1], a[2], a[3])
onePad(a::SVector{3,T}) where {T} = SVector{4,T}(a[1], a[2], a[3], one(T))
zeroPad(a::SVector{3,T}) where {T} = SVector{4,T}(a[1], a[2], a[3], zero(T))

@inline mul_then_un_pad(m::SMatrix{4,4,T1,16}, v::SVector{4,T2}) where {T1,T2} = unPad(m * v)
@inline one_pad_then_mul(m::SMatrix{4,4,T1,16}, v::SVector{3,T2}) where {T1,T2} = m * onePad(v)

function weightPoly(p1::SVector{N,T}, p2::SVector{N,T}, w1::T, w2::T) where {N,T}
  sum_weight = w1 - w2
  c1 = w1 / sum_weight
  c2 = w2 / sum_weight
  return c1 * p2 - c2 * p1
end

function make_pd_gains(time_settle::Float64, damping::Float64=1.0)
    w_n = 2 * pi / time_settle
    p = w_n^2
    d = 2 * damping * w_n
    return p, d
end

"""
Creates a random positive definite SMatrix for a size.
"""
function rand_pd(n::Int64)
    USV = svd(randn(SMatrix{n,n,Float64,n^2}))
    return Hermitian(USV.U * Diagonal(SVector(USV.S)) * USV.U')
end


text:
"""
    split_field(mesh, ξ)

Split the field `ξ` into a vector containing a `view` of each compartment.
"""
function split_field(mesh, ξ)
    ncompartment = length(mesh.points)
    npoint_cmpts = size.(mesh.points, 2)
    inds = [0; cumsum(npoint_cmpts[1:end])]
    [@view(ξ[1+inds[i]:inds[i+1]]) for i = 1:ncompartment]
end


text:
"""
Quadrature information

Most spectral-element operations require calculation of integrals and derivatives.
The `QuadType` structure stores the required information when calculating integrals
and derivatives
"""
type QuadType{T<:Number}
    Q::Int
    z::Vector{T}
    w::Vector{T}
    D::Matrix{T}
end

"""
Constructs quadrature rules for Gauss-type quadrature.
"""
function QuadType{T<:Number,QT<:QUADRATURE_TYPE}(q::Int, ::Type{QT}=GLJ, ::Type{T}=Float64, a=0, b=0)
    z = Jacobi.qzeros(QT, q, a, b, T)
    w = Jacobi.qweights(QT, z, a, b)
    D = Jacobi.qdiff(QT, z, a, b)
    QuadType(q, z, w, D)
end

"""
Creates the ith Lagrange polynomial from a vector of nodes

Explicit construction of Lagrange polynomials from the nodes.
This function builds a `Poly` object.

 * `i` Index specifying which Lagrange polynomial should be returned.
 * `z` Vector containing the nodes of the Lagrange polynomials.
"""
function lagrange_poly{T<:Number}(i::Integer, z::AbstractVector{T})
    np = length(z)

    y = zeros(T, np-1)

    for j = 1:(i-1)
        y[j] = z[j]
    end
    for j = (i+1):np
        y[j-1] = z[j]
    end
    den = one(T)
    for j = 1:(np-1)
        den = den * (z[i] - y[j])
    end
    p = poly(y)
    return p/den
end

"""
Build Lagrange basis from any nodes

Probably not the fastest algorithm. But this won't be seriously used.
"""
function QuadType{T<:Number}(x::AbstractVector{T})  # Generic numbers:
    np = length(x)
    z = zeros(T,np)
    for i = 1:np
        z[i] = x[i]
    end

    w = zeros(T, np)
    D = zeros(T, np, np)

    for i = 1:np
        p = lagrange_poly(i, z)
        pint = polyint(p)
        pder = polyder(p)
        w[i] = polyval(pint, one(T)) - polyval(pint, -one(T))
        # Compute the rows of the derivative matrix
        for k = 1:np
            D[k,i] = polyval(pder, z[k])
        end
    end
    QuadType(np, z, w, D)
end


text:
abstract type AbstractLShaped end

StochasticPrograms.num_subproblems(lshaped::AbstractLShaped) = lshaped.num_subproblems
num_cuts(lshaped::AbstractLShaped) = lshaped.data.num_cuts
num_iterations(lshaped::AbstractLShaped) = lshaped.data.iterations
tolerance(lshaped::AbstractLShaped) = lshaped.parameters.τ

# Initialization #
# ======================================================================== #
function initialize!(lshaped::AbstractLShaped)
    # Initialize progress meter
    lshaped.progress.thresh = lshaped.parameters.τ
    # Initialize subproblems
    initialize_subproblems!(lshaped, scenarioproblems(lshaped.structure, 2))
    # Prepare the master optimization problem
    prepare_master!(lshaped)
    # Initialize regularization policy
    initialize_regularization!(lshaped)
    # Finish initialization
    finish_initilization!(lshaped)
    return nothing
end
# ======================================================================== #

# Functions #
# ======================================================================== #
function prepare_master!(lshaped::AbstractLShaped)
    # Check sense first
    sense = MOI.get(lshaped.master, MOI.ObjectiveSense())
    if sense == MOI.FEASIBILITY_SENSE
        lshaped.data.no_objective = true
        # Use min-sense during L-shaped procedure
        MOI.set(lshaped.master, MOI.ObjectiveSense(), MOI.MIN_SENSE)
        F = AffineDecisionFunction{Float64}
        MOI.set(lshaped.master, MOI.ObjectiveFunction{F}(), zero(F))
    else
        # Cache the objective function
        F = MOI.get(lshaped.master, MOI.ObjectiveFunctionType())
        lshaped.data.master_objective = MOI.get(lshaped.master, MOI.ObjectiveFunction{F}())
    end
    # Initialize the required number of master variables. Use
    # MOI.VariableIndex(0) as an undef value
    append!(lshaped.master_variables, fill(MOI.VariableIndex(0), num_thetas(lshaped)))
end

function add_master_variable!(lshaped::AbstractLShaped, index::Integer)
    F = MOI.get(lshaped.master, MOI.ObjectiveFunctionType())
    master_variable = MOI.add_variable(lshaped.master)
    # Set name
    MOI.set(lshaped.master, MOI.VariableName(), master_variable,
            add_subscript("θ", index))
    # Get sense
    sense = MOI.get(lshaped.master, MOI.ObjectiveSense())
    coeff = sense == MOI.MIN_SENSE ? 1.0 : -1.0
    # Add to objective
    MOI.modify(lshaped.master, MOI.ObjectiveFunction{F}(), MOI.ScalarCoefficientChange(master_variable, coeff))
    lshaped.master_variables[index] = master_variable
    return nothing
end

function remove_cut_constraints!(lshaped::AbstractLShaped)
    # Decrease count
    lshaped.data.num_cuts -= length(lshaped.cut_constraints)
    # Remove cuts
    for ci in lshaped.cut_constraints
        if !iszero(ci.value)
            MOI.delete(lshaped.master, ci)
        end
    end
    empty!(lshaped.cut_constraints)
    return nothing
end

function restore_master!(lshaped::AbstractLShaped)
    # Remove cut constraints
    remove_cut_constraints!(lshaped)
    # Remove master variables
    for var in lshaped.master_variables
        if !iszero(var.value)
            MOI.delete(lshaped.master, var)
        end
    end
    empty!(lshaped.master_variables)
    # Remove any regularization terms
    restore_regularized_master!(lshaped)
    if lshaped.data.no_objective
        # Re-set FEASIBILITY_SENSE
        MOI.set(lshaped.master, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    else
        # Re-add original objective
        @unpack master_objective = lshaped.data
        F = typeof(master_objective)
        MOI.set(lshaped.master, MOI.ObjectiveFunction{F}(), master_objective)
    end
    return nothing
end

function filter_cuts!(lshaped::AbstractLShaped, list)
    # Nothing to do if constraints do not match
    return nothing
end

function filter_cuts!(lshaped::AbstractLShaped, list::Vector{<:CutConstraint})
    # Filter out the cut constraints
    for cut in optimizer.lshaped.cut_constraints
        i = something(findfirst(isequal(c), list), 0)
        if !iszero(i)
            MOI.deleteat!(list, i)
        end
    end
    return nothing
end

function active_model_objectives(lshaped::AbstractLShaped)
    return map(lshaped.master_variables) do var
        var.value != 0
    end
end

function update_solution!(lshaped::AbstractLShaped)
    ncols = num_decisions(lshaped.structure)
    nb = num_thetas(lshaped)
    lshaped.x .= MOI.get.(lshaped.master, MOI.VariablePrimal(), all_decisions(lshaped.decisions))
    θs = map(lshaped.master_variables) do vi
        if vi.value == 0
            -1e10
        else
            MOI.get(lshaped.master, MOI.VariablePrimal(), vi)
        end
    end
    set_model_objectives(lshaped, θs)
    return nothing
end

function decision(lshaped::AbstractLShaped, index::MOI.VariableIndex)
    i = something(findfirst(i -> i == index, all_decisions(lshaped.decisions)), 0)
    if iszero(i)
        throw(MOI.InvalidIndex(index))
    end
    return decision(lshaped)[i]
end

function evaluate_first_stage(lshaped::AbstractLShaped, x::AbstractVector)
    model = lshaped.master
    # Get objective
    @unpack master_objective = lshaped.data
    # Evaluate objective
    decisions = all_decisions(lshaped.decisions)
    obj_val = MOIU.eval_variables(master_objective) do vi
        if vi in decisions
            # Only evaluate decision
            x[vi.value]
        else
            0.0
        end
    end
    # Return value
    return obj_val
end

function current_objective_value(lshaped::AbstractLShaped)
    # Get sense
    sense = MOI.get(lshaped.master, MOI.ObjectiveSense())
    correction = sense == MOI.MIN_SENSE ? 1.0 : -1.0
    # Return sense-corrected value
    return evaluate_first_stage(lshaped, current_decision(lshaped)) +
        correction * sum(subobjectives(lshaped))
end

function calculate_estimate(lshaped::AbstractLShaped)
    # Get sense
    sense = MOI.get(lshaped.master, MOI.ObjectiveSense())
    correction = sense == MOI.MIN_SENSE ? 1.0 : -1.0
    # Return sense-corrected value
    return evaluate_first_stage(lshaped, lshaped.x) +
        correction * sum(model_objectives(lshaped))
end

function log!(lshaped::AbstractLShaped; optimal = false, status = nothing)
    @unpack Q, θ, iterations = lshaped.data
    @unpack keep, offset, indent = lshaped.parameters
    # Early termination log
    if status != nothing && lshaped.parameters.log
        lshaped.progress.thresh = Inf
        lshaped.progress.printed = true
        val = if status == MOI.INFEASIBLE
            Inf
        elseif status == MOI.DUAL_INFEASIBLE
            -Inf
        else
            0.0
        end
        ProgressMeter.update!(lshaped.progress, val,
                              showvalues = [
                                  ("$(indentstr(indent))Objective", val),
                                  ("$(indentstr(indent))Early termination", status),
                                  ("$(indentstr(indent))Number of cuts", num_cuts(lshaped)),
                                  ("$(indentstr(indent))Iterations", iterations)
                              ], keep = keep, offset = offset)
        return nothing
    end
    # Value update
    push!(lshaped.Q_history, Q)
    push!(lshaped.θ_history, θ)
    lshaped.data.iterations += 1
    # Log
    log_regularization!(lshaped)
    if lshaped.parameters.log
        current_gap = optimal ? 0.0 : gap(lshaped)
        ProgressMeter.update!(lshaped.progress, current_gap,
                              showvalues = [
                                  ("$(indentstr(indent))Objective", objective_value(lshaped)),
                                  ("$(indentstr(indent))Gap", current_gap),
                                  ("$(indentstr(indent))Number of cuts", num_cuts(lshaped)),
                                  ("$(indentstr(indent))Iterations", iterations)
                              ], keep = keep, offset = offset)
    end
    return nothing
end

function log!(lshaped::AbstractLShaped, t::Integer; optimal = false, status = nothing)
    @unpack Q,θ,iterations = lshaped.data
    @unpack keep, offset, indent = lshaped.parameters
    lshaped.Q_history[t] = Q
    lshaped.θ_history[t] = θ
    if status != nothing && lshaped.parameters.log
        lshaped.progress.thresh = Inf
        lshaped.progress.printed = true
        val = if status == MOI.INFEASIBLE
            Inf
        elseif status == MOI.DUAL_INFEASIBLE
            -Inf
        else
            0.0
        end
        ProgressMeter.update!(lshaped.progress, val,
                              showvalues = [
                                  ("$(indentstr(indent))Objective", val),
                                  ("$(indentstr(indent))Early termination", status),
                                  ("$(indentstr(indent))Number of cuts", num_cuts(lshaped)),
                                  ("$(indentstr(indent))Iterations", iterations)
                              ], keep = keep, offset = offset)
        return nothing
    end
    log_regularization!(lshaped,t)
    if lshaped.parameters.log
        current_gap = optimal ? 0.0 : gap(lshaped)
        ProgressMeter.update!(lshaped.progress, current_gap,
                              showvalues = [
                                  ("$(indentstr(indent))Objective", objective_value(lshaped)),
                                  ("$(indentstr(indent))Gap", current_gap),
                                  ("$(indentstr(indent))Number of cuts", num_cuts(lshaped)),
                                  ("$(indentstr(indent))Iterations", iterations)
                              ], keep = keep, offset = offset)
    end
    return nothing
end

function indentstr(n::Integer)
    return repeat(" ", n)
end

function check_optimality(lshaped::AbstractLShaped)
    @unpack τ = lshaped.parameters
    @unpack θ = lshaped.data
    return θ > -Inf && gap(lshaped) <= τ
end
# ======================================================================== #

# Cut functions #
# ======================================================================== #
active(lshaped::AbstractLShaped, hyperplane::AbstractHyperPlane) = active(hyperplane, decision(lshaped), tolerance(lshaped))
active(lshaped::AbstractLShaped, cut::HyperPlane{OptimalityCut}) = optimal(cut, decision(lshaped), model_objectives(lshaped)[cut.id], tolerance(lshaped))
active(lshaped::AbstractLShaped, cut::AggregatedOptimalityCut) = optimal(cut, decision(lshaped), sum(model_objectives(lshaped)[cut.ids]), tolerance(lshaped))
satisfied(lshaped::AbstractLShaped, hyperplane::AbstractHyperPlane) = satisfied(hyperplane, decision(lshaped), tolerance(lshaped))
satisfied(lshaped::AbstractLShaped, cut::HyperPlane{OptimalityCut}) = satisfied(cut, decision(lshaped), model_objectives(lshaped)[cut.id], tolerance(lshaped))
satisfied(lshaped::AbstractLShaped, cut::AggregatedOptimalityCut) = satisfied(cut, decision(lshaped), sum(model_objectives(lshaped)[cut.ids]), tolerance(lshaped))
violated(lshaped::AbstractLShaped, hyperplane::AbstractHyperPlane) = !satisfied(lshaped, hyperplane)
gap(lshaped::AbstractLShaped, hyperplane::AbstractHyperPlane) = gap(hyperplane, decision(lshaped))
gap(lshaped::AbstractLShaped, cut::HyperPlane{OptimalityCut}) = gap(cut, decision(lshaped), model_objectives(lshaped)[cut.id])
gap(lshaped::AbstractLShaped, cut::AggregatedOptimalityCut) = gap(cut, decision(lshaped), sum(model_objectives(lshaped)[cut.ids]))

function add_cut!(lshaped::AbstractLShaped, cut::AbstractHyperPlane; consider_consolidation = true, check = true)
    added = add_cut!(lshaped, cut, model_objectives(lshaped), subobjectives(lshaped), check = check)
    if consider_consolidation
        added && add_cut!(lshaped, lshaped.consolidation, cut)
    end
    return added
end

function add_cut!(lshaped::AbstractLShaped, cut::HyperPlane{OptimalityCut}, θs::AbstractVector, subobjectives::AbstractVector, Q::AbstractFloat; check = true)
    if lshaped.master_variables[cut.id].value == 0
        # Add master variable if cut is encountered for the first time
        add_master_variable!(lshaped, cut.id)
    end
    # Get the model value
    θ = θs[cut.id]
    @unpack cut_scaling = lshaped.parameters
    # Update objective
    subobjectives[cut.id] = Q
    # Check if cut gives new information
    if check && θ > -Inf && (θ + sqrt(eps()) >= Q || θ + sqrt(eps()) >= cut(lshaped.x))
        # Optimal with respect to this subproblem
        return false
    end
    # Add optimality cut
    process_cut!(lshaped, cut)
    f, set = moi_constraint(cut, lshaped.master_variables, cut_scaling)
    push!(lshaped.cut_constraints, MOI.add_constraint(lshaped.master, f, set))
    lshaped.data.num_cuts += 1
    if lshaped.parameters.debug
        push!(lshaped.cuts, cut)
    end
    return true
end
function add_cut!(lshaped::AbstractLShaped, cut::AggregatedOptimalityCut, θs::AbstractVector, subobjectives::AbstractVector, Q::AbstractFloat; check = true)
    for id in cut.ids
        if lshaped.master_variables[id].value == 0
            # Add master variable if cut is encountered for the first time
            add_master_variable!(lshaped, id)
        end
    end
    θs = θs[cut.ids]
    θ = sum(θs)
    @unpack cut_scaling = lshaped.parameters
    # Update objective
    subobjectives[cut.ids] .= Q/length(cut.ids)
    # Check if cut gives new information
    if check && θ > -Inf && (θ + sqrt(eps()) >= Q || θ + sqrt(eps()) >= cut(lshaped.x))
        # Optimal with respect to these subproblems
        return false
    end
    # Add optimality cut
    process_cut!(lshaped, cut)
    f, set = moi_constraint(cut, lshaped.master_variables, cut_scaling)
    push!(lshaped.cut_constraints, MOI.add_constraint(lshaped.master, f, set))
    lshaped.data.num_cuts += 1
    if lshaped.parameters.debug
        push!(lshaped.cuts, cut)
    end
    return true
end
add_cut!(lshaped::AbstractLShaped, cut::AnyOptimalityCut, θs::AbstractVector, subobjectives::AbstractVector, x::AbstractVector; check = true) = add_cut!(lshaped, cut, θs, subobjectives, cut(x); check = check)
add_cut!(lshaped::AbstractLShaped, cut::AnyOptimalityCut, θs::AbstractVector, subobjectives::AbstractVector; check = true) = add_cut!(lshaped, cut, θs, subobjectives, current_decision(lshaped); check = check)

function add_cut!(lshaped::AbstractLShaped, cut::HyperPlane{FeasibilityCut}, ::AbstractVector, subobjectives::AbstractVector, Q::AbstractFloat; check = true)
    # Ensure that there is no false convergence
    subobjectives[cut.id] = Q
    # Add feasibility cut
    process_cut!(lshaped, cut)
    f, set = moi_constraint(cut, lshaped.master_variables, 1.0)
    push!(lshaped.cut_constraints, MOI.add_constraint(lshaped.master, f, set))
    lshaped.data.num_cuts += 1
    if lshaped.parameters.debug
        push!(lshaped.feasibility.cuts, cut)
    end
    return true
end
add_cut!(lshaped::AbstractLShaped, cut::HyperPlane{FeasibilityCut}, θs::AbstractVector, subobjectives::AbstractVector; check = true) = add_cut!(lshaped, cut, θs, subobjectives, Inf)

function add_cut!(lshaped::AbstractLShaped, cut::HyperPlane{Infeasible}, ::AbstractVector, subobjectives::AbstractVector; check = true)
    subobjectives[cut.id] = Inf
    return true
end

function add_cut!(lshaped::AbstractLShaped, cut::HyperPlane{Unbounded}, ::AbstractVector, subobjectives::AbstractVector; check = true)
    subobjectives[cut.id] = -Inf
    return true
end


text:
"""
        destroy!(sound::SoundObject) -> Nothing

Deletes the `sound` object.
"""
function destroy! end
"""
        nchannels(sound::AbstractSound) -> Int32

Get the number of channels used by the `sound`.
"""
function nchannels end
"""
        samplerate(sound::SoundObject) -> Int32

Get the `sample rate` of the sound.

The sample rate is the number of samples played per second. The higher,
the better the quality (for example, 44100 samples/s is CD quality).
"""
function samplerate end
"""
        duration(sound::SoundObject) -> Int64

Get the total `duration` of the sound in microseconds.
"""
function duration end
"""
        nsamples(sound | buffer) -> Int32

Get the number of samples stored in the `buffer`.
"""
function nsamples end
"""
        getsamples(sound | buffer) -> Vector{Int16}

Get the array of audio samples stored in the `buffer`.

The format of the returned samples is 16 bits signed
integer (Int16). The total number of samples in
this array is given by the [`nsamples`](@ref) function.
"""
function getsamples end

"""
        play!(sound::AbstractSound) -> Nothing

Start or resume playing the `sound`.
"""
function play! end
"""
        pause!(sound::AbstractSound) -> Nothing

Pause the `sound`.

This function pauses the sound if it was playing,
otherwise (sound already paused or stopped) it has
no effect.
"""
function pause! end
"""
        stop!(sound::AbstractSound) -> Nothing

Stop playing the `sound`.

This function stops the sound if it was playing or paused,
and does nothing if it was already stopped. It also resets
the playing position (unlike [`pause`](@ref)).
"""
function stop! end

"""
        isplaying(sound::AbstractSound) -> Bool

Returns true if the `sound` is playing, `false` otherwise.
"""
function isplaying end
"""
        ispaused(sound::AbstractSound) -> Bool

Returns true if the `sound` is paused, `false` otherwise.
"""
function ispaused end
"""
        isstopped(sound::AbstractSound) -> Bool

Returns true if the `sound` is stopped, `false` otherwise.
"""
function isstopped end

"""
        timepos!(sound::AbstractSound, offset::Integer) -> Nothing

Change the current playing position of the `sound`.

The playing position can be changed when the sound is either paused or playing.
Changing the playing position when the sound is stopped has no effect, since
playing the sound will reset its position.
"""
function timepos! end
"""
        timepos(sound::AbstractSound) -> Int64

Get the current playing position of the `sound` in microseconds.
"""
function timepos end
"""
        volume!(sound::AbstractSound, volume::Real) -> Nothing

Set the `volume` of the `sound`.

The volume is a value between ``0`` (mute) and ``100`` (full volume). The default value for the volume is ``100``.
"""
function volume! end
"""
        volume(sound::AbstractSound) -> Float32[0:100]

Get the volume of the sound.
"""
function volume end
"""
        pitch!(sound::AbstractSound) -> Float32

Set the `pitch` of the `sound`.

The pitch represents the perceived fundamental frequency of a sound;
thus you can make a sound more acute or grave by changing its pitch.
A side effect of changing the pitch is to modify the playing speed
of the sound as well. The default value for the pitch is 1.
"""
function pitch! end
"""
        pitch(sound::AbstractSound) -> Float32

Get the pitch of the `sound`.
"""
function pitch end
"""
        loop!(sound::AbstractSound, loop::Bool) -> Nothing

Set whether or not the sound should loop after reaching the end.
"""
function loop! end
"""
        islooping(sound::AbstractSound) -> Bool

Returns `true` if the `sound` is in `loop` mode, `false` otherwise.
"""
function islooping end
"""
        getstatus(sound::AbstractSound) -> Symbol

Get the current status of the `sound` (stopped, paused, playing)
"""
function getstatus end


text:
@model function model_deaths(
    θ,     # [Bool] if `false`, will only compute what's needed to `observe` but not more
	::Type{TV} = Vector{Float64},
	::Type{V}  = Float64;
	predict = false,
	vectorize = true,
	context   = :inference
) where {TV, V}

    @unpack num_observations, num_total_days, num_regions, num_rt_steps, deathmodel, init_infected, num_covariates, invlink, rwscale, include_early_dynamic = θ
      # If we don't want to predict the future, we only need to compute up-to time-step `num_obs_countries[m]`
    num_time_steps = num_observations

    ############# 2.) time varying reproduction number

    R0_std     ~ truncated(Normal(.5, .2), 0.01, Inf)
    R0s        ~ filldist(truncated(Normal(0.8, R0_std), 0., 4.5), num_regions)
    σ_rt       ~ truncated(Normal(0.3*rwscale, .02*rwscale), 0, .5*rwscale)

    n = first( num_rt_steps )
    latent_Rts_z ~ filldist( CovidSurvey.RandomWalk( n ), num_regions)
    latent_Rts = CovidSurvey.rescale( latent_Rts_z, σ_rt, invlink.(R0s) )

    Rts = TV[TV(undef, num_time_steps[m]) for m in 1:num_regions]
    if num_covariates > 0
        grouped_effect ~ filldist( Laplace(0,0.2), num_covariates)
        effect_std ~ filldist( GammaMeanCv(0.1, 0.5), num_covariates)
        effects_z ~ filldist( MvNormal( num_covariates, 1.), num_regions)
        effects = [ effects_z[:,m] .* effect_std .+ grouped_effect for m in 1:num_regions]

        # grouped_effect ~ filldist( Exponential(0.2), num_covariates)
        # effect_std ~ filldist( GammaMeanCv(0.1, 0.5), num_covariates)
        # pooled_effects_distr_raw = Normal.( grouped_effect, effect_std )
        # pooled_effects_distr = arraydist( truncated.( pooled_effects_distr_raw, zero(V), V(Inf) ) )
        # effects ~ filldist( pooled_effects_distr, num_regions )

    	if θ.semiparametric
    	    CovidSurvey.semiparametric!(Rts, θ, predict, latent_Rts, R0s, σ_rt, effects)
    	else
            CovidSurvey.mixed!(Rts, θ, predict, latent_Rts, R0s, σ_rt, effects)
    	end
    else
        CovidSurvey.random_walks!(Rts, θ, predict, latent_Rts, R0s, σ_rt)
    end

    ############ 3.) infection dynamics
    # τ  ~ Exponential(100.) # `Exponential` has inverse parameterization of the one in Stan
    ys ~ arraydist(Exponential.(init_infected))
    # ys = fill(100, num_regions)

    newly_infecteds      = TV[TV(undef, num_time_steps[m]) for m in 1:num_regions]
    CovidSurvey.infections!(newly_infecteds, θ, ys, Rts)

    ########### 4.) derive observables
    μ_i2d ~ truncated(Normal(21.82, 1.01), 18, 24)
    σ_i2d ~ truncated(Normal(4.26, 5.18), 1, 14)
    # idr   = 0.8 / 100
    idr   ~ truncated(Normal(8/1000, 1/1000), 1/1000, 15/1000)
    ϕ_d   ~ truncated(Normal(50, 10), 20, Inf)

    ### day of week effect
    weekdayeffect_simplex ~ filldist(Dirichlet([7,7,7,7,7,3,1]), num_regions)
    weekdayeffect = TV[TV(undef, 7) for m in 1:num_regions]
    CovidSurvey.weekdayeffect!(weekdayeffect, weekdayeffect_simplex)

    ### holiday effect
    holidayeffect ~ filldist(Beta(1,1), num_regions)

    expected_daily_deaths    = TV[TV(undef, num_time_steps[m]) for m in 1:num_regions]
    deaths_observation_model = CovidSurvey.WeekdayHolidayObsModel2(
        deathmodel, μ_i2d, σ_i2d, idr, ϕ_d, expected_daily_deaths, weekdayeffect, holidayeffect, vectorize
    )

    # deaths_observation_model = CovidSurvey.SimpleObsModel2(deathmodel, μ_i2d, σ_i2d, idr, ϕ_d, expected_daily_deaths, vectorize)
    CovidSurvey.expected!(deaths_observation_model, newly_infecteds)

    ########### 4.) compare model to observations
    ## 4.1) observe deaths

    if context == :inference
		ℓ = CovidSurvey.logpdf(deaths_observation_model, θ.deaths)
		Turing.@addlogprob! ℓ
		return (
			newly_infecteds = newly_infecteds,
			expected_daily_deaths = expected_daily_deaths,
			Rts = Rts,
		)
	else
		if num_covariates > 0
			return (
				newly_infecteds = newly_infecteds,
				expected_daily_deaths = expected_daily_deaths,
				Rts = Rts,
				predicted_daily_deaths = CovidSurvey.prediction(deaths_observation_model), #effects = effects,
				grouped_effects = [grouped_effects for m in 1:num_regions],
				effects = effects
			)
		else
			return (
				newly_infecteds = newly_infecteds,
				expected_daily_deaths = expected_daily_deaths,
				Rts = Rts,
				predicted_daily_deaths = CovidSurvey.prediction(deaths_observation_model) #effects = effects,
				#effective_Rt = effective_Rts
			)
		end
	end
end

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
# INFORMATION
# This file demonstrates how the benchmarks for our paper were done.
# I have removed most other benchmark and plotting files from the example
# folder to keep things clean. There were scripts to retrospectively
# add and combine results from different runs, which can be avoided
# with careful planning.

# I have used `DataFrames.jl` to setup the different problem and model configurations.
# The folder `benchmark_helpers` contains helpers to work with data frames.
# Each test run is defined by a set of "features" and "dependent features".
# Dependent features depend on a subset of other features, e.g, the number of objectives
# might depend on the number of decision variables.

# There is some basic parallelization by using a `@threads` for-loop.

# ----------------------------------------------------------

# Activate benchmark environment to have all dependencies
using Pkg;
Pkg.activate(@__DIR__)

# Parse command line arguments
#
# args will have fields
# `outdir`; default = joinpath( ENV["HOME"], "MORBIT_BENCHMARKS" )
# `filename`; default = ""
# `runs-per-setting`; default = Threads.nthreads()
# `resume-from`; default = ""
include( "benchmark_helpers/utils.jl")
args = parse_commandline();

# Specific to our use case:
using Morbit;                   # our algorithm we test
using MultiObjectiveProblems;   # test problems (zdt, dtlz etc.)
using HaltonSequences           # for deterministic starting points

# Optional : Set logging level
using Logging;
logger = ConsoleLogger(stdout, Logging.Warn)
global_logger(logger);

using Dates        # for printing below

# Or logging to file; `close(io)` when done!
#=
using Logging;
using Dates;
io = open(
    joinpath(args["outdir"],
    string( Dates.format(now(), "d_u_Y__H_M_S"), "_log.txt") ),
    "w+"
);
logger = ConsoleLogger(io, Logging.Warn);
global_logger(logger);
=#

# Include helpers for DataFrames
include("benchmark_helpers/saving.jl")
include("benchmark_helpers/setting_helpers.jl");

# Preparations Done!!!

#%% ---------------------------------------------------------------------------- #

# Define "features" and their values to test for

# Small example where we simply test the performance of two model types for two problems

features = Dict(
    "method" => [:steepest_descent, :ps],
    "model" => ["cubic", "TP1", "LP1", "LP2"],
    "problem_str" => ["ZDT1", "ZDT2", "ZDT3", "DTLZ1", "DTLZ6"],
    "n_vars" => collect(2:15),
);
num_runs = args["runs-per-setting"];

# We perform `num_runs` runs for each possible combination of
# feature values. Each run should have a different starting point "x0", but the problems
# could have different box constraints. Hence the feature "x0" is **dependent** on the
# feature "problems"

# 1) define little helper to get start points from problem string
# note that its keyword arguments correspond to the feature(s)
# it depends on!
function get_test_problem(; kwargs... )
    n_vars = kwargs[:n_vars];
    prob_symb = Symbol(kwargs[:problem_str]);
    prob_type = @eval $prob_symb;
    if prob_type <: ZDT
        test_problem = prob_type( n_vars )
    elseif prob_type <: DTLZ
        k = min( n_vars - 1, 5 );
        test_problem = prob_type(
            n_vars,
            n_vars - k + 1,
            k
        )
    end
    return test_problem
end

function x0_values_function(; kwargs... )
    n_vars = kwargs[:n_vars];
    test_problem = get_test_problem(; kwargs... );
    box = constraints( test_problem );
    return [ box.lb .+ ( box.ub .- box.lb ) .* site for site in
        HaltonPoint( n_vars; length = num_runs )
    ];
end

# 2) define the actual dependent feature using a Dict and put it in a list
dependent_features = [
    Dict(
        "name" => "x0",
        "depends_on" => ["n_vars", "problem_str",],
        "values" => x0_values_function,
    ),
]

# create a table with all possible settings as rows
results = generate_all_settings( features, dependent_features );

# finally add observation columns by providing their names and types
observations = (; :n_evals => Int, :x => Vector{Float64}, :ω => Float64  );
add_observation_columns!(results, observations)

sort!(results, ["n_vars", "method", "model"]);

missing_rows = collect( 1 : size(results,1) );
# Optional: Resultume previous simulation
if length( args["resume-from"] ) > 0 && isfile( args["resume-from"] )
    old_results = load_previous_results( args["resume-from"] );
    _, missing_rows = fill_from_partial_results!( results, old_results, features, dependent_features, observations );
end

println("All in all we have to perform $(length(missing_rows)) of $(size(results,1)) tests.")

#%% ---------------------------------------------------------------------------- #

# Now setup the actual tests.
# To keep things clean we use several getter functions that take
# as keyword arguments the features defined above.

# reading the large data files in parallel does often not work,
# hence the lock for retrieving precalculated point sets with LagrangeModels
lagrange_lock = ReentrantLock();

function get_model_config(; kwargs...)
    global lagrange_lock;
    model_str = kwargs[:model];
    n_vars = kwargs[:n_vars];

    if model_str == "cubic"
        cfg = RbfConfig(
            kernel = :cubic,
            shape_parameter = 1.0,
            max_model_points = n_vars <= 10 ? (n_vars +1 )*(n_vars + 2)/2 : 2*n_vars + 1,
            θ_enlarge_1 = 2.0,
            θ_pivot = 1/4
        )
    elseif model_str == "TP1"
        cfg = TaylorConfig(
            degree = 1,
            gradients = :fdm
        )
    elseif model_str == "LP1"
        cfg = LagrangeConfig(
            degree = 1,
            optimized_sampling = true,
            Λ = 1.5,
        )
    elseif model_str == "LP2"
        cfg = LagrangeConfig(
            degree = 2,
            optimized_sampling = (n_vars <= 5),
            Λ = 1.5,
            save_path = joinpath(ENV["HOME"], "LagrangeSites", "n_vars_$(n_vars)_deg_2_Λ_15_.jld2"),
            io_lock = lagrange_lock,
        )
    end
    cfg.max_evals = n_vars * 1000;
    return cfg;
end

function get_mop(test_problem; kwargs...)
    model_cfg = get_model_config(; kwargs... );

    box = constraints( test_problem );
    objectives = get_objectives( test_problem );

    mop = MixedMOP( box.lb, box.ub );
    add_objective!( mop, objectives[1], ExactConfig(gradients=:fdm) );
    for objf in objectives[2:end]
        add_objective!( mop, objf, model_cfg );
    end
    return mop;
end

function get_algo_config( test_problem; kwargs...)
    model = kwargs[:model]
    n_vars = kwargs[:n_vars]

    AlgoConfig(;
        max_critical_loops = model == "TP1" ? 0 : 3,
        _eps_crit = .01,
        max_iter = 100,
        Δ_0 = .1,
        Δ_max = .5,
        x_tol_rel = 1e-3,
        f_tol_rel = 1e-3,
        descent_method = kwargs[:method],
        strict_acceptance_test = true,
        strict_backtracking = true,
        # some pascoletti_serafini settings
        reference_point = get_ideal_point( test_problem ) .- 1,
        ps_algo = :GN_ISRES,
        max_ps_problem_evals = 50 * (n_vars + 1),
        max_ps_polish_evals = 100 * (n_vars + 1),
        ps_polish_algo = :LD_MMA,
    )
end

# here we take arbitrary keyword arguments so that we can pass
# whole rows from the result dataframe
function perform_test(; kwargs... )

    x0 = kwargs[:x0];

    test_problem = get_test_problem(;kwargs...);
    ω_func = get_omega_function( test_problem );

    ac = get_algo_config(test_problem ;kwargs...);

    mop = get_mop(test_problem; kwargs... );

    x, _ = optimize( mop, x0; algo_config = ac)
    ω = ω_func(x);

    n_evals = maximum( o.n_evals for o ∈ mop.vector_of_objectives[2:end] );

    return (; :n_evals => n_evals, :x => x, :ω => ω);
end

# obtain functions to save results dataframe
save_func = get_save_function(args; caller_path = @__FILE__);
save_lock = ReentrantLock()     # optional, for saving when MultiThreading
save_counter = 0;
# RUN the actual tests; note the `@view` which is important to actually
# be able to store the results;

# %% start the fun
results_lock = ReentrantLock();
Threads.@threads for result_table_row in eachrow( @view(results[ missing_rows, : ]) )
    global save_counter;
    println("$(Dates.format(now(), "HH:MM:SS:sss")) : $(result_table_row)")
    try
        row_results = perform_test(; result_table_row... )

        lock(results_lock) do
            for observation in keys(observations)
                result_table_row[ observation ] = row_results[ observation ]
            end
        end

        # Option: save results every T-th iteration
        lock( save_lock ) do
            save_counter += 1;
            if save_counter % args["save-every"] == 1
                save_func( results; save_csv = false )
            end
        end
    catch e
        @error "Error for \n$(result_table_row)." exception=(e, catch_backtrace())
    end
end

#%% ---------------------------------------------------------------------------- #
# Save final results
save_func(results;save_csv = false);


text:
treehash(pkgdir::String) = Base.SHA1(Pkg.GitTools.tree_hash(pkgdir))
treehash(pkg::Module) = treehash(pkgdir(pkg))

function githttpsurl(url::AbstractString)
    m = match(LibGit2.URL_REGEX, url)
    LibGit2.git_url(scheme = "https", host = m[:host], path = m[:path])
end

create_git_cmd(config::AbstractDict; kwargs...) = create_git_cmd(pairs(config)...; kwargs...)
function create_git_cmd(config::Pair...; path = nothing)
    cmd = isnothing(path) ? ["git"] : ["git", "-C", path]
    for (n, v) in config
        push!(cmd, "-c")
        push!(cmd, "$n=$v")
    end
    Cmd(cmd)
end

gitauthurl(user, token, fullname) = "https://$(user):$(token)@github.com/$(fullname).git"

isdirty(repo_path::AbstractString) = LibGit2.isdirty(LibGit2.GitRepo(repo_path))

function repository_dispatch(
    repo::Repo,
    auth::Authorization;
    event_type::AbstractString = "dispatch",
    client_payload::Dict = Dict(),
    kwargs...,
)
    headers = Dict(
        "Accept" => "application/vnd.github.everest-preview+json",
        "Content-Type" => "application/json",
    )
    params = Dict("event_type" => event_type, "client_payload" => client_payload)
    result = GitHub.github_request(
        GitHub.DEFAULT_API,
        HTTP.post,
        "/repos/$(GitHub.name(repo))/dispatches";
        headers = headers,
        params = params,
        auth = auth,
        kwargs...,
    )
end


text:
module GibbsSeaWater

using Compat
using GibbsSeaWater_jll
using Libdl

const libgswteos = GibbsSeaWater_jll.libgswteos_path

include("gen_gswteos_h.jl")
include("gen_gswteos10.jl")
include("wrappers.jl")
include("utils.jl")

"""
wrap a gsw function with `nin` input arguments and `nout` output arguments defined
as references. All types are `Cdouble`.

Expanded functions look like these:
https://github.com/TEOS-10/GibbsSeaWater.jl/blob/9375647d97fe7cd54e11135a18e0c5447ddbed8c/src/GibbsSeaWater.jl#L17
"""
macro wr(name, nin, nout)
    return Expr(:function,
        Expr(:call,
            Symbol(:gsw_, name),
            (Expr(:(::), Symbol(:i, i), :Cdouble) for i = 1:nin)...),
        Expr(:block,
            (Expr(:(=), Symbol(:p, i), :(Ref{Cdouble}())) for i = 1:nout)...,
            Expr(:call,
                Symbol(:gsw_, name),
                (Symbol(:i, i) for i = 1:nin)..., (Symbol(:p, i) for i = 1:nout)...),
            Expr(:return,
                Expr(:tuple,
                    (Expr(:ref, Symbol(:p, i)) for i = 1:nout)...
                )
            )
        )
    )
end

#   function name                                nin nout
@wr specvol_alpha_beta 3 3
@wr specvol_first_derivatives 3 3
@wr specvol_second_derivatives 3 5
@wr specvol_first_derivatives_wrt_enthalpy 3 2
@wr specvol_second_derivatives_wrt_enthalpy 3 3
@wr rho_alpha_beta 3 3
@wr rho_first_derivatives 3 3
@wr rho_second_derivatives 3 5
@wr rho_first_derivatives_wrt_enthalpy 3 2
@wr rho_second_derivatives_wrt_enthalpy 3 3

export gsw_add_barrier
export gsw_add_mean
export gsw_adiabatic_lapse_rate_from_ct
export gsw_adiabatic_lapse_rate_ice
export gsw_alpha
export gsw_alpha_on_beta
export gsw_alpha_wrt_t_exact
export gsw_alpha_wrt_t_ice
export gsw_beta_const_t_exact
export gsw_beta
export gsw_cabbeling
export gsw_c_from_sp
export gsw_chem_potential_water_ice
export gsw_chem_potential_water_t_exact
export gsw_cp_ice
export gsw_cp_t_exact
export gsw_ct_first_derivatives
export gsw_ct_first_derivatives_wrt_t_exact
export gsw_ct_freezing
export gsw_ct_freezing_first_derivatives
export gsw_ct_freezing_first_derivatives_poly
export gsw_ct_freezing_poly
export gsw_ct_from_enthalpy
export gsw_ct_from_enthalpy_exact
export gsw_ct_from_entropy
export gsw_ct_from_pt
export gsw_ct_from_rho
export gsw_ct_from_t
export gsw_ct_maxdensity
export gsw_ct_second_derivatives
export gsw_deltasa_atlas
export gsw_deltasa_from_sp
export gsw_dilution_coefficient_t_exact
export gsw_dynamic_enthalpy
export gsw_enthalpy_ct_exact
export gsw_enthalpy_diff
export gsw_enthalpy
export gsw_enthalpy_first_derivatives_ct_exact
export gsw_enthalpy_first_derivatives
export gsw_enthalpy_ice
export gsw_enthalpy_second_derivatives_ct_exact
export gsw_enthalpy_second_derivatives
export gsw_enthalpy_sso_0
export gsw_enthalpy_t_exact
export gsw_entropy_first_derivatives
export gsw_entropy_from_ct
export gsw_entropy_from_pt
export gsw_entropy_from_t
export gsw_entropy_ice
export gsw_entropy_part
export gsw_entropy_part_zerop
export gsw_entropy_second_derivatives
export gsw_fdelta
export gsw_frazil_properties
export gsw_frazil_properties_potential
export gsw_frazil_properties_potential_poly
export gsw_frazil_ratios_adiabatic
export gsw_frazil_ratios_adiabatic_poly
export gsw_geo_strf_dyn_height
export gsw_geo_strf_dyn_height_1
export gsw_geo_strf_dyn_height_pc
export gsw_gibbs_ice
export gsw_gibbs_ice_part_t
export gsw_gibbs_ice_pt0
export gsw_gibbs_ice_pt0_pt0
export gsw_gibbs
export gsw_gibbs_pt0_pt0
export gsw_grav
export gsw_helmholtz_energy_ice
export gsw_hill_ratio_at_sp2
export gsw_ice_fraction_to_freeze_seawater
export gsw_internal_energy
export gsw_internal_energy_ice
export gsw_ipv_vs_fnsquared_ratio
export gsw_kappa_const_t_ice
export gsw_kappa
export gsw_kappa_ice
export gsw_kappa_t_exact
export gsw_latentheat_evap_ct
export gsw_latentheat_evap_t
export gsw_latentheat_melting
export gsw_linear_interp_sa_ct
export gsw_melting_ice_equilibrium_sa_ct_ratio
export gsw_melting_ice_equilibrium_sa_ct_ratio_poly
export gsw_melting_ice_into_seawater
export gsw_melting_ice_sa_ct_ratio
export gsw_melting_ice_sa_ct_ratio_poly
export gsw_melting_seaice_equilibrium_sa_ct_ratio
export gsw_melting_seaice_equilibrium_sa_ct_ratio_poly
export gsw_melting_seaice_into_seawater
export gsw_melting_seaice_sa_ct_ratio
export gsw_melting_seaice_sa_ct_ratio_poly
export gsw_nsquared
export gsw_o2sol
export gsw_o2sol_sp_pt
export gsw_pot_enthalpy_from_pt_ice
export gsw_pot_enthalpy_from_pt_ice_poly
export gsw_pot_enthalpy_ice_freezing
export gsw_pot_enthalpy_ice_freezing_first_derivatives
export gsw_pot_enthalpy_ice_freezing_first_derivatives_poly
export gsw_pot_enthalpy_ice_freezing_poly
export gsw_pot_rho_t_exact
export gsw_pressure_coefficient_ice
export gsw_pressure_freezing_ct
export gsw_pt0_cold_ice_poly
export gsw_pt0_from_t
export gsw_pt0_from_t_ice
export gsw_pt_first_derivatives
export gsw_pt_from_ct
export gsw_pt_from_entropy
export gsw_pt_from_pot_enthalpy_ice
export gsw_pt_from_pot_enthalpy_ice_poly_dh
export gsw_pt_from_pot_enthalpy_ice_poly
export gsw_pt_from_t
export gsw_pt_from_t_ice
export gsw_pt_second_derivatives
export gsw_rho_alpha_beta
export gsw_rho
export gsw_rho_first_derivatives
export gsw_rho_first_derivatives_wrt_enthalpy
export gsw_rho_ice
export gsw_rho_second_derivatives
export gsw_rho_second_derivatives_wrt_enthalpy
export gsw_rho_t_exact
export gsw_rr68_interp_sa_ct
export gsw_saar
export gsw_sa_freezing_estimate
export gsw_sa_freezing_from_ct
export gsw_sa_freezing_from_ct_poly
export gsw_sa_freezing_from_t
export gsw_sa_freezing_from_t_poly
export gsw_sa_from_rho
export gsw_sa_from_sp_baltic
export gsw_sa_from_sp
export gsw_sa_from_sstar
export gsw_sa_p_inrange
export gsw_seaice_fraction_to_freeze_seawater
export gsw_sigma0
export gsw_sigma1
export gsw_sigma2
export gsw_sigma3
export gsw_sigma4
export gsw_sound_speed
export gsw_sound_speed_ice
export gsw_sound_speed_t_exact
export gsw_specvol_alpha_beta
export gsw_specvol_anom_standard
export gsw_specvol
export gsw_specvol_first_derivatives
export gsw_specvol_first_derivatives_wrt_enthalpy
export gsw_specvol_ice
export gsw_specvol_second_derivatives
export gsw_specvol_second_derivatives_wrt_enthalpy
export gsw_specvol_sso_0
export gsw_specvol_t_exact
export gsw_sp_from_c
export gsw_sp_from_sa_baltic
export gsw_sp_from_sa
export gsw_sp_from_sk
export gsw_sp_from_sr
export gsw_sp_from_sstar
export gsw_spiciness0
export gsw_spiciness1
export gsw_spiciness2
export gsw_sr_from_sp
export gsw_sstar_from_sa
export gsw_sstar_from_sp
export gsw_t_deriv_chem_potential_water_t_exact
export gsw_t_freezing
export gsw_t_freezing_first_derivatives_poly
export gsw_t_freezing_first_derivatives
export gsw_t_freezing_poly
export gsw_t_from_ct
export gsw_t_from_pt0_ice
export gsw_thermobaric
export gsw_turner_rsubrho
export gsw_util_indx
export gsw_util_interp1q_int
export gsw_util_linear_interp
export gsw_util_sort_real
export gsw_util_xinterp1
export gsw_util_pchip_interp
export gsw_z_from_p
export gsw_p_from_z

end # module


text:
@testset "Column" begin

@testset "Constructors" begin
    v1 = [1,2,3]

    @test @inferred(Column{:A}(v1)) == Column{:A, Vector{Int}}(v1)
    @test @inferred(Column{:A, Vector{Int}}()) == Column{:A, Vector{Int}}(Vector{Int}())

    @test @Column(A = [1,2,3]) == Column{:A, Vector{Int}}(v1)
    @test (@Column A = [1,2,3]) == Column{:A, Vector{Int}}(v1)
    @test @Column(A::Vector{Int} = [1,2,3]) == Column{:A, Vector{Int}}(v1)
    @test (@Column A::Vector{Int} = [1,2,3]) == Column{:A, Vector{Int}}(v1)

    @test Column{:A, Vector{Float64}}([1,2,3]) == @Column(A=[1.0,2.0,3.0])
    @test convert(Column{:A, Vector{Float64}}, @Column(A=[1,2,3])) == @Column(A=[1.0,2.0,3.0])
end

@testset "Introspection" begin
    v1 = [1,2,3]
    col1 = Column{:A}(v1)

    @test (show(col1); println(); true)

    @test name(col1) == :A
    @test name(typeof(col1)) == :A
    @test eltype(col1) == Int
    @test eltype(typeof(col1)) == Int

    @test @inferred(rename(col1, Val{:A_new})) == @Column(A_new = v1)

    @test length(col1) == 3
    @test ncol(col1) == 1
    @test nrow(col1) == 3
    @test ndims(col1) == 1
    @test size(col1) == (3,)
    @test size(col1,1) == 3
    @test isempty(col1) == false
    @test endof(col1) == 3

end

@testset "Accessing and iterating" begin
    v1 = [1,2,3]
    col1 = Column{:A}(v1)

    @test map(x->x^2, col1) == v1.^2

    @test @inferred(getindex(col1, 2))   == 2
    @test @inferred(getindex(col1, 2:3)) == @Column A=[2,3]
    @test @inferred(getindex(col1, :))   == col1

    @test (col1_copy = copy(col1); col1_copy[3] = col1[1]; col1_copy == Column{:A}([1,2,1]))
    @test (col1_copy = copy(col1); col1_copy[3:-1:2] = col1[2:3]; col1_copy == Column{:A}([1,3,2]))
    @test (col1_copy = copy(col1); col1_copy[:] = col1[:]; col1_copy == Column{:A}([1,2,3]))

    @test (col1_copy = copy(col1); @inbounds col1_copy[3] = col1[1]; col1_copy == Column{:A}([1,2,1]))
    @test (col1_copy = copy(col1); @inbounds col1_copy[3:-1:2] = col1[2:3]; col1_copy == Column{:A}([1,3,2]))
    @test (col1_copy = copy(col1); @inbounds col1_copy[:] = col1[:]; col1_copy == Column{:A}([1,2,3]))

    @test first(col1) == 1
    @test next(col1,start(col1)) == (1, 2)
    @test last(col1) == 3

    @test (col1_copy = copy(col1); col1_copy[2] = 3; col1_copy.data == [1,3,3])
    @test (col1_copy = deepcopy(col1); col1_copy[2] = Cell{:A}(3); col1_copy.data == [1,3,3])
end

@testset "Composing and concatenating" begin
    c1 = @Cell(A = 1)
    c2 = @Cell(A = 2)
    c3 = @Cell(A = 3)
    v1 = [1,2,3]
    col1 = Column{:A}(v1)

    @test @inferred(vcat(c1)) == Column{:A}([1])
    @test @inferred(vcat(c1,c2,c3)) == col1
    @test (vcat(col1,c1)) == @Column A = [1,2,3,1] # inference fails for mixed scalar-vector vcat() (TODO fix upstream?)
    @test (vcat(c1,col1)) == @Column A = [1,1,2,3] # inference fails for mixed scalar-vector vcat() (TODO fix upstream?)
    @test @inferred(vcat(col1,col1)) == @Column A = [1,2,3,1,2,3]
end

end


text:
function modularity(g, c)
	Q = 0.
	m = 2*ne(g)
	m == 0 && return 0.
	s1 = 0
	s2 = 0
	for u in keys(g)
		for v in keys(g)
			c[u] != c[v] && continue
			s1 += haskey(g[u], v) ? 1 : 0
			s2 += length(g[u])*length(g[v])
		end
	end
	Q = s1/m - s2/m^2
	return Q
end

function range_shuffle!(r::Integer, a::AbstractVector)
    @inbounds for i=r:-1:2
        j = StatsBase.randi(i)
        a[i], a[j] = a[j], a[i]
    end
end


text:
module CropPlanner

import HTTP
import JSON
import DataFrames
using DataFramesMeta
using Pipe
import CSV
using Dates
using Unitful, Unitful.DefaultSymbols
using UnitfulUS

include("addearliestharvest!.jl")
include("addseedingdateoftransplant!.jl")
include("convertstringarray.jl")
include("converttodate.jl")
include("converttodays.jl")
include("getbasketprice.jl")
include("listavailablecropondate.jl")
include("tosiunits.jl")

export listavailablecropondate

vegdic = JSON.parse(String(read("data/vegetable_name_dictionary.json")))

# Create table of range of harvest_dates
crop_calendar = CSV.read("data/crop plan - ecoumene_crop_calendar.csv", DataFrames.DataFrame)
crop_calendar[!,:vegetable] = [get(vegdic, v, v) for v in  crop_calendar[!,:vegetable]]

rename!(crop_calendar, [ i => crop_calendar[1,i] for i in 3:ncol(crop_calendar)])
crop_calendar1 = crop_calendar[3:end,:]

possible_harvest_weeks = @pipe crop_calendar1 |>
    @where(_, :activity .== "harvest") |>
    stack(_, 3:ncol(_)) |>
    dropmissing |>
    @transform(_, week = [i.value for i in Week.(Date.(:variable, DateFormat("m/d/y")))]) |>
    @select(_, :vegetable, :activity, :week, date = :variable)

CSV.write("data/possible_harvest_weeks.csv", possible_harvest_weeks)

# Load tables
crop_notes = CSV.read("data/crop plan - jm_crop_notes.csv", DataFrames.DataFrame)
crop_notes[!,:vegetable] = [get(vegdic, v, v) for v in  crop_notes[!,:vegetable] ]
CSV.write("data/crop plan - jm_crop_notes.csv",crop_notes)

garden_plan = CSV.read("data/crop plan - jm_garden_plan.csv", DataFrames.DataFrame)
garden_plan[!,:vegetable] = [get(vegdic, v, v) for v in  garden_plan[!,:vegetable] ]
CSV.write("data/crop plan - jm_garden_plan.csv",garden_plan)

transplant_table = CSV.read("data/crop plan - jm_transplant_table.csv", DataFrames.DataFrame)
transplant_table[!,:vegetable] = [get(vegdic, v, v) for v in  transplant_table[!,:vegetable] ]
CSV.write("data/crop plan - jm_transplant_table.csv",transplant_table)

yields_and_sales = CSV.read("data/crop plan - jm_yields_and_sales.csv", DataFrames.DataFrame)
yields_and_sales[!,:vegetable] = [get(vegdic, v, v) for v in  yields_and_sales[!,:vegetable] ]
CSV.write("data/crop plan - jm_yields_and_sales.csv",yields_and_sales)

possible_harvest_weeks = CSV.read("data/possible_harvest_weeks.csv", DataFrame)

all_vegdf = DataFrame(vegetable = unique(outerjoin(crop_notes, garden_plan, transplant_table, yields_and_sales, possible_harvest_weeks, on = :vegetable,makeunique = true)[:,:vegetable]))

crop_df = outerjoin(crop_notes, garden_plan, transplant_table, yields_and_sales, on = :vegetable,makeunique = true)

crop_df_cln = @transform(crop_df,
    intensive_spacing_between_rows_in = tosiunits.(:intensive_spacing_between_rows_in),
    intensive_spacing_within_rows_in = tosiunits.(:intensive_spacing_within_rows_in),
    preferred_cultivars = convertstringarray.(:preferred_cultivars),
    dates_of_seedings = convertstringarray.(:dates_of_seeding_of_transplant),
    dates_of_transplanting = convertstringarray.(:dates_of_transplanting),
    period_of_harvest_days = convertodays.(:period_of_harvest_days),
    days_in_garden_min = convertodays.(:days_in_garden_min),
    days_in_garden_max = convertodays.(:days_in_garden_max),
    date_in = converttodate.(:date_in),
    date_out = converttodate.(:date_out)
    )

crop_df_cln[ismissing.(crop_df_cln[!,:period_of_harvest_days]), :period_of_harvest_days] .= Day(0)

addseedingdateoftransplant!(crop_df_cln)
addearliestharvest!(crop_df_cln)

show(names(crop_df_cln))
sort!(crop_df_cln, :earliest_harvest)

# Create possible_harvest_weeks for JM info
jm_harvest_range = @pipe crop_df_cln |>
    groupby(_, :vegetable) |>
    @combine(_,
        min_harvest = minimum(:earliest_harvest),
        max_harvest = maximum(:date_out)
    )

CSV.write("tmp.csv", jm_harvest_range)

end


text:
# Legit.jl -- Convert French law in LEGI format to Git & Markdown
# By: Emmanuel Raviart <emmanuel.raviart@data.gouv.fr>
#
# Copyright (C) 2015 Etalab
# https://github.com/etalab/Legit.jl
#
# The Legit.jl package is licensed under the MIT "Expat" License.

module LegitReorderCommits

using ArgParse
using Compat
using Dates
using LibGit2
using LightXML  # Not used, but needed by include("articles.jl")
import Slugify  # Not used, but needed by include("articles.jl")

include("articles.jl")

function main()
  args = parse_command_line()

  @assert ispath(args["repository"])
  @assert isdir(args["repository"])
  @assert iswritable(args["repository"])
  repository = repo_discover(args["repository"])
  @assert lookup_branch(repository, "reordered") === nothing "A branch named \"reordered\" already exists."
  @assert !is_head_unborn(repository) "HEAD points to an unborn branch (aka a branch without commit)."
  head_reference = head(repository)
  head_commit_id = target(head_reference)
  reordered_reference = lookup_ref(repository, "REORDERED")
  if reordered_reference === nothing
    reordered_reference = create_ref(repository, "REORDERED", "refs/heads/reordered")
  else
    reordered_reference = set_symbolic_target(reordered_reference, "refs/heads/reordered")
  end
  @assert lookup_branch(repository, "reordered") === nothing  # The "reordered" branch doesn't exist yet.

  commits_by_message_by_time = @compat Dict{Int64, Dict{String, Array{GitCommit}}}()
  commits_count = 1
  current_commit = lookup_commit(repository, head_commit_id)
  while true
    current_committer = committer(current_commit)
    commits_by_message = get!(commits_by_message_by_time, current_committer.time) do
      return @compat Dict{String, Array{GitCommit}}()
    end
    commits = get!(commits_by_message, message(current_commit)) do
      return GitCommit[]
    end
    push!(commits, current_commit)

    if parent_count(current_commit) == 0
      break
    end
    commits_count += 1
    current_commit = parent(current_commit, 1)
  end

  latest_commit = nothing
  latest_root_tree = nothing
  times = sort(collect(keys(commits_by_message_by_time)))
  for time in times
    commits_by_message = commits_by_message_by_time[time]
    messages = sort(collect(keys(commits_by_message)))
    for message in messages
      println("Reordering $(unix2datetime(time)) \"$message\".")
      commits = commits_by_message[message]
      section_by_git_dir_names = @compat Dict{Tuple, Union(Section, UnparsedSection)}()
      tree_builder_by_git_dir_names = @compat Dict{Tuple, GitTreeBuilder}()
      for current_commit in commits
        merge(
          repository,
          (),
          GitTree(current_commit),
          parent_count(current_commit) == 0 ? nothing : GitTree(parent(current_commit, 1)),
          latest_root_tree,
          section_by_git_dir_names,
          tree_builder_by_git_dir_names,
        )
      end

      # Generate trees (up to root tree) from tree builders.
      root_tree_id = nothing
      while !isempty(tree_builder_by_git_dir_names)
        git_dirs_names_to_build = Set(keys(tree_builder_by_git_dir_names))
        for (git_dir_names, tree_builder) in tree_builder_by_git_dir_names
          while !isempty(git_dir_names)
            git_dir_names = git_dir_names[1 : end - 1]
            pop!(git_dirs_names_to_build, git_dir_names, nothing)
          end
        end
        for git_dir_names in git_dirs_names_to_build
          tree_builder = pop!(tree_builder_by_git_dir_names, git_dir_names)
          if length(tree_builder) == 0
            tree_id = nothing
          else
            # Update README.md before building tree from tree builder.
            section = section_by_git_dir_names[git_dir_names]
            child_by_name = @compat Dict{String, Node}()
            filter!(tree_builder) do entry
              if entry.name != "README.md"
                child_by_name[entry.name] = section.child_by_name[entry.name]
              end
              return false  # Keep item in tree.
            end
            section.child_by_name = child_by_name
            blob = blob_from_buffer(repository, commonmark(section, "flat"))
            insert!(tree_builder, "README.md", Oid(blob), int(0o100644))  # FILEMODE_BLOB

            tree_id = write!(tree_builder)
          end
          if isempty(git_dir_names)
            root_tree_id = tree_id
          else
            dir_name = git_dir_names[end]
            parent_git_dir_names = git_dir_names[1 : end - 1]
            @assert !(parent_git_dir_names in git_dirs_names_to_build)
            parent_tree_builder = tree_builder_by_git_dir_names[parent_git_dir_names]
            if tree_id === nothing
              if parent_tree_builder[dir_name] !== nothing
                delete!(parent_tree_builder, dir_name)
              end
            else
              insert!(parent_tree_builder, dir_name, tree_id, int(0o40000))  # FILEMODE_TREE
            end
          end
        end
      end

      root_tree = lookup_tree(repository, root_tree_id)
      latest_commit_id = commit(repository, "REORDERED", message, author(commits[1]), committer(commits[1]),
        root_tree, (latest_commit === nothing ? [] : [latest_commit])...)
      latest_commit = lookup_commit(repository, latest_commit_id)
      latest_root_tree = root_tree
    end
  end

  # @show commits_count
  # @show length(commits_by_message_by_time)
  # @show Date(unix2datetime(times[1])) keys(commits_by_message_by_time[times[1]])
  # @show hex(times[end]) Date(unix2datetime(times[end])) keys(commits_by_message_by_time[times[end]])
end

function merge(repository::GitRepo, git_dir_names::Tuple, current_tree::GitTree, parent_tree::Union(GitTree, Nothing),
    latest_tree::Union(GitTree, Nothing), section_by_git_dir_names::Dict{Tuple, Union(Section, UnparsedSection)},
    tree_builder_by_git_dir_names::Dict{Tuple, GitTreeBuilder})
  # if Oid(first_tree) == Oid(second_tree)
  #   return first_tree
  # end
  tree_builder = get!(tree_builder_by_git_dir_names, git_dir_names) do
    return GitTreeBuilder(repository, latest_tree)
  end

  readme_entry = current_tree["README.md"]
  @assert readme_entry !== nothing
  current_section = get(section_by_git_dir_names, git_dir_names, nothing)
  if current_section === nothing
    section_by_git_dir_names[git_dir_names] = current_section = parse_section_commonmark(repository, readme_entry,"")
  else
    updated_current_section = parse_section_commonmark(repository, readme_entry, current_section.short_title)
    for (name, child_section) in updated_current_section.child_by_name
      current_section.child_by_name[name] = child_section
    end
  end
  if latest_tree !== nothing
    latest_entry = latest_tree["README.md"]
    latest_section = parse_section_commonmark(repository, latest_entry, current_section.short_title)
    for (name, child_section) in latest_section.child_by_name
      get!(current_section.child_by_name, name, child_section)  # Add section when it doesn't exist yet.
    end
  end

  for current_entry in current_tree
    current_name = current_entry.name
    if current_name == "README.md"
      continue
    end
    current_id = Oid(current_entry)
    parent_entry = parent_tree === nothing ? nothing : parent_tree[current_name]
    if parent_entry === nothing
      insert!(tree_builder, current_name, current_id, int(filemode(current_entry)))
    else
      parent_id = Oid(parent_entry)
      if current_id != parent_id
        if !isempty(git_dir_names) && (git_dir_names[end] in ("codes-en-vigueur", "codes-non-en-vigueur") ||
            ismatch(r"^\d{4}$", git_dir_names[end]))
          # Current is a document. Use current version, ignoring latest.
          insert!(tree_builder, current_name, current_id, int(filemode(current_entry)))
        else
          current_git_dir_names = tuple(git_dir_names..., current_name)
          latest_entry = latest_tree === nothing ? nothing : latest_tree[current_name]
          merge(repository, current_git_dir_names, lookup_tree(repository, current_id),
            lookup_tree(repository, parent_id),
            latest_entry === nothing ? nothing : lookup_tree(repository, Oid(latest_entry)),
            section_by_git_dir_names, tree_builder_by_git_dir_names)
        end
      end
    end
  end
  if parent_tree !== nothing
    for parent_entry in parent_tree
      parent_name = parent_entry.name
      if current_tree[parent_name] === nothing && tree_builder[parent_name] !== nothing
        delete!(tree_builder, parent_name)
      end
    end
  end
end

function parse_command_line()
  arg_parse_settings = ArgParseSettings()
  @add_arg_table arg_parse_settings begin
    "--verbose", "-v"
      action = :store_true
      help = "increase output verbosity"
    "repository"
      help = "path of Git repository to create"
      required = true
  end
  return parse_args(arg_parse_settings)
end

main()

end # module


text:
#=

  Set covering problem in Picat.

  This example is from the OPL example covering.mod
  """
  Consider selecting workers to build a house. The construction of a
  house can be divided into a number of tasks, each requiring a number of
  skills (e.g., plumbing or masonry). A worker may or may not perform a
  task, depending on skills. In addition, each worker can be hired for a
  cost that also depends on his qualifications. The problem consists of
  selecting a set of workers to perform all the tasks, while minimizing the
  cost. This is known as a set-covering problem. The key idea in modeling
  a set-covering problem as an integer program is to associate a 0/1
  variable with each worker to represent whether the worker is hired.
  To make sure that all the tasks are performed, it is sufficient to
  choose at least one worker by task. This constraint can be expressed by a
  simple linear inequality.
  """

  Solution from the OPL model:
  """
  Optimal solution found with objective: 14
  crew= {23 25 26}
  """

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function covering_opl(problem,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   # "all_solutions"=> all_solutions,
                                                            "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    num_workers = problem[:num_workers]
    qualified   = problem[:qualified]
    cost        = problem[:cost]
    num_tasks = problem[:num_tasks]

    # Which workers to hire
    @variable(model, hire[1:num_workers], Bin)
    @variable(model, 0 <= total_cost <= num_workers*maximum(cost), Int)

    scalar_product(model,cost,hire,total_cost)
    for j in 1:num_tasks
        @constraint(model, sum([hire[w] for w in qualified[j]]) >= 1)
    end

    @objective(model,Min, total_cost)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                hire_val = convert.(Integer,JuMP.value.(hire; result=sol))
                total_cost_val = convert.(Integer,JuMP.value.(total_cost; result=sol))
                println("total_cost:$total_cost_val")
                # println("hire:$hire_val")
                println("hire:", [i for i in 1:num_workers if hire_val[i] == 1 ])
                println()

            end
        end
    else
        println("status:$status")
    end

    return status
end

covering_opl_problems = Dict(
  :1 => Dict(
        :num_workers => 32,

        # The workers qualified for the works
        :qualified =>
        [
            [  1,  9, 19, 22, 25, 28, 31 ],
            [  2, 12, 15, 19, 21, 23, 27, 29, 30, 31, 32 ],
            [  3, 10, 19, 24, 26, 30, 32 ],
            [  4, 21, 25, 28, 32 ],
            [  5, 11, 16, 22, 23, 27, 31 ],
            [  6, 20, 24, 26, 30, 32 ],
            [  7, 12, 17, 25, 30, 31 ] ,
            [  8, 17, 20, 22, 23  ],
            [  9, 13, 14, 26, 29, 30, 31 ],
            [ 10, 21, 25, 31, 32 ],
            [ 14, 15, 18, 23, 24, 27, 30, 32 ],
            [ 18, 19, 22, 24, 26, 29, 31 ],
            [ 11, 20, 25, 28, 30, 32 ],
            [ 16, 19, 23, 31 ],
            [  9, 18, 26, 28, 31, 32 ]
        ],

        :num_tasks => 15,

        # cost per worker
        :cost => [1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3,
                3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9 ]
    )
)

@time covering_opl(covering_opl_problems[:1])

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
module ReferenceFEInterfacesTests

using Test
using InplaceArrays.Fields
using InplaceArrays.Polynomials
using InplaceArrays.ReferenceFEs

D = 2
T = Float64
order = 1
prebasis = MonomialBasis{D}(T,order)

polytope = QUAD
x = get_vertex_coordinates(polytope)
dofs = LagrangianDofBasis(T,x)

facedofids = [[1],[2],[3],[4],Int[],Int[],Int[],Int[],Int[]]

reffe = GenericRefFE(polytope,prebasis,dofs,facedofids)
test_reference_fe(reffe)

shapefuns = get_shapefuns(reffe)

@test evaluate(shapefuns,x) == [1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0]

end # module


text:
@testset "Minimum knapsack problem tests" begin
    @test minimumKnapsackProblem([3, 5, 6], [1, 2, 3], 5) == [0, 1, 1]
end


text:
# First-order reversible series reaction INCONSISTENT NOTATION IN TABLE, CHECK LATER
id = :PerezGalvan2017b

pL = [2.0; 1.0]
pU = [6.0; 3.0]
PerezGalvan2017b_x0(p::Vector{T}) where T = T[one(T); zero(T)]
function PerezGalvan2017b_f!(du::Vector{T}, u::Vector{T}, p::Vector{T}, t) where T
    k1 = p[1]
    k1rev = p[2]
    k2 = 2.0
    k2rev = 20.0
    du[1] = -k1*u[1] + k1rev*u[2]
    du[2] = k1*u[1] - (k1rev + k2)*u[2] + k2rev*(1.0 - u[1] - u[2])
    return
end
tspan = (0.0, 0.5)

prob = ODERelaxProb(PerezGalvan2017b_f!, tspan, PerezGalvan2017b_x0, pL, pU)
url = "https://www.sciencedirect.com/science/article/abs/pii/S0098135417300923"
source = "PerezGalvan2017"
desc = "First-order reversible series reaction"