text:
immutable MPoly{T}
    terms::OrderedDict{Vector{Int},T}
    vars::Vector{Symbol}
end

terms(p::MPoly) = p.terms

eltype{T}(::MPoly{T}) = T

vars(p::MPoly) = p.vars

nvars(p::MPoly) = length(vars(p))

zero{T}(::Type{MPoly{T}}; vars::Vector{Symbol}=Symbol[]) = MPoly{T}(OrderedDict{Vector{Int},T}(), vars)
zero{T}(p::MPoly{T}) = zero(MPoly{T}, vars=vars(p))

one{T}(::Type{MPoly{T}}; vars::Vector{Symbol}=Symbol[]) = MPoly{T}(OrderedDict(zeros(Int, length(vars)) => one(T)), vars)
one{T}(p::MPoly{T}) = one(MPoly{T}, vars=vars(p))

function generators{T}(::Type{MPoly{T}}, vars::Symbol...)
     m = eye(Int, length(vars))
     [MPoly{T}(OrderedDict(m[:,i] => one(T)), [vars...]) for i = 1:length(vars)]
end

generator{T}(::Type{MPoly{T}}, var::Symbol) =
     MPoly{T}(OrderedDict([1] => one(T)), [var])

call{T}(::Type{MPoly{T}}, var::Symbol) =
     MPoly{T}(OrderedDict([1] => one(T)), [var])

call(::Type{MPoly}, var::Symbol) =
     MPoly{Float64}(var)

promote_rule{T,U}(::Type{MPoly{T}}, ::Type{MPoly{U}}) = MPoly{promote_type(T, U)}
promote_rule{T,U}(::Type{MPoly{T}}, ::Type{U}) = MPoly{promote_type(T, U)}

function convert{T}(P::Type{MPoly{T}}, p::MPoly{T})
    p
end

function convert{T}(P::Type{MPoly{T}}, p::MPoly)
    r = zero(P, vars=vars(p))
    for (m, c) in p
        r[m] = convert(T, c)
    end
    r
end

convert{T}(::Type{MPoly{T}}, c::T) =
    MPoly{T}(OrderedDict(Int[] => c))

monomials(p::MPoly) = keys(terms(p))

getindex{T}(p::MPoly{T}, m::Vector{Int}) =
    get(terms(p), m, zero(T))

getindex(p::MPoly, m::Int...) = p[[m...]]

function setindex!{T}(p::MPoly{T}, v, m::Int...)
    v = convert(T, v)
    if method_exists(isapprox, (T,T)) && isapprox(v, zero(T))
        delete!(terms(p), [m...])
    else
        terms(p)[[m...]] = v
    end
end

setindex!{T}(p::MPoly{T}, v, m::Vector{Int}) = p[m...] = v

start(p::MPoly) = start(terms(p))
next(p::MPoly, state) = next(terms(p), state)
done(p::MPoly, state) = done(terms(p), state)

copy{T}(p::MPoly{T}) = MPoly{T}(copy(terms(p)), copy(vars(p)))

function deg(p::MPoly)
    d = 0
    for m in monomials(p)
        s = sum(m)
        if s > d
            d = s
        end
    end
    d
end

function ==(p::MPoly, q::MPoly)
    for (m, c) in p
        if !isapprox(q[m], c)
            return false
        end
    end
    for (m, c) in q
        if !isapprox(p[m], c)
            return false
        end
    end
    true
end


text:
@doc raw"""
```
binarize(Balanced(), img)
```
Binarizes the image using the balanced histogram thresholding method.

# Output

Returns the binarized image as an `Array{Gray{Bool},2}`.

# Details
In balanced histogram thresholding, one interprets a  bin as a  physical weight
with a mass equal to its occupancy count. The balanced histogram method involves
iterating the following three steps: (1) choose the midpoint bin index as a
"pivot",  (2) compute the combined weight to the left and right of the pivot bin
and (3) remove the leftmost bin if the left side is the heaviest, and the
rightmost bin otherwise. The algorithm stops when only a single bin remains. The
last bin determines the sought-after threshold with which the image is
binarized.

Let ``f_n`` (``n = 1 \ldots N``) denote the number of observations in the ``n``th
bin of the image histogram. The balanced histogram method constructs a sequence
of nested intervals

```math
[1,N] \cap \mathbb{Z} \supset I_2 \supset I_3 \supset \ldots \supset I_{N-1},
```
where for ``k = 2 \ldots N-1``
```math
I_k = \begin{cases}
   I_{k-1} \setminus \{\min \left( I_{k-1} \right) \} &\text{if } \sum_{n = \min \left( I_{k-1} \right)}^{I_m}f_n \gt   \sum_{n =  I_m + 1}^{ \max \left( I_{k-1} \right)} f_n, \\
   I_{k-1} \setminus \{\max \left( I_{k-1} \right) \} &\text{otherwise},
\end{cases}
```
and ``I_m = \lfloor \frac{1}{2}\left(  \min \left( I_{k-1} \right) +  \max \left( I_{k-1} \right) \right) \rfloor ``.
The final interval ``I_{N-1}`` consists of a single element which is the bin index
corresponding to the desired threshold.

If one interprets a bin as a physical weight with a mass equal to its occupancy
count, then each step of the algorithm can be conceptualised as removing the
leftmost or rightmost bin to "balance" the resulting histogram on a pivot. The
pivot is defined to be the midpoint between the start and end points of the
interval under consideration.

If it turns out that the single element in ``I_{N-1}`` equals ``1`` or ``N`` then
the original histogram must have a single peak and the algorithm has failed to
find a suitable threshold. In this case the algorithm will fall back to using
the `UnimodalRosin` method to select the threshold.

# Arguments

The function argument is described in more detail below.

##  `img`

An `AbstractArray` representing an image. The image is automatically converted
to `Gray` in order to construct the requisite graylevel histogram.

# Example

Binarize the "cameraman" image in the `TestImages` package.

```julia
using TestImages, ImageBinarization

img = testimage("cameraman")
img_binary = binarize(Balanced(), img)
```

# Reference

1. “BI-LEVEL IMAGE THRESHOLDING - A Fast Method”, Proceedings of the First International Conference on Bio-inspired Systems and Signal Processing, 2008. Available: [10.5220/0001064300700076](https://doi.org/10.5220/0001064300700076)
"""
function binarize(algorithm::Balanced,  img::AbstractArray{T,2}) where T <: Colorant
  img₀₁ = zeros(Gray{Bool}, axes(img))
  edges, counts = build_histogram(img,  256)
  t = find_threshold(HistogramThresholding.Balanced(), counts[1:end], edges)
  for i in CartesianIndices(img)
    img₀₁[i] = img[i] < t ? 0 : 1
  end
  img₀₁
end

ext: jl


text:
# Expected Transformation
#
# [Lib.]dnnl_f(a,b,c) -> Lib.dnnl_f(dnnl_convert(a), dnnl_convert(b), dnnl_convert(c))
#
# Pretty straightforward.
# Note that `dnnl_convert` should return valid GC tracked objects and not raw pointers
# because the results of `dnnl_convert` are passed through the
# `Base.cconvert` -> `Base.unsafe_convert` chain.
macro apicall(expr)
    # Split into a partial implementation for testing purposes.
    expr = _apicall_partial_impl(expr)
    return quote
        status = $expr
        if status != Lib.dnnl_success
            error("DNNL Failure: $status")
        end
        status
    end
end

function _apicall_partial_impl(expr)
    if expr.head != :call
        error("Only call `@apicall` on function calls")
    end

    # Prefix "Lib." in front of the function call.
    # However, sometimes the function to call is passed as a higher order function.
    # Thus, we only implicitly attach "Lib" is the function name starts with "dnnl".
    fname = expr.args[1]
    if isa(fname, Symbol)
        if startswith(string(fname), "dnnl")
            fname = :(Lib.$fname)
        else
            fname = :($(esc(fname)))
        end
    end

    # Escape and convert each of the arguments.
    args = expr.args[2:end]
    for i in eachindex(args)
        # Handle splats.
        arg = args[i]
        if isa(arg, Expr) && arg.head == :...
            args[i] = :(_dnnl_convert($(esc(arg.args[1]))...)...)
        else
            args[i] = :(dnnl_convert($(esc(args[i]))))
        end
    end
    return :($fname($(args...)))
end

# Get to the ultimate parent.
# Strategy: keep calling `parent` until `typeof(parent(x)) == typeof(x)`.
ancestor(x::AbstractArray) = ancestor(x, parent(x))
ancestor(x::AbstractArray, y::AbstractArray) = ancestor(y, parent(y))
ancestor(x::T, ::T) where {T<:AbstractArray} = x

#####
##### Runtime Arguments
#####

# In the macro below, add support for converting expressions like `:(a.b)` to just
# plain `b`
getsym(x::Symbol) = string(x)
getsym(x::QuoteNode) = getsym(x.value)
getsym(x::Expr) = getsym(last(x.args))

# Note - these definitions aren't strictly necessary, but useful for CachedArrays at
# the moment so I'm keeping them.
abstract type AccessContext end
struct Reading <: AccessContext end
struct Writing <: AccessContext end
struct Unknown <: AccessContext end

function dnnl_arg(x, y, context::AccessContext = Reading())
    return Lib.dnnl_exec_arg_t(x, dnnl_exec_arg(y, context))
end

Base.length(::Lib.dnnl_exec_arg_t) = 1
Base.iterate(x::Lib.dnnl_exec_arg_t, s = (x, nothing)) = s

function dnnl_arg(x, y::Union{AbstractArray{<:AbstractArray}, Tuple}, context::AccessContext = Reading())
    return map(enumerate(y)) do (i, _y)
        Lib.dnnl_exec_arg_t(x + (i - 1), dnnl_exec_arg(_y, context))
    end
end

function dnnl_exec_arg(y::T, context) where {T}
    return error("Define `dnnl_exec_arg` for type $(T)!")
end

# In general, we want there to have a static length.
# However, for primitives like Concat that may have many args, using a Vector makes more
# sense.
const TupleOrVector{T} = Union{NTuple{N,T},AbstractVector{T}} where {N}
struct Arguments{T<:TupleOrVector{Lib.dnnl_exec_arg_t}}
    args::T
end

# construct `Arguments` directly.
Arguments(args::Lib.dnnl_exec_arg_t...) = Arguments(args)
Base.length(a::Arguments) = length(a.args)

# Note: `Base.cconvert` should return valid Julia objects and not pointers.
# We only really use `cconvert` for the `Arguments` based types.
# Everything else needs to go through `Base.unsafe_convert` for the GC protection.
function Base.cconvert(::Type{Ptr{Lib.dnnl_exec_arg_t}}, x::Arguments{<:NTuple})
    return Ref(x.args)
end

function Base.cconvert(::Type{Ptr{Lib.dnnl_exec_arg_t}}, x::Arguments{<:AbstractVector})
    return x.args
end

Base.resize!(x::Arguments{<:AbstractVector}, i) = resize!(x.args, i)
Base.setindex!(x::Arguments{<:AbstractVector}, v, i::Integer) = setindex!(x.args, v, i)
Base.getindex(x::Arguments, i) = getindex(x.args, i)
Base.lastindex(x::Arguments) = lastindex(x.args)

append(x::Arguments{<:Tuple}, arg::Lib.dnnl_exec_arg_t) = Arguments(x.args..., arg)
append(x::Arguments{<:Tuple}, y::Arguments{<:Tuple}) = Arguments(x.args..., y.args...)
append(x::Arguments{<:AbstractVector}, arg::Lib.dnnl_exec_arg_t) = push!(x.args, arg)
append(x::Tuple{<:AbstractArray,<:Arguments}, y) = (x[1], append(x[2], y[2]))
#####
##### Construction Utilities
#####

# General transformation provided by this macro:
#
# `@dnnl_args src dst`
#
# becomes
#
# ```
# (OneDNN).make_args(
#   (OneDNN).dnnl_arg(Lib.DNNL_ARG_SRC, src, (OneDNN).Reading()),
#   (OneDNN).dnnl_arg(Lib.DNNL_ARG_DST, dst, (OneDNN).Writing()),
# )
# ```
#
# The argument symbols, (`src` and `dst` in the above example) get converted to
# `Lib.DNNL_ARG_SRC` and `Lib.DNNL_ARG_DST` respectively.
#
# We also try to automatically determine whether an object is being written or read
# based on the name of the argument.
#
# As far as I can tell, this should be pretty reliable, but we'll see.
#
# This requires that functions using this macro follow the same naming conventions as the
# OneDNN C-API, but that's probably good practice anyways.
const _R = :(Reading())
const _W = :(Writing())
const CONTEXT_MAP = [
    "DIFF_SRC" => _W,
    "DIFF_WEIGHTS" => _W,
    "DIFF_BIAS" => _W,
    "DIFF_DST" => _R,
    "SRC" => _R,
    "MULTIPLE_SRC" => _R,
    "WEIGHTS" => _R,
    "BIAS" => _R,
    "DST" => _W,
    "FROM" => _R,
    "TO" => _W,
    "SCALE_SHIFT" => _R,
    "DIFF_SFALE_SHIFT" => _W,
    "MEAN" => _R,
    "VARIANCE" => _R,
    "WORKSPACE" => _W,
]

macro dnnl_args(syms...)
    extra_arg = Any[]
    exprs = map(syms) do sym
        # Handle direct interpolation.
        if isa(sym, Expr)
            if sym.head == :$
                newarg = gensym()
                @assert length(sym.args) == 1
                return :($(esc(sym.args[1])))
            else
                error("Unknown expression: $sym")
            end
        end

        # Handle scratchpad generation.
        dnnl_arg_enum = "DNNL_ARG_$(uppercase(getsym(sym)))"
        # Determing the context
        context = nothing
        for (str, quot) in CONTEXT_MAP
            if occursin(str, dnnl_arg_enum)
                context = quot
                break
            end
        end

        if context === nothing
            error("Unknown DNNL Argument: $dnnl_arg_enum")
        end
        if length(extra_arg) == 0
            push!(extra_arg, esc(sym))
        end

        return :(dnnl_arg(Lib.$(Symbol(dnnl_arg_enum)), $(esc(sym)), $context))
    end

    # Capture one argument to forward into `execute!` if we're using `similar` for
    # scratchpad allocation.
    if SIMILAR_FOR_SCRATCHPAD
        return :((ancestor($(extra_arg[1])), make_args($(exprs...))))
    else
        return :(make_args($(exprs...)))
    end
end

# Generic path, if any of the args is a vector, then result will also be a vector.
# Otherwise, we have tuple arguments.
function make_args(args...)
    return Arguments(_hasarray(args...) ? _vcat(args) : _flatcat(args...))
end

_hasarray() = false
_hasarray(x, y...) = _hasarray(y...)
_hasarray(x::AbstractArray, y...) = true

_vcat(args; init = Lib.dnnl_exec_arg_t[]) = reduce(append!, args; init)

_flatcat(x, y...) = (x, _flatcat(y...)...)
_flatcat(x::Tuple, y...) = (_flatcat(x...)..., _flatcat(y...)...)
_flatcat() = ()

#####
##### Conversion Hooks
#####

dnnl_convert(x) = x
_dnnl_convert(x, y...) = (dnnl_convert(x), _dnnl_convert(y...)...)
_dnnl_convert(x) = (dnnl_convert(x),)

const MaybeRef{T} = Union{T,Ref{<:T}}
const MaybePtr{T} = Union{T,Ptr{<:T}}
wrap_ref(x::Ref) = x
wrap_ref(x) = Ref(x)
unwrap_ref(x::Ref) = x[]
unwrap_ref(x) = x

#####
##### Dimension Helpers
#####

expand(N, i::Tuple) = i
expand(N, i::Integer) = ntuple(_ -> i, N)

_paddims(x::Tuple, y::Tuple) = (x..., y[(end - (length(y) - length(x) - 1)):end]...)

struct Dims{N}
    kernel::NTuple{N,Int}
    strides::NTuple{N,Int}
    dilation::NTuple{N,Int}
    padding::NTuple{N,Int}
end

function _output_size(dims::Dims, sz::Int, i::Int)
    @unpack kernel, strides, dilation, padding = dims
    numerator = sz + 2 * padding[i] - ((kernel[i] - 1) * (dilation[i] + 1) + 1)
    return div(numerator, strides[i]) + 1
end

function output_size(sz::NTuple{N,Int}, dims::Dims{M}) where {N,M}
    return _paddims(ntuple(i -> _output_size(dims, sz[i], i), Val(M)), sz)
end

function tuple_replace(x::NTuple{N,Int}, v, ::Val{I}) where {N,I}
    return ntuple(i -> (i == I) ? v : x[i], Val(N))
end

#####
##### Kind Types
#####

# Move the primitive kind enums into the type domain to give the Julia compiler better
# type-level information when the return type depends on things like inference vs traiing.

abstract type AbstractKind end
abstract type AbstractForwardKind <: AbstractKind end
abstract type AbstractBackwardKind <: AbstractKind end

struct Inference <: AbstractForwardKind end
struct Training <: AbstractForwardKind end
struct Backward <: AbstractBackwardKind end

kind(::Inference) = Lib.dnnl_forward_inference
kind(::Training) = Lib.dnnl_forward_training
kind(::Backward) = Lib.dnnl_backward

dnnl_convert(x::AbstractKind) = kind(x)


text:
vnp_year = 2020
vnp_fname = "/home/psomkuti/VNP43C1_2020_full_compressed.h5"

function create_VNP_SD_from_locbounds(lon_min::Real, lat_min::Real,
                                      lon_max::Real, lat_max::Real;
                                      vnp_fname::String=vnp_fname)

    vnp_h5 = h5open(vnp_fname, "r")

    vnp_lons = collect(-180:0.05:180)
    vnp_lats = collect(-60:0.05:90)
    vnp_variables = vnp_h5["names"][:]
    vnp_fill_values = vnp_h5["fill_values"][:]
    vnp_year = 2020
    vnp_times = DateTime(vnp_year, 1, 1) .+ (p -> Dates.Day(min(p, 366))).(0:size(vnp_h5["VNP43C1"])[3] - 1)

    # Trim down to our Antarctica limit
    if lat_min < -60
        lat_min = -59.99999
    end

    idx_lon_min = searchsortedfirst(vnp_lons, lon_min) - 1
    idx_lon_max = searchsortedfirst(vnp_lons, lon_max) - 1

    idx_lat_min = 3001 - searchsortedfirst(vnp_lats, lat_max) + 1
    idx_lat_max = 3001 - searchsortedfirst(vnp_lats, lat_min) + 1

    #println(idx_lon_min, ": ", vnp_lons[idx_lon_min])
    #println(idx_lon_max, ": ", vnp_lons[idx_lon_max])
    #println(idx_lat_min, ": ", vnp_lats[end:-1:1][idx_lat_max])
    #println(idx_lat_max, ": ", vnp_lats[end:-1:1][idx_lat_min])

    # Load data and flip along lat dimension (original data is in decreasing lats)
    data = vnp_h5["VNP43C1"][idx_lat_min:idx_lat_max, idx_lon_min:idx_lon_max,:,:][end:-1:1,:,:,:]

    close(vnp_h5)

    return VNPData(
        data,
        vnp_lons[idx_lon_min:idx_lon_max+1],
        vnp_lats[end:-1:1][idx_lat_min:idx_lat_max+1][end:-1:1],
        vnp_times,
        vnp_year,
        vnp_variables,
        vnp_fill_values
    )

end

function sample_vnp_data(loctime::GeolocationTime, var_list, vnp_sd::VNPData)

    if !(typeof(var_list) <: Vector)
         var_list = [var_list]
    end

    for var in var_list
        if !(var in vnp_sd.names)
            @error "Sorry! Variable \"$(var)\" not preseng in VNPData object."
            return nothing
        end
        @debug "Variable check OK - all supplied names are found in VNP SD."
    end

    if ((loctime.loc.lon < vnp_sd.lon_bounds[1]) |
        (loctime.loc.lon > vnp_sd.lon_bounds[end]))
        @error "Sorry! Location $(loctime) is outside of longitude bounds!"
        return nothing
    end

    if ((loctime.loc.lat < vnp_sd.lat_bounds[1]) |
        (loctime.loc.lat > vnp_sd.lat_bounds[end]))
        @error "Sorry! Location $(loctime) is outside of latitude bounds!"
        return nothing
    end

    # Find spatial indices
    idx_x = searchsortedfirst(vnp_sd.lon_bounds, loctime.loc.lon) - 1
    idx_y = searchsortedfirst(vnp_sd.lat_bounds, loctime.loc.lat) - 1

    #idx_xplus = min(idx_x + 1, length(SD.lon_bounds))
    #idx_xminus = max(idx_x - 1, 1)
    #x_fac = (loctime.loc.lon - SD.lon_bounds[idx_x]) / (SD.lon_bounds[idx_xplus] - SD.lon_bounds[idx_x])

    #idx_y = searchsortedfirst(SD.lat_bounds, loctime.loc.lat) - 1
    #idx_yplus = min(idx_y + 1, length(SD.lat_bounds))
    #idx_yminus = max(idx_y - 1, 1)
    #y_fac = (loctime.loc.lat - SD.lat_bounds[idx_y]) / (SD.lat_bounds[idx_yplus] - SD.lat_bounds[idx_y])

    #center_weight = 0.25
    #w = (1.0 - center_weight) / 2.0
    # Spatial weights using x_fac and y_fac
    #sp_w = [0       (1-y_fac) * w      0;
    #        (1-x_fac) * w    center_weight      (x_fac) * w;
    #        0        (y_fac) * w       0]

    # Create new time based on the year of the surface dataset
    # (just replacing the year according to the dataset year)
    newtime = DateTime(
        vnp_sd.data_year,
        Dates.month(loctime.time),
        Dates.day(loctime.time),
        Dates.hour(loctime.time),
        Dates.minute(loctime.time),
        Dates.second(loctime.time)
    )

    idx_t = searchsortedfirst(vnp_sd.time_bounds, newtime) - 1

    idx_var = searchsortedfirst.(Ref(vnp_sd.names), var_list)

    return vnp_sd.data[idx_y, idx_x, idx_t, idx_var]

    #idx_tplus = min(idx_t + 1, length(SD.time_bounds) - 1)
    #t_fac = (newtime - SD.time_bounds[idx_t]) / (SD.time_bounds[idx_tplus] - SD.time_bounds[idx_t])

    #println(idx_tplus, " / ", length(SD.time_bounds))
    #println(loctime)
    #println(idx_t, ", ", idx_tplus, ", factor: ", t_fac)
    #println(idx_xminus, ", ", idx_x, ", ", idx_xplus, ", factor: ", x_fac)
    #println(idx_yminus, ", ", idx_y, ", ", idx_yplus, ", factor: ", y_fac)

    #sampled = SD.data[idx_t][idx_xminus:idx_xplus, idx_yminus:idx_yplus]

    #output = tuple(
    #    [round(sum(((p -> p[i]).(sampled)) .* sp_w)) for i in 1:SD.length]...
    #)

end

function calculate_reflectance(loctime::GeolocationTime, sza::Real, band::String, vnp_sd::VNPData)

    # Stitch together variable names
    var_list = ["BRDF_Albedo_Parameter$(i)_$(band)" for i in 1:3]

    g0_iso = 1.0
    g1_iso = 0.0
    g2_iso = 0.0

    g0_vol = -0.007574
    g1_vol = -0.070987
    g2_vol =  0.307588

    g0_geo = -1.284909
    g1_geo = -0.166314
    g2_geo =  0.041840

    f_iso, f_vol, f_geo = sample_vnp_data(loctime, var_list, vnp_sd)

    # Replace fill values by NaNs

    if f_iso == 32767
        f_iso = NaN
    end

    if f_geo == 32767
        f_geo = NaN
    end

    if f_vol == 32767
        f_vol = NaN
    end

    # Compute black-sky albedo
    sza2 = deg2rad(sza)
    sza3 = deg2rad(sza)

    # 1 / 1000.0 factor is needed to convert Int-valued VIIRS data into
    # weights.

    # White sky albedo
    #return (f_iso * 1.0 + f_vol * 0.189184 + f_geo * (-1.377622)) * 0.001

    # Black sky albedo
    return f_iso * 0.001 *  (g0_iso + g1_iso * sza2 + g2_iso * sza3) +
        f_vol * 0.001 * (g0_vol + g1_vol * sza2 + g2_vol * sza3) +
        f_geo * 0.001 * (g0_geo + g1_geo * sza2 + g2_geo * sza3)
end

function calculate_NDVI(loctime::GeolocationTime, sza::Real)

    R_nir, R_vis = calculate_reflectance(loctime, sza, ["M7", "M5"])

    if R_nir == 0.0
        return 0.0
    elseif R_vis == 0.0
        return 0.0
    else
        return (R_nir - R_vis) / (R_nir + R_vis)
    end

end

function calculate_NIRv(loctime::GeolocationTime, sza::Real)

    R_nir = calculate_reflectance(loctime, sza, "NIR")
    NDVI = calculate_NDVI(loctime, sza)

    return R_nir * NDVI
end


text:
using MieSeries

a = 1.0
ϵ = 1.0
μ = 1.0
sphere = PECSphere(a, ϵ, μ)

ω   = 2π
k   = ω * √(ϵ * μ)
pol = :x
exc = excite(sphere, ω, pol)

Θ = linspace(0.0,   π, 10)
Φ = linspace(0.0, 0.0,  1)
P = [ [cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)] for θ in Θ for ϕ in Φ]

N = length(P)
E = zeros(N)
H = zeros(N)

for (i,p) in enumerate(P)
    Ex, Ey, Ez, Hx, Hy, Hz = farfield(sphere, exc, p)
    E[i] = norm([Ex,Ey,Ez])
    H[i] = norm([Ex,Ey,Ez])
end

# Define the grid on which to compute the scattered field
x = linspace(-2.0, 2.0, 3)
z = linspace(-4.0, 4.0, 400)
P = [ [ξ, 0, ζ] for ζ in z for ξ in x ]

# Compute the scattered field
Ex = zeros(Complex128, length(P))
Ey = zeros(Complex128, length(P))
Ez = zeros(Complex128, length(P))
for (i,p) in enumerate(P)
    norm(p) < a && continue
    ex, ey, ez, _ = scatteredfield(sphere, exc, p)
    Ex[i] = ex + exp(im*k*p[3])
    Ey[i] = ey
    Ez[i] = ex
end
Ex = reshape(Ex, length(x), length(z))
Ey = reshape(Ey, length(x), length(z))
Ez = reshape(Ez, length(x), length(z))

# Plot the field
A = real(sqrt(Ex.*conj(Ex)))

# using Plots
# n = div(length(x),2) + 1
# plot(A[n,:], m=:point)
# savefig("mie.svg")


text:
# core functionalities
# the suffix ♂(\mars) means this value is 0-based index.

# grid graph
GCoptimizationGridGraph(width::T, height::T, numLabels::Integer) where {T <: Integer} = GCoptimizationGridGraph(GCOSiteID(width), GCOSiteID(height), GCOLabelID(numLabels))
setSmoothCostVH(gco::Cxx.CppPtr, smoothArray::Vector{T}, vCosts::Vector{T},
                hCosts::Vector{T}) where {T <: Integer} = setSmoothCostVH(gco, pointer(convert(Vector{GCOEnergyTermType}, smoothArray)),
                                                          pointer(convert(Vector{GCOEnergyTermType}, vCosts)),
                                                          pointer(convert(Vector{GCOEnergyTermType}, hCosts)))

# general graph
GCoptimizationGeneralGraph(numSites::Integer, numLabels::Integer) = GCoptimizationGeneralGraph(GCOSiteID(numSites), GCOLabelID(numLabels))
setNeighbors(gco::Cxx.CppPtr, site1♂::T, site2♂::T) where {T <: Integer} = setNeighbors(gco, GCOSiteID(site1♂), GCOSiteID(site2♂))
setNeighbors(gco::Cxx.CppPtr, site1♂::T, site2♂::T, weight::Integer) where {T <: Integer} = setNeighbors(gco, GCOSiteID(site1♂), GCOSiteID(site2♂), GCOEnergyTermType(weight))

# expansion
expansion(gco::Cxx.CppPtr, maxIterationNum::Integer) = expansion(gco, Cint(maxIterationNum))
alpha_expansion(gco::Cxx.CppPtr, alphaLabel♂::Integer) = alpha_expansion(gco, GCOLabelID(alphaLabel♂))

# swap
swap(gco::Cxx.CppPtr, maxIterationNum::Integer) = swap(gco, Cint(maxIterationNum))
alpha_beta_swap(gco::Cxx.CppPtr, alphaLabel♂::T, betaLabel♂::T) where {T <: Integer} = alpha_beta_swap(gco, GCOLabelID(alphaLabel♂), GCOLabelID(betaLabel♂))
alpha_beta_swap(gco::Cxx.CppPtr, alphaLabel♂::T, betaLabel♂::T,
                alphaSites♂::Vector{N}, alphaSize::N, betaSites♂::Vector{N},
                betaSize::N) where {T <: Integer,N <: Integer} = alpha_beta_swap(gco, GCOLabelID(alphaLabel♂), GCOLabelID(betaLabel♂), pointer(convert(Vector{GCOSiteID}, alphaSites♂)),
                                                     GCOSiteID(alphaSize), pointer(convert(Vector{GCOSiteID}, betaSites♂), GCOSiteID(betaSize)))

# setDataCost
setDataCost(gco::Cxx.CppPtr, dataArray::Vector{T}) where {T <: Integer} = setDataCost(gco, pointer(convert(Vector{GCOEnergyTermType}, dataArray)))
setDataCost(gco::Cxx.CppPtr, site♂::Integer, label♂::Integer, energy::Integer) = setDataCost(gco, GCOSiteID(site♂), GCOLabelID(label♂), GCOEnergyTermType(energy))
setDataCost(gco::Cxx.CppPtr, label♂::Integer, costs::GCOSparseDataCost, count::Integer) = setDataCost(gco, GCOLabelID(label♂), pointer(costs), GCOSiteID(count))

# setSmoothCost
setSmoothCost(gco::Cxx.CppPtr, label1♂::T, label2♂::T, energy::Integer) where {T <: Integer} = setSmoothCost(gco, GCOLabelID(label1♂), GCOLabelID(label2♂), GCOEnergyTermType(energy))
setSmoothCost(gco::Cxx.CppPtr, smoothArray::Vector{T}) where {T <: Integer} = setSmoothCost(gco, pointer(convert(Vector{GCOEnergyTermType}, smoothArray)))

# setLabelCost
setLabelCost(gco::Cxx.CppPtr, cost::Integer) = setLabelCost(gco, GCOEnergyTermType(cost))
setLabelCost(gco::Cxx.CppPtr, costArray::Vector{T}) where {T <: Integer} = setLabelCost(gco, pointer(convert(Vector{GCOEnergyTermType}, costArray)))
setLabelSubsetCost(gco::Cxx.CppPtr, labels♂::Vector{T}, numLabels::Integer, cost::Integer) where {T <: Integer} = setLabelSubsetCost(gco, pointer(convert(Vector{GCOLabelID}, labels♂)), GCOLabelID(numLabels), GCOEnergyTermType(cost))

# whatLabel
whatLabel(gco::Cxx.CppPtr, site♂::Integer) = whatLabel(gco, GCOSiteID(site♂))
whatLabel(gco::Cxx.CppPtr, start♂::T, count::T, labeling::Vector{N}) where {T <: Integer,N <: Integer} = whatLabel(gco, GCOSiteID(start♂), GCOSiteID(count), pointer(convert(Vector{GCOLabelID}, labeling)))

# setLabel
setLabel(gco::Cxx.CppPtr, site♂::Integer, label♂::Integer) = setLabel(gco, GCOSiteID(site♂), GCOLabelID(label♂))

# setLabelOrder
setLabelOrder(gco::Cxx.CppPtr, isRandom::Bool) = setLabelOrder(gco, GCObool(isRandom))
setLabelOrder(gco::Cxx.CppPtr, order::Vector{T}, size::Integer) where {T <: Integer} = setLabelOrder(gco, pointer(convert(Vector{GCOLabelID}, order)), GCOLabelID(size))

# setVerbosity
setVerbosity(gco::Cxx.CppPtr, level::Integer) = setVerbosity(gco, Cint(level))


text:
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# robust_uncertainty.jl
#
# Computes the Value at Risk for a data-driven uncertainty set; see
# "Data-Driven Robust Optimization" (Bertsimas 2013), section 6.1 for
# details. Closed-form expressions for the optimal value are available.
#############################################################################

using JuMP

R = 1
d = 3
𝛿 = 0.05
ɛ = 0.05
N = ceil((2+2log(2/𝛿))^2) + 1

Γ1(𝛿,N) = (R/sqrt(N))*(2+sqrt(2*log(1/𝛿)))
Γ2(𝛿,N) = (2R^2/sqrt(N))*(2+sqrt(2*log(2/𝛿)))

μhat = rand(d)
M = rand(d,d)
Σhat = 1/(d-1)*(M-ones(d)*μhat')'*(M-ones(d)*μhat')

m = Model()

@variable(m, Σ[1:d,1:d], SDP)
@variable(m, u[1:d])
@variable(m, μ[1:d])

@constraint(m, norm(μ-μhat) <= Γ1(𝛿/2,N))
@constraint(m, vecnorm(Σ-Σhat) <= Γ2(𝛿/2,N))

A = [(1-ɛ)/ɛ (u-μ)';
     (u-μ)     Σ   ]
@SDconstraint(m, A >= 0)

c = randn(d)
@objective(m, Max, dot(c,u))

solve(m)

object = getobjectivevalue(m)
exact = dot(μhat,c) + Γ1(𝛿/2,N)*norm(c) + sqrt((1-ɛ)/ɛ)*sqrt(dot(c,(Σhat+Γ2(𝛿/2,N)*eye(d,d))*c))

println("objective value:  $(object)")
println("error from exact: $(abs(exact-object))")


text:
"""
	module BoundaryLayers
"""
module BoundaryLayers

export PML
export cPML
export noBL
export AbstractBL

using ...Defaults
using ..Points
using Formatting

abstract type AbstractBL{SIDE} end
abstract type AbstractRealBL{SIDE} <: AbstractBL{SIDE} end
abstract type AbstractComplexBL{SIDE} <: AbstractBL{SIDE} end

for (TBL,ATBL) ∈ zip((:PML,:cPML,:noBL),
					 (:AbstractComplexBL,:AbstractComplexBL,:AbstractRealBL))
	@eval begin
		mutable struct $TBL{SIDE} <: $ATBL{SIDE}
			depth::Float64
		    start::Float64
		    stop::Float64

			function $TBL{SIDE}(depth,start,stop) where SIDE
				if !($TBL <: noBL)
					depth>0 || throw(ArgumentError("cannot have $($TBL) with 0 depth"))
				end
				return new{SIDE}(depth,start,stop)
			end
		end

		$TBL{SIDE}(depth) where SIDE = $TBL{SIDE}(depth,NaN,NaN)

		function Base.show(io::IO,tbl::$TBL{SIDE}) where SIDE
			print(io,$TBL," (SIDE=",SIDE, ", depth=",fmt("1.1f",tbl.depth),")")
		end

		(tbl::$TBL)(x::Real,y...) = tbl(Point(x,y...))

		(bl::$TBL{SIDE})(;depth=bl.depth, start=bl.start, stop=bl.stop) where SIDE = $TBL{SIDE}(depth,start,stop)
	end
end
noBL{SIDE}() where SIDE = noBL{SIDE}(0,NaN,NaN)

(bl::PML{SIDE})(p::Point) where {SIDE} = conductivity_profile(p[ceil(Int,SIDE/2)],bl.start,bl.stop)
(bl::cPML{SIDE})(p::Point) where {SIDE} = conj_conductivity_profile(p[ceil(Int,SIDE/2)],bl.start,bl.stop)
(bl::noBL)(p::Point) = 0.0

"""
	struct PML{SIDE}(depth) -> pml

absorbing PML along dimension `SIDE` with `depth`. The location of the PML is later determined
by associating it with a shape in a call to `Boundary`.

see also: [`Boundary`](@ref), [`cPML`](@ref), [`noBL`](@ref)

---

	(::PML)(x,y...) -> σ

evaluate the PML conductivity `σ` at the point `x,y...` in local frame.
"""
PML

"""
	struct cPML{SIDE}(depth) -> pml

amplifying conjugate PML along dimension `SIDE` with `depth`. The location of the cPML is later determined
by associating it with a shape in a call to `Boundary`.

see also: [`Boundary`](@ref), [`PML`](@ref), [`noBL`](@ref)

---

	(::cPML)(x,y...) -> σ

evaluate the PML conductivity `σ` at the point `x,y...` in local frame
"""
cPML

"""
	struct noBL{SIDE}(depth) -> nobl

lack of boundary layer along dimension `SIDE`. `depth` is not used, but it there for consistency.
The location is later determined by associating it with a shape in a call to `Boundary`.

see also: [`Boundary`](@ref), [`PML`](@ref), [`cPML`](@ref)

---

	(::noBL)(x,y...) -> σ

returns `σ` = 0. This is here for consistency.
"""
noBL

Base.conj(pml::PML{SIDE}) where SIDE = cPML{SIDE}(pml.depth,pml.start,pml.stop)
Base.conj(cpml::cPML{SIDE}) where SIDE = PML{SIDE}(cpml.depth,cpml.start,cpml.stop)
Base.conj(nbl::noBL) = nbl
"""
	conj(pml) -> cpml
	conj(cpml) -> pml
	conj(nobl) -> nobl

conjugate a boundary layer
"""
conj

# used inside PML. This is the place to modify the PML profile
function conductivity_profile(z::Real,start::Real,stop::Real)
	if (start≤stop && start≤z) || (start≥stop && start≥z)
		depth = abs(stop-start)
	    β = 2*sqrt(depth)
	    α = -(1/4)*exp(complex(0,SCALING_ANGLE))*log(EXTINCTION)/depth/((2exp(β)-β)/(2β*expm1(β))-1/β^2)
        u = abs(z-start)/depth
        return α*u*expm1(β*u)/expm1(β)
    else
        return complex(0.0)
    end
end
function integrated_conductivity_profile(z::Real,start::Real,stop::Real)
	if (start≤stop && start≤z) || (start≥stop && start≥z)
		depth = abs(stop-start)
		β = 2*sqrt(depth)
		α = -(1/4)*exp(complex(0,SCALING_ANGLE))*log(EXTINCTION)/depth/((2exp(β)-β)/(2β*expm1(β))-1/β^2)
		u = abs(z-start)/depth
		return (α*depth/expm1(β))*(u*exp(β*u)-(1+1/β)*expm1(β*u)/β)
	else
		return complex(0.0)
	end
end

# used inside cPML
conj_conductivity_profile(z,start,stop) = -conj(conductivity_profile(z,start,stop))
conj_integrated_conductivity_profile(z,start,stop) = -conj(integrated_conductivity_profile(z,start,stop))

end # module


text:
module DiffOpt

using Random
using LinearAlgebra
using SparseArrays
using IterativeSolvers: lsqr
using JuMP

import LazyArrays
import BlockDiagonals

using MathOptInterface
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities

using MathOptSetDistances
const MOSD = MathOptSetDistances

using MatrixOptInterface
const MatOI = MatrixOptInterface

const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex

const SUPPORTED_OBJECTIVES = Union{
    MOI.SingleVariable,
    MOI.ScalarAffineFunction{Float64},
    MOI.ScalarQuadraticFunction{Float64}
}

const SUPPORTED_SCALAR_SETS = Union{
    MOI.GreaterThan{Float64},
    MOI.LessThan{Float64},
    MOI.EqualTo{Float64},
    MOI.Interval{Float64}
}

const SUPPORTED_SCALAR_FUNCTIONS = Union{
    MOI.SingleVariable,
    MOI.ScalarAffineFunction{Float64}
}

const SUPPORTED_VECTOR_FUNCTIONS = Union{
    MOI.VectorOfVariables,
    MOI.VectorAffineFunction{Float64},
}

const SUPPORTED_VECTOR_SETS = Union{
    MOI.Zeros,
    MOI.Nonpositives,
    MOI.Nonnegatives,
    MOI.SecondOrderCone,
    MOI.PositiveSemidefiniteConeTriangle,
}

include("utils.jl")
include("conic_diff.jl")
include("quadratic_diff.jl")
include("diff_opt.jl")
include("moi_wrapper.jl")
include("jump_moi_overloads.jl")

export diff_optimizer

end # module


text:
@testset "Pointer" begin
    a = [1,2,3,4]
    pt = Pointer(a)
    @test pt[3] == 3
    pt = Pointer(a, 2)
    @test pt[0] == 2
    @test pt[1] == 3
    pt = pt << 1
    @test pt[0] == 3
    pt[0] = 5
    @test pt[0] == 5
    try
        b = Slice(a, 3:2)
        @test 1 == 0
    catch
        @test 1 == 1
    end
end


text:
function __init__()
    @require TaylorModels = "314ce334-5f6e-57ae-acf6-00b6e903104a" include("init_TaylorModels.jl")
    @require IntervalMatrices = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9" include("init_IntervalMatrices.jl")
end


text:
abstract type Feature end

"""
"""
series(::Dataset, ::Feature; kwargs...) = nothing


text:
module Utils

    export rel_err, l2_err

    include("utils/diag_utils.jl")

    export @define, @reexport, @dec128

    include("utils/macro_utils.jl")

    export compensated_summation, L2norm, l2norm

    include("utils/sum_utils.jl")

    export istriustrict, istrilstrict
    export simd_copy_xy_first!, simd_copy_yx_first!, simd_copy_yx_second!,
           simd_copy_yx_first_last!,
           simd_axpy!, simd_aXbpy!, simd_abXpy!, simd_mult!

    include("utils/matrix_utils.jl")

end


text:
module ChaosQuantified

export Rqa, DelayEmbedding, entropy, conditional_entropy

include("helpers.jl")
include("entropy.jl")
include("Rqa.jl")

end # module


text:
using Revise
using Caesar, RoME, IncrementalInference, DistributedFactorGraphs
using DelimitedFiles
using Test

# Visualization
using Colors, Gadfly
using RoMEPlotting, KernelDensityEstimatePlotting

#posFile = joinpath(dirname(pathof(Caesar)),"..","test","testdata","test_array_positions.csv");
posFile = joinpath(ENV["HOME"],"data","sas","test_array_positions.csv");
posData = readdlm(posFile,',',Float64,'\n')

posData[:,1] .+= 13.0
posData[:,2] .-= 40.0

# gps covariance
rtkCov = Matrix(Diagonal([0.1;0.1].^2))

## build the factor graph
N = 100
fg = initfg()

#landmark
addVariable!(fg, :l1, Point2)
# setVal!(fg, :l1, zeros(2,N))

global pp = nothing
for i in 1:5
    sym = Symbol("x$i")
    addVariable!(fg, sym, Point2)
    setVal!(fg, sym, zeros(2,N))
    pp = PriorPoint2(MvNormal(posData[i,1:2], rtkCov) )
    addFactor!(fg, [sym;], pp)
    # manual init for now
    #setValKDE!(getVariable(fg,sym), kde!(getSample(pp, N)[1]))
    #getData(fg,sym).initialized = true
end

#dataFile = joinpath(dirname(pathof(Caesar)),"..","test","testdata","test_array_waveforms.csv");
dataFile = joinpath(ENV["HOME"],"data","sas","test_array_waveforms.csv");

csvWaveData = readdlm(dataFile,',',Float64,'\n')
csvWaveData = csvWaveData[:,1:5]
#
cfgFile = joinpath(ENV["HOME"],"data","sas","SAS2D.yaml");
cfgd=loadConfigFile(cfgFile)
chirpFile = joinpath(ENV["HOME"],"data","sas","chirp250.txt");

sas2d = prepareSAS2DFactor(5, csvWaveData, rangemodel=:Correlator, cfgd=cfgd, chirpFile=chirpFile)

f = addFactor!(fg, [:l1;:x1;:x2;:x3;:x4;:x5], sas2d, autoinit=false)
#, threadmodel=SingleThreaded )

import Base.convert
using JSON2

################################
#### Saving and Loading
################################

# Save it
saveFolder = "/tmp/sasDFG"
saveDFG(fg, saveFolder)

fg2 = initfg();
loadDFG(saveFolder, IncrementalInference,fg2)
tree = solveTree!(fg)

retDFG.solverParams = SolverParams()
#### Comparing
# Enable debugging everywhere
ENV["JULIA_DEBUG"] = "all"
#https://stackoverflow.com/questions/53548681/how-to-enable-debugging-messages-in-juno-julia-editor
d = getData(f).fnc.usrfnc!
dataRet = getData(getFactor(retDFG, f.label)).fnc.usrfnc!
@test Caesar.compare(d, dataRet)

@time predL1 = IncrementalInference.approxConv(fg, ls(fg, :l1)[1], :l1, N=N)s


text:
function get_fourbar(;
    timestep=0.01,
    gravity=[0.0; 0.0; -9.81],
    model=:fourbar,
    spring=0.0,
    damper=1.0,
    T=Float64)

    path = joinpath(@__DIR__, "../deps/$(String(model)).urdf")
    mech = Mechanism(path, false, T,
        gravity=gravity,
        timestep=timestep,
        spring=spring,
        damper=damper)

    return mech
end

function initialize_fourbar!(mechanism::Mechanism;
    angle=0.0,
    angular_velocity=szeros(2)) where T

    zero_velocity!(mechanism)
    set_minimal_coordinates_velocities!(mechanism, get_joint(mechanism, :jointb1);
        xmin=[ -angle, angular_velocity[1]])
    set_minimal_coordinates_velocities!(mechanism, get_joint(mechanism, :joint12);
        xmin=[+2angle, 0])
    set_minimal_coordinates_velocities!(mechanism, get_joint(mechanism, :jointb3);
        xmin=[ +angle, angular_velocity[2]])
    set_minimal_coordinates_velocities!(mechanism, get_joint(mechanism, :joint34);
        xmin=[-2angle, 0])
end


text:
function twodShape(x,r,s,w)
#= """
  twodShape - computes test functions and derivatives on a triangular element
              given the element coordinates and quadrature points.

              ! Note: optimized for straight-sided elements.  Use
              ! `twod_shapeiso' for isoparametric elements.

              local unknown numbering follows the conventions

            3             3             3             3
            |\            |\            |\            87
            | \     or    6 5     or    6 5     or    | \   x is 10
            |  \          |  \          |7 \          9x 6
            1---2         1-4-2         1-4-2         14-52

             P1            P2      Crouzeux-Raviart    P3

  Author: Jeff Borggaard, Virginia Tech
          part of FEMfunctions.jl

  Licensing:
     This code is distributed under the MIT license.

  Usage:
  ```julia
    xg,wg,ϕ,ϕ_x,ϕ_y = twodShape(x,r,s,w)
  ```

  Arguments:
  - `x`: Coordinates of the element nodes
  - `(r,s)`: Coordinates of Gauss points in unit triangle
  - `w`: Gauss weights associated with (r,s)

  Output arguments:
  - `xg`: Coordinates of Gauss points in the element
  - `wg`: Gauss weights scaled by the element Jacobian
  - `ϕ`: Value of element shape functions at xg
  - `ϕ_x`:
  - `ϕ_y`: First spatial derivatives of ϕ

  Changes:
      Crouzeux-Raviart element implemented by Alexander Hay, 2007
""" =#

  n   = size(x,1)
  ng  = size(r,1)

  xg  = zeros(Float64,ng,2)
  ϕ   = zeros(Float64,ng,n)
  ϕ_x = zeros(Float64,ng,n)
  ϕ_y = zeros(Float64,ng,n)

  wg  = zeros(Float64,ng)

  o   = ones(Float64,size(r))

  # Compute (r,s) -> (x,y) transformation for straight-sided elements
  c0 =  x[1,:]
  c1 = -x[1,:] + x[2,:]
  c2 = -x[1,:]          + x[3,:]

  xg[:,1] = c0[1]*o .+ c1[1]*r .+ c2[1]*s
  xr      = c1[1]*o
  xs      = c2[1]*o

  xg[:,2] = c0[2]*o .+ c1[2]*r .+ c2[2]*s
  yr      = c1[2]*o
  ys      = c2[2]*o

  # Compute the Jacobian of the (r,s) -> (x,y) transformation
  jac = xr.*ys .- yr.*xs
  wg  = jac.*w

  rx  = ys./jac
  sx  =-yr./jac
  ry  =-xs./jac
  sy  = xr./jac

  # Compute shape function and derivatives at Gauss points
  if n == 3
    ϕ[:,1] = 1.0*o - r  - s
    ϕ[:,2] =         r
    ϕ[:,3] =              s

    ϕ_x[:,1] =       -rx - sx
    ϕ_x[:,2] =        rx
    ϕ_x[:,3] =             sx

    ϕ_y[:,1] =      - ry - sy
    ϕ_y[:,2] =        ry
    ϕ_y[:,3] =             sy

  elseif n == 6
    ϕ[:,1] = 1.0*o - 3.0*r - 3.0*s + 2.0*r.*r + 4.0*r.*s + 2.0*s.*s
    ϕ[:,2] =       - 1.0*r         + 2.0*r.*r
    ϕ[:,3] =               - 1.0*s                       + 2.0*s.*s
    ϕ[:,4] =         4.0*r         - 4.0*r.*r - 4.0*r.*s
    ϕ[:,5] =                                    4.0*r.*s
    ϕ[:,6] =                 4.0*s            - 4.0*r.*s - 4.0*s.*s

    ϕ_x[:,1] = ( -3.0*o + 4.0*r + 4.0*s ).*rx + ( -3.0*o + 4.0*r + 4.0*s ).*sx
    ϕ_x[:,2] = ( -1.0*o + 4.0*r         ).*rx
    ϕ_x[:,3] =                                  ( -1.0*o         + 4.0*s ).*sx
    ϕ_x[:,4] = (  4.0*o - 8.0*r - 4.0*s ).*rx + (        - 4.0*r         ).*sx
    ϕ_x[:,5] = (                  4.0*s ).*rx + (          4.0*r         ).*sx
    ϕ_x[:,6] = (                - 4.0*s ).*rx + (  4.0*o - 4.0*r - 8.0*s ).*sx

    ϕ_y[:,1] = ( -3.0*o + 4.0*r + 4.0*s ).*ry + ( -3.0*o + 4.0*r + 4.0*s ).*sy
    ϕ_y[:,2] = ( -1.0*o + 4.0*r         ).*ry
    ϕ_y[:,3] =                                  ( -1.0*o         + 4.0*s ).*sy
    ϕ_y[:,4] = (  4.0*o - 8.0*r - 4.0*s ).*ry + (        - 4.0*r         ).*sy
    ϕ_y[:,5] = (                  4.0*s ).*ry + (          4.0*r         ).*sy
    ϕ_y[:,6] = (                - 4.0*s ).*ry + (  4.0*o - 4.0*r - 8.0*s ).*sy

  elseif n == 7
    ϕ[:,1] = (1.0*o-r-s).*(2.0*(1.0-r-s)-1.0) +  3.0*(1.0-r-s).*r.*s
    ϕ[:,2] = r.*(2.0*r-1.0)                   +  3.0*(1.0-r-s).*r.*s
    ϕ[:,3] = s.*(2.0*s-1.0)                   +  3.0*(1.0-r-s).*r.*s
    ϕ[:,4] = 4.0*(1.0-r-s).*r                 - 12.0*(1.0-r-s).*r.*s
    ϕ[:,5] = 4.0*r.*s                         - 12.0*(1.0-r-s).*r.*s
    ϕ[:,6] = 4.0*s.*(1.0-r-s)                 - 12.0*(1.0-r-s).*r.*s
    ϕ[:,7] = 27.0*(1.0-r-s).*r.*s

    ϕ_r = zeros(Float64,ng,n)
    ϕ_r[:,1] = -3.0 + 4.0*r + 7.0*s - 6.0*r.*s - 3.0*(s.^2)
    ϕ_r[:,2] = -1.0 + 4.0*r + 3.0*s - 6.0*r.*s - 3.0*(s.^2)
    ϕ_r[:,3] =                3.0*s - 6.0*r.*s - 3.0*(s.^2)
    ϕ_r[:,4] =  4.0 - 8.0*r -16.0*s +24.0*r.*s +12.0*(s.^2)
    ϕ_r[:,5] =              - 8.0*s +24.0*r.*s +12.0*(s.^2)
    ϕ_r[:,6] =              -16.0*s +24.0*r.*s +12.0*(s.^2)
    ϕ_r[:,7] =               27.0*s -54.0*r.*s -27.0*(s.^2)

    ϕ_s = zeros(Float64,ng,n)
    ϕ_s[:,1] = -3.0 + 7.0*r + 4.0*s - 6.0*r.*s - 3.0*(r.^2)
    ϕ_s[:,2] =        3.0*r         - 6.0*r.*s - 3.0*(r.^2)
    ϕ_s[:,3] = -1.0 + 3.0*r + 4.0*s - 6.0*r.*s - 3.0*(r.^2)
    ϕ_s[:,4] =      -16.0*r         +24.0*r.*s +12.0*(r.^2)
    ϕ_s[:,5] =      - 8.0*r         +24.0*r.*s +12.0*(r.^2)
    ϕ_s[:,6] =  4.0 -16.0*r - 8.0*s +24.0*r.*s +12.0*(r.^2)
    ϕ_s[:,7] =       27.0*r         -54.0*r.*s -27.0*(r.^2)

    ϕ_x = ϕ_r.*rx + ϕ_s.*sx
    ϕ_y = ϕ_r.*ry + ϕ_s.*sy

  elseif ( n==10 )
    # Compute shape function and derivatives at Gauss points
    t = ones(size(r)) - r - s

    ϕ = zeros(Float64,ng,n);
    ϕ[:, 1] =  4.5*(t-1.0/3.0).*(t-2.0/3.0).*t
    ϕ[:, 2] =  4.5*(r-1.0/3.0).*(r-2.0/3.0).*r
    ϕ[:, 3] =  4.5*(s-1.0/3.0).*(s-2.0/3.0).*s
    ϕ[:, 4] =  9.0*r-22.5*r.*(r+s)+13.5*r.*(r+s).^2
    ϕ[:, 5] = -4.5*r+18.0*r.^2+4.5*r.*s-13.5*r.^3-13.5*(r.^2).*s
    ϕ[:, 6] = -4.5*r.*s+13.5*(r.^2).*s
    ϕ[:, 7] = -4.5*r.*s+13.5*r.*s.^2
    ϕ[:, 8] = -4.5*s+4.5*r.*s+18.0*s.^2-13.5*r.*s.^2-13.5*s.^3
    ϕ[:, 9] =  9.0*s-22.5*r.*s-22.5*s.^2+13.5*s.*(r+s).^2
    ϕ[:,10] = 27.0*r.*s-27.0*(r.^2).*s-27.0*r.*s.^2

    ϕ_r = zeros(Float64,ng,n);
    ϕ_r[:, 1] = -5.5 +18.0*(r+s)-13.5*(r+s).^2
    ϕ_r[:, 2] =  1.0 -9.0*r+13.5*r.^2
#   ϕ_r[:, 3] =  zeros(Float64,ng,1)
    ϕ_r[:, 4] =  9.0 -45.0*r-22.5*s+40.5*r.^2+54.0*r.*s+13.5*s.^2
    ϕ_r[:, 5] = -4.5 +36.0*r+4.5*s-40.5*r.^2-27.0*r.*s
    ϕ_r[:, 6] = -4.5*s+27.0*r.*s
    ϕ_r[:, 7] = -4.5*s+13.5*s.*s
    ϕ_r[:, 8] =  4.5*s-13.5*s.*s
    ϕ_r[:, 9] =-22.5*s+27.0*(r.*s+s.*s)
    ϕ_r[:,10] = 27.0*(s - 2.0*r.*s-s.*s)

    ϕ_s = zeros(Float64,ng,n)
    ϕ_s[:, 1] = -5.5+18.0*(r+s)-13.5*(r+s).^2
#   ϕ_s[:, 2] = zeros(Float64,ng,1)
    ϕ_s[:, 3] = 1.0-9.0*s+13.5*s.^2
    ϕ_s[:, 4] = -22.5*r+27.0*(r.^2+r.*s)
    ϕ_s[:, 5] = 4.5*r-13.5*r.^2
    ϕ_s[:, 6] =-4.5*r+13.5*r.^2
    ϕ_s[:, 7] =-4.5*r+27.0*r.*s
    ϕ_s[:, 8] =-4.5+4.5*r+36.0*s-27.0*r.*s-40.5*s.^2
    ϕ_s[:, 9] = 9.0-22.5*r-45.0*s+13.5*r.^2+54.0*r.*s+40.5*s.^2
    ϕ_s[:,10] = 27.0*(r - r.^2 - 2.0*r.*s)

    ϕ_x = ϕ_r.*rx + ϕ_s.*sx
    ϕ_y = ϕ_r.*ry + ϕ_s.*sy

#   else
#     warning('element not supported')
  end

  return xg, wg, ϕ, ϕ_x, ϕ_y

end


text:
# ---
# title: 812. Largest Triangle Area
# id: problem812
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/largest-triangle-area/description/>
# hidden: true
# ---
#
# You have a list of points in the plane. Return the area of the largest
# triangle that can be formed by any 3 of the points.
#
#
#
#     **Example:**
#     Input: points = [[0,0],[0,1],[1,0],[0,2],[2,0]]
#     Output: 2
#     Explanation:
#     The five points are show in the figure below. The red triangle is the largest.
#
#
# ![](https://s3-lc-upload.s3.amazonaws.com/uploads/2018/04/04/1027.png)
#
# **Notes:**
#
#   * `3 <= points.length <= 50`.
#   * No points will be duplicated.
#   *  `-50 <= points[i][j] <= 50`.
#   * Answers within `10^-6` of the true value will be accepted as correct.
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
const DEFAULT_REGISTRY = "General"

# Skip these packages when testing packages
const skip_lists = Dict{String,Vector{String}}()

"""
    prepare_registry([name])

Prepare the given registry for use by this package and update it if `update` is true.
`name` must correspond to an existing stanza in the `deps/Registries.toml` file.
"""
function prepare_registry(name=DEFAULT_REGISTRY; update::Bool=false)
    reg = read_registries()[name]
    regspec = RegistrySpec(name = name, url = reg["url"], uuid = UUID(reg["uuid"]))

    # clone and update the registry
    if !any(existing_regspec -> existing_regspec.name == name, Pkg.Types.collect_registries())
        Pkg.Types.clone_or_cp_registries([regspec])
    elseif update
        Pkg.Registry.update(name)
    end

    # read some metadata
    skip_lists[name] = get(reg, "skip", String[])

    return
end

"""
    read_pkgs([pkgs::Vector{String}]; [registry::String])

Read all packages from a registry and return them as a vector of tuples containing the
package name and registry, its UUID, and a path to it. If `pkgs` is given, only collect
packages matching the names in `pkgs`
"""
function read_pkgs(pkg_names::Vector{String}=String[]; registry=DEFAULT_REGISTRY)
    pkg_names = Set(pkg_names)
    want_all = isempty(pkg_names)

    pkg_data = []
    regpath = registry_path(registry)
    open(joinpath(regpath, "Registry.toml")) do io
        for (_uuid, pkgdata) in Pkg.Types.read_registry(joinpath(regpath, "Registry.toml"))["packages"]
            uuid = UUID(_uuid)
            name = pkgdata["name"]

            if !want_all
                name in pkg_names || continue
                delete!(pkg_names, name)
            end

            path = abspath(regpath, pkgdata["path"])
            push!(pkg_data, (name=name, uuid=uuid, path=path, registry=registry))
        end
    end

    if !want_all && !isempty(pkg_names)
        @warn """did not find the following packages in the $registry registry:\n $("  - " .* join(pkg_names, '\n'))"""
    end

    return pkg_data
end


text:
#=
@testset "combine" begin
    for (v1,v2,v3) in ((OneTo(10), OneTo(9), OneTo(10)),
                       (OneTo(10), 11:15, OneTo(15)),
                       (OneTo(10), 10:15, OneTo(15)),
                       (OneTo(10), 9:15, OneTo(15)),
                       (OneTo(10), -1:15, -1:15),
                       (OneTo(10), 2:9, 1:10),
                       (OneTo(10), 12:15, vcat(OneTo(10), 12:15)),
                       (OneTo(10), 1.5:10.5, sort(vcat(OneTo(10), 1.5:10.5))),
                       (1:5, 1:5, 1:5),
                       (1:6, 4:10, 1:10),
                       (1:5, 6:10, 1:10),
                       (1:5, 7:10, vcat(1:5, 7:10)),
                       (1:7, 6:10, 1:10),
                       (1:11, 6:10, 1:11),
                       (1:1:10, 1:1:9, 1:1:10),
                       (1:1.5:10, 1:1:9, sort(vcat(1:1.5:10, 1:1:10))),
                      )

        @testset "$(v1)-$(v2)" begin
            @test combine(v1, v2) == v3
            @test combine(v2, v1) == v3
        end
        rv1 = reverse(v1)
        @testset "$(rv1)-$(v2)" begin
            @test combine(rv1, v2) == reverse(v3)
            @test combine(v2, rv1) == combine(v2, reverse(rv1))
        end
        #=
        @testset "$(v4)-$(v5)" begin
            @test combine(v4, v5) == v6
            @test combine(v4, v5) == v6
        end
        =#
    end
end
=#


text:
chk = copy(io)
[chk[k] = edit_with(copy(chk[k]), Rename(:value, k)) for k in [:supply, :use]]

#   interm(yr,j(jc_use),"use") = use(yr,"interm",jc_use);
#   basicva(yr,j(jc_use),"use") = use(yr,"basicvalueadded",jc_use);
#   valueadded(yr,j(jc_use),"use") = use(yr,"valueadded",jc_use);
#   output(yr,j(jc_use),"use") = use(yr,"industryoutput",jc_use);
#   totint(yr,i(ir_use),"use") = use(yr,ir_use,"totint");
#   totaluse(yr,i(ir_use),"use") = use(yr,ir_use,"totaluse");
chk[:interm]     = chk[:use] |> @filter(.&(_.i == "interm",          _.j in set[:j]))    |> DataFrame
chk[:basicva]    = chk[:use] |> @filter(.&(_.i == "basicvalueadded", _.j in set[:j]))    |> DataFrame
chk[:valueadded] = chk[:use] |> @filter(.&(_.i == "valueadded",      _.j in set[:j]))    |> DataFrame
chk[:output]     = chk[:use] |> @filter(.&(_.i == "industryoutput",  _.j in set[:j]))    |> DataFrame
chk[:totint]     = chk[:use] |> @filter(.&(_.i in set[:i],           _.j == "totint"))   |> DataFrame
chk[:totaluse]   = chk[:use] |> @filter(.&(_.i in set[:i],           _.j == "totaluse")) |> DataFrame

#   basicsupply(yr,i(ir_supply),'supply') = supply(yr,ir_supply,"BasicSupply");
#   tsupply(yr,i(ir_supply),'supply') = supply(yr,ir_supply,"Supply");
chk[:basicsupply] = chk[:supply] |> @filter(.&(_.i in set[:i], _.j == "basicsupply")) |> DataFrame
chk[:tsupply]     = chk[:supply] |> @filter(.&(_.i in set[:i], _.j == "supply"))      |> DataFrame

# Total intermediate inputs (purchasers' prices)
#   interm(yr,j,'id0') = sum(i,id0(yr,i,j));
#   interm(yr,j,"chk") = interm(yr,j,'id0') - interm(yr,j,"use");
chk[:interm][!, :id0] .= sum_over(chk[:id0], :i)
chk[:interm][!, :chk] .= chk[:interm][!, :id0] - chk[:interm][!, :use]

# Basic value added (purchasers' prices)
#   basicva(yr,j,"va0") = sum(va,va0(yr,va,j));
#   basicva(yr,j,"chk") = basicva(yr,j,"use") - basicva(yr,j,"va0");
chk[:basicva][!, :va0] .= sum_over(chk[:va0], :va)
chk[:basicva][!, :chk] .= chk[:basicva][!, :use] - chk[:basicva][!, :va0]

# Value added (purchaser's prices)
#   valueadded(yr,j,"va0+ts0") = sum(va,va0(yr,va,j)) + ts0(yr,"taxes",j) - ts0(yr,"subsidies",j);
#   valueadded(yr,j,"chk") = valueadded(yr,j,"use") - valueadded(yr,j,"va0+ts0");
chk[:valueadded][!, :va0_ts0] .=
    sum_over(chk[:va0], :va) +
    chk[:ts0][chk[:ts0][:,:i] .== "taxes", :value] -
    chk[:ts0][chk[:ts0][:,:i] .== "subsidies", :value]
chk[:valueadded][!, :chk] .= chk[:valueadded][!, :use] - chk[:valueadded][!, :va0_ts0]

# Check on total taxes
#   taxtotal(yr,"ts_subsidies") = sum(j, ts0(yr,"subsidies",j));
#   taxtotal(yr,"ts_taxes") = sum(j, ts0(yr,"taxes",j));
#   taxtotal(yr,"s0") = sum(i,sbd0(yr,i));
#   taxtotal(yr,"t0+duty") = sum(i,tax0(yr,i)+duty0(yr,i));
chk[:taxtotal] = DataFrame(yr = set[:yr])
chk[:taxtotal][!, :ts_subsidies] .= sum_over(chk[:ts0][chk[:ts0][:,:i] .== "subsidies", :], :j)
chk[:taxtotal][!, :ts_taxes] .= sum_over(chk[:ts0][chk[:ts0][:,:i] .== "taxes", :], :j)
chk[:taxtotal][!, :s0] .= sum_over(chk[:sbd0], :i)
chk[:taxtotal][!, :t0_duty] .= sum_over(chk[:tax0], :i) + sum_over(chk[:duty0], :i)

# Total industry output (basic prices)
#   output(yr,j,"id0+va0") = sum(va,va0(yr,va,j)) + sum(i,id0(yr,i,j));
#   output(yr,j,'ys0') = sum(i,ys0(yr,j,i));
#   output(yr,j,"chk") = output(yr,j,"id0+va0") - output(yr,j,"use");
#   output(yr,j,"chk-ys0") = output(yr,j,"id0+va0") - output(yr,j,'ys0');
chk[:output][!, :id0_va0] .= sum_over(chk[:va0], :va) + sum_over(chk[:id0], :i)
chk[:output][!, :ys0]     .= sum_over(chk[:ys0], :i)
chk[:output][!, :chk]     .= chk[:output][:, :id0_va0] - chk[:output][:, :use]
chk[:output][!, :chk_ys0] .= chk[:output][:, :id0_va0] - chk[:output][:, :ys0]

# Total intermediate use (purchasers' prices)
#   totint(yr,i,'id0') = sum(j,id0(yr,i,j));
#   totint(yr,i,"chk") = totint(yr,i,"use") - totint(yr,i,'id0');
chk[:totint][!, :id0] .= sum_over(chk[:id0], :j)
chk[:totint][!, :chk] .= chk[:totint][:, :use] - chk[:totint][:, :id0]

# Total use of commodities (purchasers' prices)
#   totaluse(yr,i,"id0+fd0") = sum(j,id0(yr,i,j)) + sum(fd,fd0(yr,i,fd)) + x0(yr,i);
#   totaluse(yr,i,"chk") = totaluse(yr,i,"use") - totaluse(yr,i,"id0+fd0");
chk[:totaluse][!, :id0_fd0] .= sum_over(chk[:id0], :j) + sum_over(chk[:fd0], :fd) + chk[:x0][:,:value]
chk[:totaluse][!, :chk]     .= chk[:totaluse][:, :use] - chk[:totaluse][:, :id0_fd0]

# Basic supply
#   basicsupply(yr,i,'ys0') = sum(j,ys0(yr,j,i));
#   basicsupply(yr,i,"chk") = basicsupply(yr,i,'supply') - basicsupply(yr,i,'ys0');
chk[:basicsupply][!, :ys0] .= sum_over(chk[:ys0], :j)
chk[:basicsupply][!, :chk] .= chk[:basicsupply][:, :supply] - chk[:basicsupply][:, :ys0]

# Total supply
#   tsupply(yr,i,"ys0+...") = sum(j,ys0(yr,j,i)) + m0(yr,i) + mrg0(yr,i) +
#       trn0(yr,i) + duty0(yr,i) + tax0(yr,i) - sbd0(yr,i);
#   tsupply(yr,i,"totaluse") = sum(j,id0(yr,i,j)) + sum(fd,fd0(yr,i,fd)) + x0(yr,i);
#   tsupply(yr,i,"chk") = tsupply(yr,i,'supply') - tsupply(yr,i,"totaluse");
#   tsupply(yr,i,"supply-use") = tsupply(yr,i,"ys0+...") - tsupply(yr,i,"totaluse");
chk[:tsupply][!, :ys0_etc] .= sum_over(chk[:ys0], :j) + chk[:m0][:,:value] + chk[:mrg0][:,:value] +
    chk[:trn0][:,:value] + chk[:duty0][:,:value] + chk[:tax0][:,:value] - chk[:sbd0][:,:value]
chk[:tsupply][!, :totaluse] .= sum_over(chk[:id0], :j) + sum_over(chk[:fd0], :fd) + chk[:x0][:,:value]
chk[:tsupply][!, :chk]        .= chk[:tsupply][:, :supply] - chk[:tsupply][:, :totaluse]
chk[:tsupply][!, :supply_use] .= chk[:tsupply][:, :ys0_etc] - chk[:tsupply][:, :totaluse]

# Check on accounting identities.
#   details(yr,i,"y0") = y0(yr,i);
#   details(yr,i,"m0") = m0(yr,i) + duty0(yr,i);
#   details(yr,i,"mrg+trn") = mrg0(yr,i) + trn0(yr,i);
#   details(yr,i,"tax-sbd") = tax0(yr,i)-sbd0(yr,i);
#   details(yr,i,'id0') = sum(j, id0(yr,i,j));
#   details(yr,i,"fd0") = sum(fd,fd0(yr,i,fd));
#   details(yr,i,"x0") = x0(yr,i);
#   details(yr,i,"balance") = y0(yr,i)+m0(yr,i)+duty0(yr,i) + tax0(yr,i)-sbd0(yr,i)
#       - sum(j, id0(yr,i,j)) - sum(fd,fd0(yr,i,fd)) - x0(yr,i)
# 		+ (mrg0(yr,i) + trn0(yr,i));
chk[:details] = DataFrame(permute((yr = set[:yr], i = set[:i])))
chk[:details][!,:y0] .= chk[:y0][:,:value]
chk[:details][!,:m0] .= chk[:m0][:,:value] + chk[:duty0][:,:value]
chk[:details][!,:mrg_trn] .= chk[:mrg0][:,:value] + chk[:trn0][:,:value]
chk[:details][!,:tax_sbd] .= chk[:tax0][:,:value] - chk[:sbd0][:,:value]
chk[:details][!,:id0] .= sum_over(chk[:id0], :j)
chk[:details][!,:fd0] .= sum_over(chk[:fd0], :fd)
chk[:details][!,:x0]  .= chk[:x0][:,:value]
chk[:details][!,:balance] .=
    chk[:y0][:,:value] + chk[:m0][:,:value] + chk[:duty0][:,:value] + chk[:tax0][:,:value] - chk[:sbd0][:,:value] -
    sum_over(chk[:id0], :j) - sum_over(chk[:fd0], :fd) - chk[:x0][:,:value] +
    chk[:mrg0][:,:value] + chk[:trn0][:,:value];

# Report of margin producing sectors.
    # imrginfo(yr,"a0",imrg) = a0(yr,imrg);
    # imrginfo(yr,"tax0",imrg) = tax0(yr,imrg);
    # imrginfo(yr,"sbd0",imrg) = sbd0(yr,imrg);
    # imrginfo(yr,"y0",imrg) = y0(yr,imrg);
    # imrginfo(yr,"x0",imrg) = x0(yr,imrg);
    # imrginfo(yr,"m0",imrg) = m0(yr,imrg);
    # imrginfo(yr,"duty0",imrg) = duty0(yr,imrg);
    # imrginfo(yr,m,imrg) = md0(yr,m,imrg);
chk[:imrg] = fill_zero((yr = set[:yr],))

# Only keep keys unique to the "chk" Dictionary.
chk = Dict(k => chk[k] for k in setdiff(keys(chk), keys(io)))

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
"""
# OptimizeModel

Create a OptimizeModel.

### Required arguments
```julia
* `name::AbstractString`        : Name for the model
* `model::AbstractString`       : Stan model source
```

### Optional arguments
```julia
* `n_chains::Vector{Int64}=[4]`        : Optionally updated in stan_sample()
* `seed::RandomSeed`                     : Random seed settings
* `output::Output`              : File output options
* `init::Init`                         : Default interval bound for parameters
* `tmpdir::AbstractString`             : Directory where output files are stored
* `output_base::AbstractString`        : Base name for output files
* `exec_path::AbstractString`          : Path to cmdstan executable
* `data_file::vector{AbstractString}`  : Path to per chain data file
* `init_file::Vector{AbstractString}`  : Path to per chain init file
* `cmds::Vector{Cmd}`                  : Path to per chain init file
* `sample_file::Vector{String}         : Path to per chain samples file
* `log_file::Vector{String}            : Path to per chain log file
* `diagnostic_file::Vector{String}    : Path to per chain diagnostic file
* `summary=true`                       : Create computed stan summary
* `printsummary=true`                  : Show computed stan summary
* `sm::StanRun.StanModel`              : StanRun.StanModel
* `method::Optimize`                   : Fix Optimize
```

"""
mutable struct OptimizeModel <: CmdStanModels
  @shared_fields_stanmodels
  method::Optimize
end

function OptimizeModel(
  name::AbstractString,
  model::AbstractString,
  n_chains=[4];
  method = Optimize(),
  seed = StanBase.RandomSeed(),
  init = StanBase.Init(),
  output = StanBase.Output(),
  tmpdir = mktempdir())

  !isdir(tmpdir) && mkdir(tmpdir)

  StanBase.update_model_file(joinpath(tmpdir, "$(name).stan"), strip(model))
  sm = StanModel(joinpath(tmpdir, "$(name).stan"))

  output_base = StanRun.default_output_base(sm)
  exec_path = StanRun.ensure_executable(sm)

  stan_compile(sm)

  OptimizeModel(name, model, n_chains, seed, init, output,
    tmpdir, output_base, exec_path, String[], String[],
    Cmd[], String[], String[], String[], false, false, sm, method)
end

function optimize_show(io::IO, m::OptimizeModel, compact::Bool)
  println(io, "  name =                    \"$(m.name)\"")
  println(io, "  n_chains =                $(StanBase.get_n_chains(m))")
  println(io, "  output =                  Output()")
  println(io, "    file =                    \"$(split(m.output.file, "/")[end])\"")
  println(io, "    diagnostics_file =        \"$(split(m.output.diagnostic_file, "/")[end])\"")
  println(io, "    refresh =                 $(m.output.refresh)")
  println(io, "  tmpdir =                  \"$(m.tmpdir)\"")
  optimize_show(io, m.method, compact)
end

show(io::IO, m::OptimizeModel) = optimize_show(io, m, false)


text:
####
# запустить: ]dev TestFilter

using TestFilter
using DSP
using DelimitedFiles
using Plots

x = readdlm("examples/ecg.csv")[:, 1]

f = DSPFilter()
fcut = 30
fs = 257
init(f, Lowpass(2 * fcut / fs), Butterworth(2))

y = similar(x)
run_filter(f, y, x)

plot(x)
plot!(y)


text:
import Base.zero

zero(::Type{Tuple{UInt16,UInt16,UInt16}}) = (zero(UInt16),zero(UInt16),zero(UInt16))
zero(::Type{Tuple{UInt32,UInt32,UInt32,UInt32}}) = (zero(UInt32),zero(UInt32),zero(UInt32),zero(UInt32))

zero(::Type{hsa_dim3_t}) = hsa_dim3_t(0,0,0)
zero(::Type{hsa_isa_t}) = hsa_isa_t(0)
zero(::Type{hsa_ext_control_directives_t}) =
hsa_ext_control_directives_t(
0,0,0,0,0,0,0,
zero(Array_3_UInt64),
zero(hsa_dim3_t), 0,
zero(Array_75_UInt8)
)

import Base.convert
function convert(::Type{hsa_dim3_t}, values :: Vector{UInt32})
    if length(values) != 3
        error("invalid input argument)")
    end

    hsa_dim3_t(values...)
end

function status_string(status::hsa_status_t)
    string = Ptr{UInt8}[0]

    err = ccall((:hsa_status_string,libhsa),hsa_status_t,(hsa_status_t,Ptr{Ptr{UInt8}}),status,string)

    # manually test for error
    if err & HSA_STATUS_ERROR == HSA_STATUS_ERROR
        # throw regular error to avoid recursion to HSAError
        error("could not retrieve status string")
    end

    return bytestring(string[1])
end


text:
"""
  bar_plot(OperationModelResults)

This function plots a bar plot for the generators in each variable within
the results variables dictionary, and makes a bar plot for all of the variables.

#Examples

results = solve_op_model!(OpModel)
bar_plot(results)

generates a bar plot for each variable,
and one bar plot of all the variables

"""

function bar_plot(res::OperationModelResults)

  key_name = string.(collect(keys(res.variables)))

  for i in 1:length(key_name)

    variable_bar = get_bar_plot_data(res, key_name[i])
    p = RecipesBase.plot(variable_bar, key_name[i])
    display(p)

  end

  bar_gen = get_bar_gen_data(res)
  p2 = RecipesBase.plot(bar_gen)
  display(p2)

end

"""

  stack_plot(OperationModelResults)

This function plots a stack plot for the generators in each variable within
the results variables dictionary, and makes a stack plot for all of the variables.

#Examples

results = solve_op_model!(OpModel)
stack_plot(results)

generates a stack plot for each variable,
and one stack plot of all the variables

"""

function stack_plot(res::OperationModelResults)

  key_name = string.(collect(keys(res.variables)))

  for i in 1:length(key_name)

    variable_stack = get_stacked_plot_data(res, key_name[i])
    p3 = RecipesBase.plot(variable_stack, key_name[i])
    display(p3)

  end

  stacked_gen = get_stacked_generation_data(res)
  p4 = RecipesBase.plot(stacked_gen)
  display(p4)

end
"""
    Recipe StackedPlot(results::StackedArea, variable::String)

This is a Plots series recipe expecting to use Plots and the gr() backend.
RecipesBase.@Recipe function StackedPlot(res::StackedArea)

This recipe is called when plot() takes in data of struct StackedArea.
It creates a stacked area plot. To overlay a scatter plot or line, add
another series:

  Recipesbase.@series begin
    seriestype := (plot type, such as :scatter)
    x, y
  end

  plot attributes can be changed and spacing for the x-tick can be changed.
  For example, if a plot is a 7 day period, xtick can be changed to
  xtick := time[1]:Dates.Day(1):time[n-1]

"""

RecipesBase.@recipe function StackedPlot(results::StackedArea, variable::String)

  time = convert.(Dates.DateTime,results.time_range)
  n = length(time)
  data = results.data_matrix
  z = cumsum(data, dims = 2)
  # Plot attributes
  grid := false
  title := variable
  label := results.labels
  legend := :topleft
  alpha := 0.6
  seriescolor := [:lightblue :darkorange :lightgreen :red :turquoise :blue :orange]
  time_interval = Dates.Hour(convert(Dates.DateTime,time[n])-convert(Dates.DateTime,time[1]))
  xlabel := "$time_interval"
  ylabel := "Generation (MW)"
  xtick := time[1]:Dates.Hour(12):time[n-1]

    #create filled polygon
    sy = vcat(z[:,1],zeros(n-1))
    sx = [time[1:n-1]; reverse(time[1:n-1])]

     for c=1:size(z,2)
      if c !== 1
        sy = hcat(sy,vcat(z[:,c],reverse(z[:,c-1])))
      end
    end

  RecipesBase.@series begin

    seriestype := :shape
    sx, sy

  end

end

RecipesBase.@recipe function StackedGeneration(res::StackedGeneration)

  time = convert.(Dates.DateTime,res.time_range)
  n = length(time)
  data = res.data_matrix
  z = cumsum(data, dims = 2)
  # Plot Attributes
  grid := false
  title := "Generation Type"
  alpha := 0.6
  seriescolor := [:lightblue :orange :lightgreen :red :turquoise]
  label := res.labels
  legend := :bottomright
  time_interval = Dates.Hour(convert(Dates.DateTime,time[n])-convert(Dates.DateTime,time[1]))
  xlabel := "$time_interval"
  ylabel := "Generation (MW)"
  xtick := time[1]:Dates.Hour(12):time[n-1]

  # Create filled polygon
 sy = vcat(z[:,1],zeros(n-1))
 sx = [time[1:n-1]; reverse(time[1:n-1])]

for c=1:size(z,2)
   if c !== 1
     sy = hcat(sy,vcat(z[:,c],reverse(z[:,c-1])))
   end
end

  RecipesBase.@series begin

  seriestype := :shape
   sx, sy

  end

end

RecipesBase.@recipe function BarPlot(res::BarPlot, variable::String)

  time = convert.(Dates.DateTime,res.time_range)
 n = length(time)
 data_point = res.bar_data
 data = [data_point; data_point]
 z = cumsum(data, dims = 2)
 # Plot Attributes
 grid := false
 title := variable
 seriestype := :shape
 label := res.labels
 start_time = time[1]
 time_interval = Dates.Hour(convert(Dates.DateTime,time[n])-convert(Dates.DateTime,time[1]))
 xlabel := "$time_interval, $start_time"
 ylabel := "Generation(MW)"
 alpha := 0.6
 seriescolor := [:lightblue :orange :lightgreen :red :turquoise :blue]
 xlims := (1, 8)
 xticks := false
 n = 2

    # Create filled polygon

  for c=1:size(z,2)
   sx = [[4,5]; [5,4]]
   sy = vcat(z[:,c], c==1 ? zeros(n) : reverse(z[:,c-1]))
   RecipesBase.@series sx, sy
 end

end

RecipesBase.@recipe function BarGen(res::BarGeneration)

  time = convert.(Dates.DateTime,res.time_range)
  n = 2
  data_point = res.bar_data
  data = [data_point; data_point]
  z = cumsum(data, dims = 2)
   # Plot Attributes
   grid := false
   title := "Generation Type"
   seriestype := :shape
   label := res.labels
   start_time = time[1]
   xticks := false
   alpha := 0.6
   seriescolor := [:lightblue :orange :lightgreen :red :turquoise :blue]
   xlims := (1, 8)

   for c=1:size(z,2)
    sx = [[4,5]; [5,4]]
    sy = vcat(z[:,c], c==1 ? zeros(n) : reverse(z[:,c-1]))
    RecipesBase.@series sx, sy
  end
end


text:
module BenchmarkTooling

using Reexport

@reexport using BenchmarkTools

#=
    time_tolerance: The noise tolerance for the benchmark's time estimate, as a percentage.
        This is utilized after benchmark execution, when analyzing results.
        Defaults to BenchmarkTools.DEFAULT_PARAMETERS.time_tolerance = 0.05
    memory_tolerance: The noise tolerance for the benchmark's memory estimate, as a percentage.
        This is utilized after benchmark execution, when analyzing results.
        Defaults to BenchmarkTools.DEFAULT_PARAMETERS.memory_tolerance = 0.01.
=#

BenchmarkTools.DEFAULT_PARAMETERS.time_tolerance = 0.005
BenchmarkTools.DEFAULT_PARAMETERS.memory_tolerance = 0.001

end # module


text:
using Test
import MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities

include("../model.jl")

mock = MOIU.MockOptimizer(Model{Float64}())
config = MOIT.TestConfig()

@testset "QP" begin
    MOIU.set_mock_optimize!(mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [4/7, 3/7, 6/7]))
    MOIT.qp1test(mock, config)
    MOIT.qp2test(mock, config)
    MOIU.set_mock_optimize!(mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1/4, 3/4]),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 0.0]))
    MOIT.qp3test(mock, config)
end
@testset "QCP" begin
    MOIU.set_mock_optimize!(mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1/2, 7/4], MOI.FEASIBLE_POINT))
    MOIT.qcp1test(mock, config)
    MOIU.set_mock_optimize!(mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [√2], MOI.FEASIBLE_POINT))
    MOIT.qcp2test(mock, config)
    MOIU.set_mock_optimize!(mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [√2], MOI.FEASIBLE_POINT))
    MOIT.qcp3test(mock, config)
    MOIU.set_mock_optimize!(mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 1.0], MOI.FEASIBLE_POINT))
    MOIT.qcp4test(mock, config)
end
@testset "Non-convex QCP" begin
    MOIU.set_mock_optimize!(mock,
                            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [4.0, 1.0], MOI.FEASIBLE_POINT))
    MOIT.ncqcp1test(mock, config)
    MOIU.set_mock_optimize!(mock,
                            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [2.0, 2.0], MOI.FEASIBLE_POINT))
    MOIT.ncqcp2test(mock, config)
end
@testset "SOCP" begin
    MOIU.set_mock_optimize!(mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1/2, 1/2, 1/√2]))
    MOIT.socp1test(mock, config)
end


text:
function make_function_dataframe(symbols::Array{Symbol,1}, function_array::Union{Array{Any,1},Array{Function,1},Array{NamedArray,1},Array{Union{NamedArray,Function},1}}, types::Array{Symbol,1}) where T<:Real
    if sum(isa.(function_array, Union{Real,Function,NamedArray})) < length(function_array) error("Function array can only contain reals, Functions or NamedArrays") end
    return DataFrame(name = symbols,
                   func = function_array,
                   type = types)
end


text:
using VectorLogging
using Logging
using Test

@testset "VectorLogging.jl" begin
    logger = VectorLogger()
    with_logger(logger) do
        @info "First log message"
        @warn("Second log message, with payload", payload="foo")
    end
    # Firsdt log entry:
    @test logger.log[1].message == "First log message"
    @test logger.log[1].level == Logging.Info
    # Second log entry:
    @test logger.log[2].message == "Second log message, with payload"
    @test logger.log[2].level == Logging.Warn
    @test logger.log[2].keys[:payload] == "foo"
end


text:
#-------------------------------------
#* Model of a packed column stage
#-------------------------------------
type packedStage
	packedStage()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number Of Components"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Resistance coefficient on the liquid load",
				:Default=>1,
				:Hidden=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"surface area per packing volume",
				:Unit=>"m^2/m^3",
				:Hidden=>true
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Column diameter",
				:Hidden=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Void fraction of packing, (m^3 void space/m^3 packed bed)",
				:Hidden=>true
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Height of packing",
				:Hidden=>true
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of Stages",
				:Default=>3,
				:Hidden=>true
			)),
			heat_rate(Dict{Symbol,Any}(
				:Brief=>"Rate of heat supply",
				:Hidden=>true
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"The Height Equivalent to a Theoretical Plate",
				:Hidden=>true
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Column Sectional Cross Area",
				:Hidden=>true
			)),
			volume(Dict{Symbol,Any}(
				:Brief=>"Total Volume of the tray",
				:Hidden=>true
			)),
			constant(Dict{Symbol,Any}(
				:Brief=>"Pi Number",
				:Default=>3.14159265,
				:Symbol=>"\\pi",
				:Hidden=>true
			)),
			acceleration(Dict{Symbol,Any}(
				:Brief=>"Gravity Acceleration",
				:Default=>9.81,
				:Hidden=>true
			)),
			flow_mol(Dict{Symbol,Any}(
				:Brief=>"Low Flow",
				:Default=>1E-6,
				:Hidden=>true
			)),
			pressure(Dict{Symbol,Any}(
				:Brief=>"Low Pressure",
				:Default=>1E-6,
				:Hidden=>true
			)),
			flow_mol(Dict{Symbol,Any}(
				:Brief=>"No Flow",
				:Default=>0,
				:Hidden=>true
			)),
			fill(molweight(Dict{Symbol,Any}(
				:Brief=>"Component Mol Weight",
				:Hidden=>true
			)),(NComp)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Vapour Flow",
				:Valid=>["on", "off"],
				:Default=>"on",
				:Hidden=>true
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Feed stream",
				:Symbol=>"_{in}",
				:Protected=>true
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet liquid stream",
				:Symbol=>"_{inL}",
				:Protected=>true
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet vapour stream",
				:Symbol=>"_{inV}",
				:Protected=>true
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Outlet liquid stream",
				:Symbol=>"_{outL}",
				:Protected=>true
			)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"Outlet vapour stream",
				:Symbol=>"_{outV}",
				:Protected=>true
			)),
			fill(mol(Dict{Symbol,Any}(
				:Brief=>"Molar Holdup in the tray",
				:Default=>0.01,
				:Lower=>-0.000001,
				:Upper=>100,
				:Protected=>true
			)),(NComp)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Molar liquid holdup",
				:Default=>0.01,
				:Lower=>0,
				:Upper=>100,
				:Protected=>true
			)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Molar vapour holdup",
				:Default=>0.01,
				:Lower=>0,
				:Upper=>100,
				:Protected=>true
			)),
			energy(Dict{Symbol,Any}(
				:Brief=>"Total Energy Holdup on tray",
				:Default=>-500,
				:Protected=>true
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Liquid Molar Volume",
				:Protected=>true
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Vapour Molar volume",
				:Protected=>true
			)),
			viscosity(Dict{Symbol,Any}(
				:Brief=>"Liquid dynamic viscosity",
				:DisplayUnit=>"kg/m/s",
				:Protected=>true
			)),
			dens_mass(Dict{Symbol,Any}(
				:Brief=>"Liquid mass density",
				:Protected=>true
			)),
			dens_mass(Dict{Symbol,Any}(
				:Brief=>"Vapour mass density",
				:Protected=>true
			)),
			velocity(Dict{Symbol,Any}(
				:Brief=>"volume flow rate of liquid, m^3/m^2/s",
				:Lower=>-10,
				:Upper=>1000,
				:Protected=>true
			)),
			velocity(Dict{Symbol,Any}(
				:Brief=>"volume flow rate of vapor, m^3/m^2/s",
				:Lower=>-10,
				:Upper=>1000,
				:Protected=>true
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Area occupied by the liquid",
				:Default=>0.001,
				:Upper=>10,
				:Protected=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Column holdup",
				:Unit=>"m^3/m^3",
				:Default=>0.01,
				:Upper=>10,
				:Protected=>true
			)),
			pressure(Dict{Symbol,Any}(
				:Brief=>"Stage Pressure drop",
				:Protected=>true
			)),
			[
				:(deltaP/HETP = LiquidResistanceCoeff *( 0.5*AreaPerPackingVolume + 2/ColumnInternalDiameter) * 1/((PackingVoidFraction-hl)^3) * (uV^2) *rhoV),
				:(InletVapour.F = zero_flow),
				:(diff(E) = (Inlet.F*Inlet.h + InletLiquid.F*InletLiquid.h + InletVapour.F*InletVapour.h- OutletLiquid.F*OutletLiquid.h - OutletVapour.F*OutletVapour.h + HeatOnStage )),
				:(E = ML*OutletLiquid.h + MV*OutletVapour.h - OutletLiquid.P*V),
				:(V*PackingVoidFraction= ML*vL + MV*vV),
				:(uV * (V*PackingVoidFraction/HETP - Al) = InletVapour.F * vV),
				:(hl*V*PackingVoidFraction = ML*vL),
				:(hl^3 = (12/Gconst) * AreaPerPackingVolume^2 * (miL/rhoL) * uL),
				:(Al = ML*vL/HETP),
				:(diff(M)=Inlet.F*Inlet.z + InletLiquid.F*InletLiquid.z + InletVapour.F*InletVapour.z- OutletLiquid.F*OutletLiquid.z - OutletVapour.F*OutletVapour.z),
				:(M = ML*OutletLiquid.z + MV*OutletVapour.z),
				:(sum(OutletLiquid.z)= 1.0),
				:(sum(OutletLiquid.z)= sum(OutletVapour.z)),
				:(vL = PP.LiquidVolume(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z)),
				:(vV = PP.VapourVolume(OutletVapour.T, OutletVapour.P, OutletVapour.z)),
				:(PP.LiquidFugacityCoefficient(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z)*OutletLiquid.z = PP.VapourFugacityCoefficient(OutletVapour.T, OutletVapour.P, OutletVapour.z)*OutletVapour.z),
				:(OutletVapour.T = OutletLiquid.T),
				:(OutletLiquid.P = OutletVapour.P),
				:(rhoL = PP.LiquidDensity(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z)),
				:(rhoV = PP.VapourDensity(InletVapour.T, InletVapour.P, InletVapour.z)),
				:(miL = PP.LiquidViscosity(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z)),
				:(uL * Al = OutletLiquid.F * vL),
				:(deltaP = InletVapour.P - OutletVapour.P),
			],
			[
				"Pressure drop and Vapor flow, Billet (4-58)","Vapour Flow","Energy Balance","Energy Holdup","Geometry Constraint","Volume flow rate of vapor, m^3/m^2/s","Liquid holdup","Liquid velocity as a function of liquid holdup, Billet (4-27)","Area occupied by the liquid","Component Molar Balance","Molar Holdup","Mol Fraction Normalisation","Mol Fraction Constraint","Liquid Volume","Vapour Volume","Chemical Equilibrium","Thermal Equilibrium","Mechanical Equilibrium","Liquid Density","Vapour Density","Liquid viscosity","Volume flow rate of liquid, m^3/m^2/s","Pressure Drop",
			],
			[:PP,:NComp,:LiquidResistanceCoeff,:AreaPerPackingVolume,:ColumnInternalDiameter,:PackingVoidFraction,:HeightOfPacking,:Number_Stages,:HeatOnStage,:HETP,:ColumnArea,:V,:Pi,:Gconst,:low_flow,:low_pressure,:zero_flow,:Mw,:VapourFlow,],
			[:Inlet,:InletLiquid,:InletVapour,:OutletLiquid,:OutletVapour,:M,:ML,:MV,:E,:vL,:vV,:miL,:rhoL,:rhoV,:uL,:uV,:Al,:hl,:deltaP,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	LiquidResistanceCoeff::positive
	AreaPerPackingVolume::DanaReal
	ColumnInternalDiameter::length
	PackingVoidFraction::DanaReal
	HeightOfPacking::length
	Number_Stages::DanaInteger
	HeatOnStage::heat_rate
	HETP::length
	ColumnArea::area
	V::volume
	Pi::constant
	Gconst::acceleration
	low_flow::flow_mol
	low_pressure::pressure
	zero_flow::flow_mol
	Mw::Array{molweight}
	VapourFlow::DanaSwitcher
	Inlet::stream
	InletLiquid::stream
	InletVapour::stream
	OutletLiquid::liquid_stream
	OutletVapour::vapour_stream
	M::Array{mol}
	ML::mol
	MV::mol
	E::energy
	vL::volume_mol
	vV::volume_mol
	miL::viscosity
	rhoL::dens_mass
	rhoV::dens_mass
	uL::velocity
	uV::velocity
	Al::area
	hl::positive
	deltaP::pressure
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export packedStage
function set(in::packedStage)
	Mw = PP.MolecularWeight()
	 ColumnArea = 0.25*Pi*ColumnInternalDiameter^2
	 HETP = HeightOfPacking/Number_Stages
	 V = HETP * ColumnArea
	 low_pressure = 1E-4 * "atm"
	 low_flow = 1E-6 * "kmol/h"
	 zero_flow = 0 * "kmol/h"

end
function setEquationFlow(in::packedStage)
	let switch=VapourFlow
		if VapourFlow==InletVapour.F < low_flow
			set(switch,"off")
		end
		if VapourFlow==deltaP > low_pressure
			set(switch,"on")
		end
		if switch=="on"
			addEquation(1)
		elseif switch=="off"
			addEquation(2)
		end
	end
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
	addEquation(12)
	addEquation(13)
	addEquation(14)
	addEquation(15)
	addEquation(16)
	addEquation(17)
	addEquation(18)
	addEquation(19)
	addEquation(20)
	addEquation(21)
	addEquation(22)
	addEquation(23)
end
function atributes(in::packedStage,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="Complete model of a packed column stage."
	fields[:Info]="== Specify ==
* the Feed stream
* the Liquid inlet stream
* the Vapour inlet stream
* the stage pressure drop (deltaP)

== Initial ==
* the plate temperature (OutletLiquid.T)
* the liquid molar holdup ML
* (NoComps - 1) OutletLiquid compositions
"
	drive!(fields,_)
	return fields
end
packedStage(_::Dict{Symbol,Any})=begin
	newModel=packedStage()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
module OutlierDetectionBenchmark
    using JLD
    using BenchmarkTools
    using Random: MersenneTwister

    export generate_report_single,
           generate_report_comparison,
           prepare_benchmarks,
           run_benchmarks

    include("benchmarks.jl")
end


text:
using Test
using JustSayIt
using JustSayIt.API
using PyCall
import JustSayIt: DEFAULT_MODEL_NAME, TYPE_MODEL_NAME, MODELDIR_PREFIX, DEFAULT_NOISES, AUDIO_ELTYPE, COMMAND_RECOGNIZER_ID
import JustSayIt: init_jsi, finalize_jsi, recognizer, noises, reader, start_reading, stop_reading, read_wav, set_default_streamer, reset_all, next_partial_recognition, next_recognition, next_token, _is_next, is_next, _are_next, are_next

# Test setup
const SAMPLEDIR_CMD     = joinpath("samples", "commands")
const SAMPLEDIR_SILENCE = joinpath("samples", "silence")

commands = Dict("help"      => Help.help,
                "type"      => Keyboard.type,
                "ma"        => Mouse.click_left,
                "select"    => Mouse.press_left,
                "okay"      => Mouse.release_left,
                "middle"    => Mouse.click_middle,
                "right"     => Mouse.click_right,
                "double"    => Mouse.click_double,
                "triple"    => Mouse.click_triple,
                "email"     => Email.email,
                "internet"  => Internet.internet,
                "copy"      => (Key.ctrl, 'c'),
                "cut"       => (Key.ctrl, 'x'),
                "paste"     => (Key.ctrl, 'v'),
                "undo"      => (Key.ctrl, 'z'),
                "redo"      => (Key.ctrl, Key.shift, 'z'),
                "page up"   => Key.page_up,
                "page down" => Key.page_down,
                );
modeldirs = Dict(DEFAULT_MODEL_NAME => joinpath(MODELDIR_PREFIX, "vosk-model-small-en-us-0.15"),
                 TYPE_MODEL_NAME    => joinpath(MODELDIR_PREFIX, "vosk-model-small-en-us-0.15"))
init_jsi(commands, modeldirs, DEFAULT_NOISES)

samples = Dict("help"      => read_wav(joinpath(SAMPLEDIR_CMD, "help.wav")),
               "type"      => read_wav(joinpath(SAMPLEDIR_CMD, "type.wav")),
               "ma"        => read_wav(joinpath(SAMPLEDIR_CMD, "ma.wav")),
               "select"    => read_wav(joinpath(SAMPLEDIR_CMD, "select.wav")),
               "okay"      => read_wav(joinpath(SAMPLEDIR_CMD, "okay.wav")),
               "middle"    => read_wav(joinpath(SAMPLEDIR_CMD, "middle.wav")),
               "right"     => read_wav(joinpath(SAMPLEDIR_CMD, "right.wav")),
               "double"    => read_wav(joinpath(SAMPLEDIR_CMD, "double.wav")),
               "triple"    => read_wav(joinpath(SAMPLEDIR_CMD, "triple.wav")),
               "email"     => read_wav(joinpath(SAMPLEDIR_CMD, "email.wav")),
               "internet"  => read_wav(joinpath(SAMPLEDIR_CMD, "internet.wav")),
               "copy"      => read_wav(joinpath(SAMPLEDIR_CMD, "copy.wav")),
               "cut"       => read_wav(joinpath(SAMPLEDIR_CMD, "cut.wav")),
               "paste"     => read_wav(joinpath(SAMPLEDIR_CMD, "paste.wav")),
               "undo"      => read_wav(joinpath(SAMPLEDIR_CMD, "undo.wav")),
               "redo"      => read_wav(joinpath(SAMPLEDIR_CMD, "redo.wav")),
               "page up"   => read_wav(joinpath(SAMPLEDIR_CMD, "page_up.wav")),
               "page down" => read_wav(joinpath(SAMPLEDIR_CMD, "page_down.wav")),
               );
_05 = read_wav(joinpath(SAMPLEDIR_SILENCE, "silence_501ms.wav"))
_1  = read_wav(joinpath(SAMPLEDIR_SILENCE, "silence_1001ms.wav"))
_2  = read_wav(joinpath(SAMPLEDIR_SILENCE, "silence_2001ms.wav"))
_5  = read_wav(joinpath(SAMPLEDIR_SILENCE, "silence_5001ms.wav"))

twoword_cmds    = ["page up", "page down"]
singleword_cmds = [cmd for cmd in keys(commands) if cmd ∉ twoword_cmds]

@testset "$(basename(@__FILE__))" begin
    @testset "1. dynamic recognizers" begin
        valid_input = ["world", "universe"]
        @test isa(recognizer(valid_input, noises(DEFAULT_MODEL_NAME)), PyObject)
        @test isa(recognizer(valid_input, noises(TYPE_MODEL_NAME); modelname=TYPE_MODEL_NAME), PyObject)
    end;
    @testset "2. single-word recognitions ($cmd)" for cmd in singleword_cmds
        sample = samples[cmd]
        id     = cmd
        start_reading([sample; _2]; id=id)
        set_default_streamer(reader, id)
        @testset "partial" begin
            text, is_partial_result, has_timed_out = next_partial_recognition(recognizer(COMMAND_RECOGNIZER_ID))
            @test split(text)[end] == cmd
            @test is_partial_result
            @test !has_timed_out
        end;
        @testset "full" begin
            text = next_recognition(recognizer(COMMAND_RECOGNIZER_ID))
            @test split(text)[end] == cmd
        end;
        stop_reading(id=id)
    end;
    @testset "3. two-word recognitions ($cmd)" for cmd in twoword_cmds
        sample = samples[cmd]
        id     = cmd
        start_reading([sample; _2]; id=id)
        set_default_streamer(reader, id)
        @testset "full" begin
            next_partial_recognition(recognizer(COMMAND_RECOGNIZER_ID))
            text = next_recognition(recognizer(COMMAND_RECOGNIZER_ID))
            @test join(split(text)[end-1:end], " ") == cmd
        end;
        stop_reading(id=id)
    end;
    recognizer(COMMAND_RECOGNIZER_ID).Reset()
    reset_all(hard=true)
    @testset "4. token buffering - max speed ($cmd)" for cmd in singleword_cmds
        sample = samples[cmd]
        id     = cmd
        start_reading([sample; _05]; id=id)
        set_default_streamer(reader, id)
        @testset "_is_next" begin
            @test _is_next(cmd, recognizer(COMMAND_RECOGNIZER_ID), noises(DEFAULT_MODEL_NAME); use_partial_recognitions=true, ignore_unknown=false)
        end;
        # NOTE: The recogniser resets implied by dynamic recognisers in the following tests seem to break these tests.
        # @testset "is_next" begin
        #     @test is_next(cmd; use_max_speed=true, ignore_unknown=false)
        # end;
        # @testset "is_next (multi-token check)" begin
        #     @test is_next(["cat", cmd, "dog"]; use_max_speed=true, ignore_unknown=false)
        # end;
        @testset "next_token" begin
            token = next_token(recognizer(COMMAND_RECOGNIZER_ID), noises(DEFAULT_MODEL_NAME); use_partial_recognitions=true, ignore_unknown=false)
            @test token == cmd
        end;
        stop_reading(id=id)
    end;
    @testset "5. token buffering - max accuracy ($cmd)" for cmd in singleword_cmds
        sample = samples[cmd]
        id     = cmd
        start_reading([sample; _2]; id=id)
        set_default_streamer(reader, id)
        @testset "_is_next" begin
            @test _is_next(cmd, recognizer(COMMAND_RECOGNIZER_ID), noises(DEFAULT_MODEL_NAME); use_partial_recognitions=false, ignore_unknown=false)
        end;
        # NOTE: The recogniser resets implied by dynamic recognisers in the following tests seem to break these tests.
        # @testset "is_next" begin
        #     @test is_next(cmd; use_max_speed=false, ignore_unknown=false)
        # end;
        # @testset "is_next (multi-token check)" begin
        #     @test is_next(["cat", cmd, "dog"]; use_max_speed=false, ignore_unknown=false)
        # end;
        @testset "next_token" begin
            token = next_token(recognizer(COMMAND_RECOGNIZER_ID), noises(DEFAULT_MODEL_NAME); use_partial_recognitions=false, ignore_unknown=false)
            @test token == cmd
        end;
        stop_reading(id=id)
    end;
    @testset "6. token buffering (two-word retrieval) - max accuracy ($cmd)" for cmd in twoword_cmds
        sample = samples[cmd]
        id     = cmd
        start_reading([sample; _2]; id=id)
        set_default_streamer(reader, id)
        @testset "_are_next" begin
            tokens = String.(split(cmd))
            is_match, match = _are_next(tokens, recognizer(COMMAND_RECOGNIZER_ID), noises(DEFAULT_MODEL_NAME); use_partial_recognitions=false, ignore_unknown=false)
            @test is_match
            @test join(match, " ") == cmd
        end;
        @testset "are_next" begin
            tokens = String.(split(cmd))
            is_match, match = are_next(tokens; use_max_speed=false, ignore_unknown=false)
            @test is_match
            @test join(match, " ") == cmd
        end;
        @testset "are_next (additional tokens for check)" begin
            tokens = String.(split("cat "*cmd*" dog"))
            is_match, match = are_next(tokens; use_max_speed=false, ignore_unknown=false)
            @test is_match
            @test join(match, " ") == cmd
        end;
        @testset "next_token" begin
            token  = next_token(recognizer(COMMAND_RECOGNIZER_ID), noises(DEFAULT_MODEL_NAME); use_partial_recognitions=false, ignore_unknown=false)
            tokens = token * " " * next_token(recognizer(COMMAND_RECOGNIZER_ID), noises(DEFAULT_MODEL_NAME); use_partial_recognitions=false, ignore_unknown=false)
            @test tokens == cmd
        end;
        stop_reading(id=id)
    end;
end;

# Test tear down
finalize_jsi()


text:
export t_ind, invert_∆l

# Below, t_ind returns
# - a tuple if the first argument is a tuple of tuple, and
# - an SVec if the first argument is a tuple of vector.

# Below, earlier methods delegate actions to later methods.  Maybe they need to be
# implemented separately for speed?

# t_ind(t::Tuple3{T}, ind::Tuple3{Int}) where {T} = (t[ind[1]], t[ind[2]], t[ind[3]])
# t_ind(t::Tuple32{T}, i::Int, j::Int, k::Int) where {T} = (t[1][i], t[2][j], t[3][k])  # i, j, k = 1 or 2

# From a tuple of two tuples 1 and 2, each with K entries, construct a tuple (with K entries)
# whose kth entry is the kth entry of either tuple 1 or tuple 2.
# E.g., t_ind(((0.1,0.2,0.3), (1.0,2.0,3.0)), 1, 2, 1) = (0.1, 2.0, 0.3)
@inline t_ind(t::Tuple23, i₁₂::T, j₁₂::T, k₁₂::T) where {T<:Union{GridType,Sign,Integer}} = t_ind(t, (i₁₂,j₁₂,k₁₂))  # i₁₂, j₁₂, k₁₂ = 1 or 2
@inline t_ind(t::Tuple2{NTuple{K}}, ind₁₂::CartesianIndex{K}) where {K} = t_ind(t, ind₁₂.I)  # ind₁₂.I: NTuple{K,Int}
@inline t_ind(t::Tuple2{NTuple{K}}, ind₁₂::SVec{K,T}) where {K,T<:Union{GridType,Sign,Integer}} = t_ind(t, ind₁₂.data)
@inline t_ind(t::Tuple2{NTuple{K}}, ind₁₂::NTuple{K,T}) where {K,T<:Union{GridType,Sign,Integer}} =  # ind₁₂[k] = 1 or 2
    map((t₁,t₂,i) -> Int(i)==1 ? t₁ : t₂, t[1], t[2], ind₁₂)  # NTuple{K}

# From a tuple of two SVecs 1 and 2, each with K entries, construct an SVec (with K
# entries) whose kth entry is the kth entry of either SVec 1 or SVec 2.
# E.g., t_ind((SVec(0.1,0.2,0.3), SVec(1.0,2.0,3.0)), 1, 2, 1) = SVec(0.1, 2.0, 0.3)
@inline t_ind(t::Tuple2{SVec{3}}, i₁₂::T, j₁₂::T, k₁₂::T) where {T<:Union{GridType,Sign}} = t_ind(t, SVec(i₁₂,j₁₂,k₁₂))
# @inline t_ind(t::Tuple2{SVec{K}}, ind₁₂::CartesianIndex{K}) where {K} = t_ind(t, ind₁₂.I)
@inline t_ind(t::Tuple2{SVec{K}}, ind₁₂::NTuple{K,T}) where {K,T<:Union{GridType,Sign}} = t_ind(t, SVec(ind₁₂))
@inline t_ind(t::Tuple2{SVec{K}}, ind₁₂::SVec{K,T}) where {K,T<:Union{GridType,Sign}} =  # ind₁₂[k] = 1 or 2
    map((t₁,t₂,i) -> Int(i)==1 ? t₁ : t₂, t[1], t[2], ind₁₂)  # SVec{K}

# From a tuple of K vectors, construct a vector with K entries whose kth entry is
# taken from the kth vector.  Which entry to take from the kth vector is specified
# by indices.
# E.g., t_ind(([0.1,0.2,0.3,0.4], [1.0,2.0,3.0,4.0], [10.0,20.0,30.0,40.0]), 3, 1, 4) = SVec(0.3, 1.0, 40.0)
@inline t_ind(t::Tuple3{AbsVec}, i::Int, j::Int, k::Int) = t_ind(t, SVec(i,j,k))
@inline t_ind(t::Tuple2{AbsVec}, i::Int, j::Int) = t_ind(t, SVec(i,j))
@inline t_ind(t::NTuple{K,AbsVec}, ind::CartesianIndex{K}) where {K} = t_ind(t, ind.I)
@inline t_ind(t::NTuple{K,AbsVec}, ind::NTuple{K,Int}) where {K} = t_ind(t, SVec(ind))
@inline t_ind(t::NTuple{K,AbsVec}, ind::SVec{K,Int}) where {K} = map((tₖ,iₖ) -> tₖ[iₖ], SVec(t), ind)  # SVec{K}

function invert_∆l(∆l::Tuple2{Tuple{Vararg{AbsVecNumber}}})  # ∆l[PR][k]: grid point spacings at primal grid point locations on k-axis
    ∆lprim, ∆ldual = ∆l
    ∆lprim⁻¹ = map(v->(1 ./ v), ∆lprim)
    ∆ldual⁻¹ = map(v->(1 ./ v), ∆ldual)

    return (∆lprim⁻¹, ∆ldual⁻¹)
end

function movingavg(l::AbsVec{T}) where {T<:Number}
    # Return (l[1:end-1] + l[2:end]) / 2

    n = length(l)
    if n ≤ 1  # n = 0 or 1
        lmov = float(T)[]
    else  # n ≥ 2
        lmov = Vector{float(T)}(undef, n-1)
        for i = 1:n-1
            lmov[i] = (l[i] + l[i+1]) / 2
        end
    end

    return lmov
end

function movingmin(l::AbsVec{T}) where {T<:Number}
    # Return min.([l; Inf], [Inf; l])

    n = length(l)
    if n == 0
        lmov = [Inf]
    else  # n ≥ 1
        lmov = Vector{T}(undef, n+1)
        lmov[1] = l[1]
        lmov[n+1] = l[n]
        for i = 2:n
            lmov[i] = min(l[i-1], l[i])
        end
    end

    return lmov
end

# Numerical differentiation by forward difference.
f′fwd(f::Function, xₙ::Number, fₙ::Number=f(xₙ); h::Real=abs(xₙ)*Base.rtoldefault(Float)) = (f(xₙ+h) - fₙ) / h

function newtsol(x₀::Number, f::Function, f′::Function=(x,fₙ)->f′fwd(f,x,fₙ); rtol::Real=Base.rtoldefault(Float), atol::Real=eps(Float))
    # Solve f(x) = 0 using the Newton-Armijo method.
    # x₀: initial guess
    # f: function of x
    # f′: derivative of f at x.  If f′ can be evaluated using the value of f(x),
    #   write f′ such that it takes f(x) as the 2nd argument

    isconverged = true  # true if solution converged; false otherwise
    maxit = 100    # maximum iterations
    maxitls = 20   # maximum iterations inside the line search
    α = 1e-4
    maxs = 1 / Base.rtoldefault(Float)  # maximum step size
    perturbls = eps(Float)^0.75  # ≈ 1e-12, between sqrt(eps) and eps; some perturbation allowed in line search

    # Initialize.
    n = 0
    xₙ = float(x₀)
    T = typeof(xₙ)
    fₙ::T = f(xₙ)
    τf = rtol*abs(fₙ) + atol
    rx₀ = abs(x₀)
    has2ndarg = hasmethod(f′, Tuple2{Number})

    # Perform the Newton method.
    # @info "fₙ = $fₙ"
    while abs(fₙ) > τf
        # @info "n = $n"

        # abs(xₙ/x₀) ≤ 1e3 || (isconverged = false; break)
        # abs(xₙ/x₀) ≤ 1e3 || throw(ErrorException("Solution xₙ = $xₙ has diverged from x₀ = $x₀."))

        λ = 1.
        nls = 0  # line search iteration counter

        f′ₙ::T = has2ndarg ? f′(xₙ,fₙ) : f′(xₙ)

        # Avoid too large Newton steps.
        s = -fₙ/f′ₙ
        # @info "fₙ = $fₙ, f′ₙ = $f′ₙ, xₙ = $xₙ, s = $s"
        abs(s) ≤ maxs || (isconverged = false; break)
        # abs(s) ≤ maxs || throw(ErrorException("Newton step s = $s is larger than maximum step size $maxs."))
        xₙ₊₁ = xₙ + λ*s
        fₙ₊₁ = f(xₙ₊₁)

        # Perform the line search to determine λ.  The stopping criterion does not
        # have perturbls ≈ 1e-12 on the RHS, but I guess this kind of perturbation
        # allows update in xₙ even in the situation where line search is supposed
        # to fail.
        while abs(fₙ₊₁) ≥ (1 - α*λ) * abs(fₙ) + perturbls
            # @info "nls = $nls, fₙ₊₁ = $(fₙ₊₁)"
            λ /= 2
            xₙ₊₁ = xₙ + λ*s
            fₙ₊₁::T = f(xₙ₊₁)
            nls += 1

            # Too many iteration steps in line search
            nls ≤ maxitls || (isconverged = false; break)
            # nls ≤ maxitls || throw(ErrorException("Line search fails in $nls iteration steps."))
        end

        # Step accepted; continue the Newton method.
        xₙ = xₙ₊₁
        n += 1

        # Too many iteration steps in Newton's method.
        n ≤ maxit || (isconverged = false; break)
        # n ≤ maxit || throw(ErrorException("Newton method fails to converge in $n iteration steps."))

        fₙ = f(xₙ)
        # @info "fₙ = $fₙ"
    end
    # @info "Newton done"

    return xₙ, isconverged
end

# isapprox_kernel(x, y, rtol, atol) = norm2diff(x,y) ≤ atol + rtol * max(norm2(x), norm2(y))
# isapprox(x, y; rtol::Real=Base.rtoldefault(Float), atol::Real=eps(Float)) = isapprox_kernel(x, y, rtol, atol)

# # Compare x and y with respect to some large number L.  Useful when x or y is zero.
# # It is OK to use ≈ or ≉ for positive numbers (like ∆).
isapprox_wrt(x, y, L::Number, rtol::Real=Base.rtoldefault(Float)) = norm(x-y) ≤ rtol * abs(L)
# isapprox_wrt(x, y, L::Real, rtol::Real=Base.rtoldefault(Float)) = norm2diff(x,y) ≤ rtol * L

# # When this is enabled to support isapprox for arrays of tuples, gen_sublprim1d fails
# # at the line where curr[1:2] ≈ prev[end-1:end].  Don't understand why, because
# # the line uses isapprox for arrays of numbers.
# function isapprox(x::AbsArr{NTuple{N,T}}, y::AbsArr{NTuple{N,T}}; rtol::Real=Base.rtoldefault(Float), atol::Real=eps(Float)) where {N,T<:Number}
#     return isapprox_kernel(x, y, rtol, atol)
# end

function consolidate_similar!(l::AbsVecNumber, L::Number)
    # Consolidate approximately equal points.  Assume l is sorted.

    # Original code:
    # ind_unique = (≉).(@view(lprim0[1:end-1]), @view(lprim0[2:end]))  # length(ind_unique) = length(lprim0) - 1
    # push!(ind_unique, true)
    # lprim0 = lprim0[ind_unique]  # [lprim0[1:end-1][original ind_unique]; lprim0[end]]

    if !isempty(l)
        ind_sim = findall(isapprox_wrt.(@view(l[1:end-1]), @view(l[2:end]), L))

        # Below, l[end] must be always included as a unique entry.  There are two cases: where
        # l[end-1] is included as a unique entry or not.  If l[end-1] is included, it must be
        # very different from l[end], so we must include l[end].  If l[end-1] is not included,
        # it must be approximately l[end], but because l[end-1] is not included, we must include
        # l[end].  In other words, l[end] represents the last set of approximately equal entries
        # in any cases.
        deleteat!(l, ind_sim)
        @assert !isapprox_wrt(l[end], l[end-1], L)  # last entry is always unique
    end
end

function findsim(l1::AbsVecNumber, l2::AbsVecNumber, L::Number)
    ind = Int[]
    for i = 1:length(l1)
        v1 = l1[i]
        for v2 = l2
            isapprox_wrt(v1, v2, L) && push!(ind, i)
        end
    end

    return ind
end


text:
include("/home/astukalov/projects/oesc_jl/test/runtests.jl")

sets = [Set([:a]), Set([:b]), Set([:a :b]), Set([:a, :b, :c]), Set([:d, :e, :f]), Set([:a, :b, :c, :d, :e, :f])]
active = Set([:a, :b, :d])

include("/home/astukalov/projects/oesc_jl/oesc.jl")
params = OESC.CoverParams(0.1, 0.01, 0.01);
sm = OESC.SetMosaic(sets);
msm = OESC.mask(sm, Set([:a,:b,:e,:c]));
problem = OESC.CoverProblem(msm, params);
opt_model = OESC.opt_model(problem)

using JuMP
solve(opt_model)
getValue(getVar(opt_model, :w))

res = OESC.optimize(problem)

using Gadfly
plot(x=linspace(0.0,1.0,100), y=linspace(0.0,1.0,100),
     z=(x,y) -> OESC.score(problem, 1.0 - [x,y]),
     Geom.contour(levels=50))

OESC.score(problem, res.set_weights)

sm = OESC.SetMosaic(Set{Symbol}[Set(Symbol[:a])])
msm = OESC.mask(sm, Set{Symbol}([:a]))

en_problem = OESC.CoverProblem(msm, OESC.CoverParams(0.5, 0.5, 0.7));
plot(x -> OESC.score(en_problem, 1.0 - [x]), 0.0, 1.0)

dis_problem = OESC.CoverProblem(msm, OESC.CoverParams(0.5, 0.5, 0.4));
plot(x -> OESC.score(dis_problem, 1.0 - [x]), 0.0, 1.0)

en_res = OESC.optimize(en_problem)
@fact en_res.set_weights --> roughly(Float64[1.0], atol=1E-8)

sm = OESC.SetMosaic([Set(Symbol[:a, :b]), Set(Symbol[:b, :c]), Set(Symbol[:a, :b, :c])])
sm_ab = OESC.mask(sm, Set(Symbol[:a]))

problem_ab_lowp = OESC.CoverProblem(sm_ab, OESC.CoverParams(0.1, 0.1, 0.01))
@fact OESC.nsets(problem_ab_lowp) --> 3
@fact OESC.ntiles(problem_ab_lowp) --> 3

res_ab_lowp = OESC.optimize(problem_ab_lowp)
@fact res_ab_lowp.set_weights --> roughly(Float64[0.0, 0.0, 0.0], atol=1E-4)

problem_ab = OESC.CoverProblem(sm_ab, OESC.CoverParams(0.4, 0.1, 0.3))
res_ab = OESC.optimize(problem_ab)
@fact res_ab.set_weights --> roughly(Float64[1.0, 0.0, 0.0], atol=1E-4)

sm = OESC.SetMosaic([Set(Symbol[:a]), Set(Symbol[:b]),
                      Set(Symbol[:c]), Set(Symbol[:a, :b, :c])])
sm_ab = OESC.mask(sm, Set(Symbol[:a, :b]))

# low prior probability to select sets, high probability to miss active element, so select abc
problem_ab = OESC.CoverProblem(sm_ab, OESC.CoverParams(0.5, 0.1, 0.1))
OESC.score(problem_ab, [1.0, 1.0, 0.0])

using Gadfly
plot(x=linspace(0.0,1.0,100), y=linspace(0.0,1.0,100),
     z=(x,y) -> OESC.score(problem_ab, Float64[x,y,0.001]),
     Geom.contour(levels=50))

include("/home/astukalov/projects/oesc/oesc.jl")
sm = OESC.SetMosaic([Set(Symbol[:a]), Set(Symbol[:b]),
                      Set(Symbol[:c]), Set(Symbol[:a, :b, :c])])
sm_ab = OESC.mask(sm, Set(Symbol[:a, :b]))

cover_etor = OESC.CoverEnumerator(sm_ab, OESC.CoverParams(0.5, 0.1, 0.1))
cover_coll = collect(cover_etor)
DataTable(cover_coll)

cover_etor2 = OESC.CoverEnumerator(sm_ab, OESC.CoverParams(0.1, 0.1, 0.2))
cover_coll2 = collect(cover_etor2)

cover_results[(:MF, 1, :-, 3796)]

@where(cover_results_df, :prob .> 0.8)

ProfileView.view()

clu_go_bp = OESC.mask(go_bp, exp_comp_clusters["2+"][7])
OESC.nsets(clu_go_bp)
OESC.ntiles(clu_go_bp)

cover_etor = OESC.CoverEnumerator(clu_go_bp, OESC.CoverParams(0.5, 0.1, 0.1))
cover_coll = collect(cover_etor)