text:
"""
    TmuxDisplays

Using `tmux` as Julia's external displays.

API:

* `tmuxdisplay([object])`
* `TmuxDisplays.split_window(; ...)`
* `TmuxDisplays.new_window(; ...)`

They all return a `TmuxDisplays` (`<: AbstractDisplay`).
"""
module TmuxDisplays

export tmuxdisplay

using Base: Filesystem

struct TmuxDisplay <: AbstractDisplay
    pane_id::String
    pane_pid::Int
    waiter_fifo::Base.PipeEndpoint
    waiter_path::String
    tty::Base.TTY
end

Base.IO(tmux::TmuxDisplay) = tmux.tty

function Base.display(tmux::TmuxDisplay, mime::MIME"text/plain", x)
    io = IO(tmux)
    println(io)
    println(io)
    printstyled(io, "_"^displaysize(io)[2]; color = :blue)
    println(io)
    println(io, ' ')
    # print(io, "\e[3J")
    print(io, "\033c")
    show(io, mime, x)
    return
end

Base.display(tmux::TmuxDisplay, x) = display(tmux, MIME"text/plain"(), x)

Base.isopen(tmux::TmuxDisplay) = success(`tmux has-session -t $(tmux.pane_id)`)

function silent_run(cmd, msg)
    io = IOBuffer()
    success(pipeline(cmd; stdout=io, stderr=io)) ||
        @debug(
            "$cmd failed to run. ($msg)",
            output = Text(chomp(String(take!(io))))
        )
end

function Base.close(tmux::TmuxDisplay)
    silent_run(`tmux kill-pane -t $(tmux.pane_id)`, "already closed?")
    close(tmux.waiter_fifo)
    close(tmux.tty)
    cleanup_tmux(tmux)
end

function Base.show(io::IO, ::MIME"text/plain", tmux::TmuxDisplay)
    print(io, "TmuxDisplay: pane ", tmux.pane_id, ' ')
    if isopen(tmux)
        printstyled(io, "(open)"; color=:green)
    else
        printstyled(io, "(closed)"; color=:red)
    end
end

struct Lockable{T,L}
    value::T
    lock::L
end

Lockable(value) = Lockable(value, ReentrantLock())

Base.lock(f, lockable::Lockable) = lock(() -> f(lockable.value), lockable.lock)

const DISPLAYS = Lockable(Dict{String,TmuxDisplay}())
const DEFAULT_DISPLAY = Lockable(Ref{Union{Nothing,TmuxDisplay}}(nothing))

# Wrap C function?
mkfifo(path::AbstractString) = run(`mkfifo $path`)

function tmpfifo()
    path = tempname() * "_TmuxDisplays"
    mkfifo(path)
    # Retry when fail?
    return path
end

open_pipe(path::AbstractString) =
    Base.PipeEndpoint(fd(Filesystem.open(path, Filesystem.JL_O_RDONLY)))

make_tty(path::AbstractString) = Base.TTY(fd(Filesystem.open(path, Filesystem.JL_O_WRONLY)))

as_pane_id(tmux::TmuxDisplay) = tmux.pane_id
as_pane_id(id::AbstractString) = id
as_pane_id(x) = error("Cannot interpret as tmux pane:\n", x)

"""
    TmuxDisplays.split_window() :: TmuxDisplay

Use `tmux split-window` to create a new pane that is used as an
external display.

# Example

```julia
display(TmuxDisplays.split_window(), Text("hello"))
```

# Keyword Arguments
- `horizontal::Bool = false`: Split pane horizontally if `true`.
- `size::Integer`: Number of lines (for vertical split) or cells (for
  horizontal split).
- `target::Union{TmuxDisplay,AbstractString}`: Tmux pane to be split.
  If it is a string, it is passed to `-t` option of `split-window` as-is.
- `focus::Bool = false`: Focus newly created pane if `true`.
"""
function split_window(;
    horizontal = false,
    size = nothing,
    target = get(ENV, "TMUX_PANE", nothing),
    focus = false,
)
    cmd = `tmux split-window`
    if horizontal
        cmd = `$cmd -h`
    end
    if size !== nothing
        cmd = `$cmd -l $size`
    end
    if target !== nothing
        cmd = `$cmd -t $(as_pane_id(target))`
    end
    if !focus
        cmd = `$cmd -d`
    end

    waiter_path = tmpfifo()
    cmd = `
    $cmd
    -P
    -F '#{pane_id} #{pane_pid} #{pane_tty}'
    "sleep 2147483647d > $waiter_path"
    `
    return tmux_run(cmd, waiter_path)
end

"""
    TmuxDisplays.new_window(; [focus]) -> d::TmuxDisplay

Open a new window in the current tmux session.

# Example

```julia
display(TmuxDisplays.new_window(), Text("hello"))
```

# Keyword Arguments
- `focus::Bool = true`: Focus newly created pane if `true` (default).
- `name::String = "jl"`: Name of the window.
"""
function new_window(;
    # TODO: add support for other options
    focus = true,
    name = "jl",
)
    cmd = `tmux new-window`
    if !focus
        cmd = `$cmd -d`
    end

    waiter_path = tmpfifo()
    cmd = `
    $cmd
    -n $name
    -P
    -F '#{pane_id} #{pane_pid} #{pane_tty}'
    "sleep 2147483647d > $waiter_path"
    `
    return tmux_run(cmd, waiter_path)
end

function  tmux_run(cmd, waiter_path)
    pane_id, pane_pid, pane_tty = split(read(cmd, String))
    tmux = TmuxDisplay(
        pane_id,
        parse(Int, pane_pid),
        open_pipe(waiter_path),
        waiter_path,
        make_tty(pane_tty),
    )
    lock(DISPLAYS) do d
        d[tmux.pane_id] = tmux
    end
    @debug "Created: $(sprint(show, "text/plain", tmux))"

    # Cleanup `DISPLAYS`
    @async tmux_waiter(tmux)  # unstructured concurrency...

    return tmux
end

function tmux_waiter(tmux)
    try
        @debug "Waiting for $(tmux.waiter_path) to be closed."
        read(tmux.waiter_fifo)
    catch err
        @error("`read(tmux.waiter_fifo)` failed", exception = (err, catch_backtrace()))
    finally
        @debug "Closing: $(sprint(show, "text/plain", tmux))"
        try
            close(tmux)
        catch err
            @error("`close(tmux)` failed", exception = (err, catch_backtrace()))
        end
    end
end

function cleanup_tmux(tmux)
    try
        rm(tmux.waiter_path, force = true)
        lock(DISPLAYS) do d
            pop!(d, tmux.pane_id, nothing)
        end
    catch err
        @error("`cleanup_tmux` failed.", exception = (err, catch_backtrace()))
    end
end

"""
    tmuxdisplay() :: TmuxDisplay

Return the default `TmuxDisplay` instance.  Create one if it does not
exist.
"""
function tmuxdisplay()
    cleanupall()
    lock(DEFAULT_DISPLAY) do ref
        d = ref[]
        if d === nothing || !isopen(d)
            ref[] = d = split_window()
        end
        return d
    end
end

"""
    tmuxdisplay(x) :: TmuxDisplay

Display `x` in a separate pane of current tmux session.  Return the
`TmuxDisplay` used for showing `x`.
"""
function tmuxdisplay(x)
    tmux = tmuxdisplay()
    display(tmux, x)
    return tmux
end

"""
    TmuxDisplays.closeall()

Close all `TmuxDisplay`s.
"""
function closeall()
    lock(DEFAULT_DISPLAY) do ref
        ref[] = nothing
    end
    lock(DISPLAYS) do d
        while !isempty(d)
            (_, tmux) = pop!(d)
            close(tmux)
        end
    end
end

# TODO: get rid of this once I'm fine with `tmux_waiter`
function cleanupall()
    lock(DISPLAYS) do d
        for (pane_id, tmux) in collect(d)
            if !isopen(tmux)
                @error "Unclosed pane found" tmux
                close(tmux)
            end
        end
    end
end

function __init__()
    atexit(closeall)
end

end # module


text:
module DZOptimization

# export LineSearchFunctor, BFGSOptimizer, step!, find_saturation_threshold,
#     L2RegularizationWrapper, L2GradientWrapper

export LineSearchFunctor, L2RegularizationWrapper, L2GradientWrapper,
    step!, GradientDescentOptimizer, BFGSOptimizer, LBFGSOptimizer

using LinearAlgebra: mul!

########################################################## UNSAFE LINEAR ALGEBRA

@inline unsafe_sqrt(x::Float32) = Base.sqrt_llvm(x)
@inline unsafe_sqrt(x::Float64) = Base.sqrt_llvm(x)
@inline unsafe_sqrt(x::T) where {T} = sqrt(x)

@inline half(::Type{T}) where {T} = one(T) / (one(T) + one(T))

@inline NULL_CONSTRAINT(_...) = true

@inline function norm2(x::AbstractArray{T,N}) where {T,N}
    result = zero(real(T))
    @simd for i = 1 : length(x)
        @inbounds result += abs2(x[i])
    end
    return result
end

@inline function norm(x::AbstractArray{T,N}) where {T,N}
    return unsafe_sqrt(float(norm2(x)))
end

@inline function normalize!(x::AbstractArray{T,N}) where {T,N}
    a = inv(norm(x))
    @simd ivdep for i = 1 : length(x)
        @inbounds x[i] *= a
    end
    return x
end

function normalize_columns!(A::AbstractMatrix{T}) where {T}
    m, n = size(A)
    @inbounds for j = 1 : n
        norm_sq = zero(real(T))
        @simd ivdep for i = 1 : m
            norm_sq += abs2(A[i,j])
        end
        inv_norm = inv(unsafe_sqrt(norm_sq))
        @simd ivdep for i = 1 : m
            A[i,j] *= inv_norm
        end
    end
    return A
end

@inline function dot(v::AbstractArray{T,N},
                     w::AbstractArray{T,N}, n::Int) where {T,N}
    result = zero(T)
    @simd ivdep for i = 1 : n
        @inbounds result += conj(v[i]) * w[i]
    end
    return result
end

@inline function dot(v::AbstractArray{T,N},
                     w::AbstractMatrix{T}, n::Int, j::Int) where {T,N}
    result = zero(T)
    @simd ivdep for i = 1 : n
        @inbounds result += conj(v[i]) * w[i,j]
    end
    return result
end

@inline function identity_matrix!(A::AbstractMatrix{T}) where {T}
    m, n = size(A)
    for j = 1 : n
        @simd ivdep for i = 1 : m
            @inbounds A[i,j] = ifelse(i == j, one(T), zero(T))
        end
    end
    return A
end

@inline function negate!(v::AbstractArray{T,N},
                         w::AbstractArray{T,N}, n::Int) where {T,N}
    @simd ivdep for i = 1 : n
        @inbounds v[i] = -w[i]
    end
    return v
end

@inline function add!(v::AbstractArray{T,N},
                      w::AbstractArray{T,N}, n::Int) where {T,N}
    @simd ivdep for i = 1 : n
        @inbounds v[i] += w[i]
    end
    return v
end

@inline function add!(v::AbstractArray{T,N}, c::T,
                      w::AbstractArray{T,N}, n::Int) where {T,N}
    @simd ivdep for i = 1 : n
        @inbounds v[i] += c * w[i]
    end
    return v
end

@inline function add!(v::AbstractArray{T,N}, c::T,
                      w::AbstractMatrix{T}, n::Int, j::Int) where {T,N}
    @simd ivdep for i = 1 : n
        @inbounds v[i] += c * w[i,j]
    end
    return v
end

@inline function scalar_mul!(v::AbstractArray{T,N}, c::T) where {T,N}
    @simd ivdep for i = 1 : length(v)
        @inbounds v[i] *= c
    end
    return v
end

@inline linear_view(x::Array{T,N}) where {T,N} =
    reshape(view(x, ntuple(_ -> Colon(), N)...), length(x))

######################################################### LINE SEARCH ALGORITHMS

@inline function _qls_best(fb::T, x1::T, f1::T,
                           x2::T, f2::T, x3::T, f3::T) where {T}
    xb = zero(T)
    if f1 < fb; xb, fb = x1, f1; end
    if f2 < fb; xb, fb = x2, f2; end
    if f3 < fb; xb, fb = x3, f3; end
    return (xb, fb)
end

@inline function _qls_minimum_high(f0::T, f1::T, f2::T) where {T}
    q1 = f1 + f1
    q2 = q1 + q1
    q3 = f0 + f0
    q4 = f2 + f0
    q5 = q1 - q4
    return (q2 - q3 - q4) / (q5 + q5)
end

@inline function _qls_minimum_low(f0::T, f1::T, f2::T)::T where {T}
    q1 = f2 + f2
    q2 = q1 + q1
    q3 = f0 + f0
    q4 = f0 + f1
    q5 = q4 - q1
    q6 = q5 + q5
    return (q4 + q3 - q2) / (q6 + q6)
end

function quadratic_line_search(f::F, f0::T, x1::T,
                               max_increases::Int=10,
                               max_decreases::Int=100) where {F,T}
    # TODO: In principle, we could make this work for f0 == +Inf.
    if !isfinite(f0) || !isfinite(x1)
        return zero(T), f0
    end
    f1 = f(x1)
    while isnan(f1)
        x1 = half(T) * x1
        if iszero(x1)
            return zero(T), f0
        end
        f1 = f(x1)
    end
    if f1 < f0
        num_increases = 0
        while true
            x2 = x1 + x1
            f2 = f(x2)
            num_increases += 1
            if (f2 >= f1) || isnan(f2) || (num_increases > max_increases)
                x3 = x1 * _qls_minimum_high(f0, f1, f2)
                return _qls_best(f0, x1, f1, x2, f2, x3, f(x3))
            else
                x1, f1 = x2, f2
            end
        end
    else
        num_decreases = 0
        while true
            x2 = half(T) * x1
            f2 = f(x2)
            num_decreases += 1
            if isnan(f2)
                return (zero(T), f0)
            end
            if (f2 <= f0) || (num_decreases > max_decreases)
                x3 = x1 * _qls_minimum_low(f0, f1, f2)
                return _qls_best(f0, x2, f2, x1, f1, x3, f(x3))
            else
                x1, f1 = x2, f2
            end
        end
    end
end

############################################################ LINE SEARCH FUNCTOR

struct LineSearchFunctor{F,C,T,N}
    objective_function::F
    constraint_function!::C
    current_point::Array{T,N}
    new_point::Array{T,N}
    step_direction::Array{T,N}
end

@inline function (lsf::LineSearchFunctor{F,C,T,N})(
                  step_size::T) where {F,C,T,N}
    x0, x1, dx = lsf.current_point, lsf.new_point, lsf.step_direction
    @simd ivdep for i = 1 : length(x0)
        @inbounds x1[i] = x0[i] - step_size * dx[i]
    end
    lsf.constraint_function!(x1) && return lsf.objective_function(x1)
    return typemax(T)
end

######################################################## REGULARIZATION WRAPPERS

struct L2RegularizationWrapper{F,T}
    wrapped_function::F
    lambda::T
end

struct L2GradientWrapper{G,T}
    wrapped_gradient!::G
    lambda::T
end

function (wrapper::L2RegularizationWrapper{F,T})(
          x::AbstractArray{T,N}) where {F,T,N}
    return wrapper.wrapped_function(x) + half(T) * wrapper.lambda * norm2(x)
end

function (wrapper::L2GradientWrapper{G,T})(
          g::AbstractArray{T,N}, x::AbstractArray{T,N}) where {G,T,N}
    wrapper.wrapped_gradient!(g, x)
    lambda = wrapper.lambda
    @simd ivdep for i = 1 : length(g)
        @inbounds g[i] += lambda * x[i]
    end
    return g
end

############################################################### GRADIENT DESCENT

struct GradientDescentOptimizer{F,G,C,T,N}
    objective_function::F
    gradient_function!::G
    constraint_function!::C
    iteration_count::Array{Int,0}
    has_converged::Array{Bool,0}
    current_point::Array{T,N}
    current_objective_value::Array{T,0}
    current_gradient::Array{T,N}
    delta_point::Array{T,N}
    delta_gradient::Array{T,N}
    last_step_size::Array{T,0}
    _line_search_functor::LineSearchFunctor{F,C,T,N}
end

function GradientDescentOptimizer(objective_function::F,
                                  gradient_function!::G,
                                  initial_point::Array{T,N},
                                  initial_step_size::T) where {F,G,T,N}
    return GradientDescentOptimizer(
        objective_function, gradient_function!,
        NULL_CONSTRAINT, initial_point, initial_step_size)
end

function GradientDescentOptimizer(objective_function::F,
                                  gradient_function!::G,
                                  constraint_function!::C,
                                  initial_point::Array{T,N},
                                  initial_step_size::T) where {F,G,C,T,N}
    iteration_count = fill(0)
    has_converged = fill(false)
    current_point = copy(initial_point)
    constraint_success = constraint_function!(current_point)
    @assert constraint_success
    initial_objective_value = objective_function(current_point)
    @assert !isnan(initial_objective_value)
    current_objective_value = fill(initial_objective_value)
    current_gradient = similar(initial_point)
    gradient_function!(current_gradient, current_point)
    delta_point = zero(initial_point)
    delta_gradient = zero(initial_point)
    last_step_size = fill(initial_step_size)
    _scratch_space = similar(initial_point)
    _line_search_functor = LineSearchFunctor{F,C,T,N}(
        objective_function, constraint_function!,
        current_point, _scratch_space, current_gradient)
    return GradientDescentOptimizer{F,G,C,T,N}(
        objective_function,
        gradient_function!,
        constraint_function!,
        iteration_count,
        has_converged,
        current_point,
        current_objective_value,
        current_gradient,
        delta_point,
        delta_gradient,
        last_step_size,
        _line_search_functor)
end

function step!(opt::GradientDescentOptimizer{F,G,C,T,N}) where {F,G,C,T,N}

    if !opt.has_converged[]

        point = opt.current_point
        delta_point = opt.delta_point
        gradient = opt.current_gradient
        delta_gradient = opt.delta_gradient

        n = length(point)
        @assert n == length(delta_point)
        @assert n == length(gradient)
        @assert n == length(delta_gradient)

        step_size = opt.last_step_size[]
        gradient_norm = norm(gradient)
        next_step_size, next_obj = quadratic_line_search(
            opt._line_search_functor,
            opt.current_objective_value[],
            step_size / gradient_norm)

        if next_obj < opt.current_objective_value[]

            # Accept gradient descent step.
            opt.current_objective_value[] = next_obj
            opt.last_step_size[] = next_step_size * gradient_norm
            opt.iteration_count[] += 1

            # Update point and gradient.
            negate!(delta_point, point, n)
            negate!(delta_gradient, gradient, n)
            add!(point, -next_step_size, gradient, n)
            constraint_success = opt.constraint_function!(point)
            @assert constraint_success
            opt.gradient_function!(gradient, point)
            add!(delta_point, point, n)
            add!(delta_gradient, gradient, n)

        else
            opt.has_converged[] = true
        end

    end
    return opt
end

########################################################################### BFGS

@enum StepType begin
    NullStep
    GradientDescentStep
    BFGSStep
end

struct BFGSOptimizer{F,G,C,T,N}
    objective_function::F
    gradient_function!::G
    constraint_function!::C
    iteration_count::Array{Int,0}
    has_converged::Array{Bool,0}
    current_point::Array{T,N}
    current_objective_value::Array{T,0}
    current_gradient::Array{T,N}
    delta_point::Array{T,N}
    delta_gradient::Array{T,N}
    last_step_size::Array{T,0}
    last_step_type::Array{StepType,0}
    approximate_inverse_hessian::Matrix{T}
    next_step_direction::Array{T,N}
    _scratch_space::Array{T,N}
    _gradient_line_search_functor::LineSearchFunctor{F,C,T,N}
    _bfgs_line_search_functor::LineSearchFunctor{F,C,T,N}
end

function BFGSOptimizer(objective_function::F,
                       gradient_function!::G,
                       initial_point::Array{T,N},
                       initial_step_size::T) where {F,G,T,N}
    return BFGSOptimizer(
        objective_function, gradient_function!,
        NULL_CONSTRAINT, initial_point, initial_step_size)
end

function BFGSOptimizer(objective_function::F,
                       gradient_function!::G,
                       constraint_function!::C,
                       initial_point::Array{T,N},
                       initial_step_size::T) where {F,G,C,T,N}
    iteration_count = fill(0)
    has_converged = fill(false)
    current_point = copy(initial_point)
    constraint_success = constraint_function!(current_point)
    @assert constraint_success
    initial_objective_value = objective_function(current_point)
    @assert !isnan(initial_objective_value)
    current_objective_value = fill(initial_objective_value)
    current_gradient = similar(initial_point)
    gradient_function!(current_gradient, current_point)
    delta_point = zero(initial_point)
    delta_gradient = zero(initial_point)
    last_step_size = fill(initial_step_size)
    last_step_type = fill(NullStep)
    approximate_inverse_hessian = Matrix{T}(undef,
        length(initial_point), length(initial_point))
    identity_matrix!(approximate_inverse_hessian)
    next_step_direction = copy(current_gradient)
    _scratch_space = similar(initial_point)
    _gradient_line_search_functor = LineSearchFunctor{F,C,T,N}(
        objective_function, constraint_function!,
        current_point, _scratch_space, current_gradient)
    _bfgs_line_search_functor = LineSearchFunctor{F,C,T,N}(
        objective_function, constraint_function!,
        current_point, _scratch_space, next_step_direction)
    return BFGSOptimizer{F,G,C,T,N}(
        objective_function,
        gradient_function!,
        constraint_function!,
        iteration_count,
        has_converged,
        current_point,
        current_objective_value,
        current_gradient,
        delta_point,
        delta_gradient,
        last_step_size,
        last_step_type,
        approximate_inverse_hessian,
        next_step_direction,
        _scratch_space,
        _gradient_line_search_functor,
        _bfgs_line_search_functor)
end

function BFGSOptimizer(::Type{T},
                       opt::BFGSOptimizer{F,G,C,U,N}) where {F,G,C,U,T,N}
    return BFGSOptimizer(T, opt.objective_function, opt.gradient_function!,
                         opt.constraint_function!, opt)
end

function BFGSOptimizer(::Type{T},
                       objective_function::F,
                       gradient_function!::G,
                       constraint_function!::C,
                       opt::BFGSOptimizer{F2,G2,C2,T2,N}
                       ) where {F,G,C,T,F2,G2,C2,T2,N}
    current_point = T.(opt.current_point)
    constraint_success = constraint_function!(current_point)
    @assert constraint_success
    initial_objective_value = objective_function(current_point)
    @assert !isnan(initial_objective_value)
    current_gradient = similar(current_point)
    gradient_function!(current_gradient, current_point)
    approximate_inverse_hessian = T.(opt.approximate_inverse_hessian)
    next_step_direction = similar(current_point)
    mul!(linear_view(next_step_direction),
         approximate_inverse_hessian,
         linear_view(current_gradient))
    _scratch_space = similar(current_point)
    _gradient_line_search_functor = LineSearchFunctor{F,C,T,N}(
        objective_function, constraint_function!,
        current_point, _scratch_space, current_gradient)
    _bfgs_line_search_functor = LineSearchFunctor{F,C,T,N}(
        objective_function, constraint_function!,
        current_point, _scratch_space, next_step_direction)
    return BFGSOptimizer{F,G,C,T,N}(
        objective_function,
        gradient_function!,
        constraint_function!,
        fill(opt.iteration_count[]),
        fill(false),
        current_point,
        fill(initial_objective_value),
        current_gradient,
        T.(opt.delta_point),
        T.(opt.delta_gradient),
        fill(opt.last_step_size[]),
        fill(opt.last_step_type[]),
        approximate_inverse_hessian,
        next_step_direction,
        _scratch_space,
        _gradient_line_search_functor,
        _bfgs_line_search_functor)
end

function update_inverse_hessian!(
        inv_hess::Matrix{T}, step_size::T, step_direction::Array{T,N},
        delta_gradient::Array{T,N}, scratch_space::Array{T,N}) where {T,N}

    n = length(step_direction)
    @assert n == length(delta_gradient)
    @assert n == length(scratch_space)
    @assert (n, n) == size(inv_hess)

    overlap = dot(step_direction, delta_gradient, n)
    scalar_mul!(step_direction, inv(overlap))
    mul!(linear_view(scratch_space), inv_hess, linear_view(delta_gradient))
    delta_norm = step_size * overlap + dot(delta_gradient, scratch_space, n)

    @inbounds for j = 1 : n
        sj = step_direction[j]
        tj = scratch_space[j]
        @simd ivdep for i = 1 : n
            inv_hess[i, j] += (
                delta_norm * (step_direction[i] * sj)
                - (scratch_space[i] * sj + step_direction[i] * tj))
        end
    end

    return inv_hess
end

function step!(opt::BFGSOptimizer{F,G,C,T,N}) where {F,G,C,T,N}

    if !opt.has_converged[]

        point = opt.current_point
        delta_point = opt.delta_point
        gradient = opt.current_gradient
        delta_gradient = opt.delta_gradient
        bfgs_direction = opt.next_step_direction

        n = length(point)
        @assert n == length(delta_point)
        @assert n == length(gradient)
        @assert n == length(delta_gradient)
        @assert n == length(bfgs_direction)

        # We actually run two independent line searches on each iteration:
        # one in the quasi-Newton (Hessian)^-1 * (-gradient) direction, and one
        # in the raw (-gradient) direction. The idea is that the BFGS algorithm
        # attempts to incrementally build an approximation to the inverse
        # Hessian. However, if we spend a long time in one region of search
        # space, then quickly move to a different region, our previous Hessian
        # approximation may become unhelpful. By performing this "competitive
        # line search" procedure, we can detect whether this occurs and reset
        # the Hessian when necessary.

        # Use the previous step size as the initial guess for line search.
        step_size = opt.last_step_size[]

        # Launch line search in raw (-gradient) direction.
        grad_norm = norm(gradient)
        grad_step_size, grad_obj = quadratic_line_search(
            opt._gradient_line_search_functor,
            opt.current_objective_value[],
            step_size / grad_norm)

        # Launch line search in BFGS (Hessian)^-1 * (-gradient) direction.
        bfgs_norm = norm(bfgs_direction)
        bfgs_step_size, bfgs_obj = quadratic_line_search(
            opt._bfgs_line_search_functor,
            opt.current_objective_value[],
            step_size / bfgs_norm)

        if bfgs_obj < opt.current_objective_value[] && !(bfgs_obj > grad_obj)

            # Accept BFGS step.
            opt.current_objective_value[] = bfgs_obj
            opt.last_step_size[] = bfgs_step_size * bfgs_norm
            opt.last_step_type[] = BFGSStep
            opt.iteration_count[] += 1

            # Update point and gradient.
            negate!(delta_point, point, n)
            negate!(delta_gradient, gradient, n)
            add!(point, -bfgs_step_size, bfgs_direction, n)
            constraint_success = opt.constraint_function!(point)
            @assert constraint_success
            opt.gradient_function!(gradient, point)
            add!(delta_point, point, n)
            add!(delta_gradient, gradient, n)

            # Update inverse Hessian approximation using delta_gradient.
            update_inverse_hessian!(opt.approximate_inverse_hessian,
                                    -bfgs_step_size, bfgs_direction,
                                    delta_gradient, opt._scratch_space)

            # Compute next step direction using approximate inverse Hessian.
            mul!(linear_view(bfgs_direction),
                 opt.approximate_inverse_hessian,
                 linear_view(gradient))

        elseif grad_obj < opt.current_objective_value[]

            # Accept gradient descent step.
            opt.current_objective_value[] = grad_obj
            opt.last_step_size[] = grad_step_size * grad_norm
            opt.last_step_type[] = GradientDescentStep
            opt.iteration_count[] += 1

            # Update point and gradient.
            negate!(delta_point, point, n)
            negate!(delta_gradient, gradient, n)
            add!(point, -grad_step_size, gradient, n)
            constraint_success = opt.constraint_function!(point)
            @assert constraint_success
            opt.gradient_function!(gradient, point)
            add!(delta_point, point, n)
            add!(delta_gradient, gradient, n)

            # Reset approximate inverse Hessian to the identity matrix.
            identity_matrix!(opt.approximate_inverse_hessian)

            # Reset next step direction to gradient.
            @simd ivdep for i = 1 : n
                @inbounds bfgs_direction[i] = gradient[i]
            end

        else
            opt.has_converged[] = true
        end

    end
    return opt
end

######################################################################### L-BGFS

struct LBFGSOptimizer{F,G,C,T,N}
    objective_function::F
    gradient_function!::G
    constraint_function!::C
    iteration_count::Array{Int,0}
    has_converged::Array{Bool,0}
    current_point::Array{T,N}
    current_objective_value::Array{T,0}
    current_gradient::Array{T,N}
    delta_point::Array{T,N}
    delta_gradient::Array{T,N}
    last_step_size::Array{T,0}
    next_step_direction::Array{T,N}
    _alpha_history::Vector{T}
    _rho_history::Vector{T}
    _delta_point_history::Matrix{T}
    _delta_gradient_history::Matrix{T}
    _line_search_functor::LineSearchFunctor{F,C,T,N}
end

function LBFGSOptimizer(objective_function::F,
                        gradient_function!::G,
                        initial_point::Array{T,N},
                        initial_step_size::T,
                        history_length::Int=10) where {F,G,T,N}
    return LBFGSOptimizer(
        objective_function, gradient_function!,
        NULL_CONSTRAINT, initial_point, initial_step_size, history_length)
end

function LBFGSOptimizer(objective_function::F,
                        gradient_function!::G,
                        constraint_function!::C,
                        initial_point::AbstractArray{T,N},
                        initial_step_size::T,
                        history_length::Int=10) where {F,G,C,T,N}
    @assert history_length > 0
    iteration_count = fill(0)
    has_converged = fill(false)
    current_point = copy(initial_point)
    constraint_success = constraint_function!(current_point)
    @assert constraint_success
    initial_objective_value = objective_function(current_point)
    @assert !isnan(initial_objective_value)
    current_objective_value = fill(initial_objective_value)
    current_gradient = similar(initial_point)
    gradient_function!(current_gradient, current_point)
    delta_point = zero(initial_point)
    delta_gradient = zero(initial_point)
    last_step_size = fill(initial_step_size)
    next_step_direction = scalar_mul!(copy(current_gradient),
        initial_step_size / norm(current_gradient))
    _alpha_history = zeros(T, history_length)
    _rho_history = zeros(T, history_length)
    _delta_point_history = zeros(T, length(initial_point), history_length)
    _delta_gradient_history = zeros(T, length(initial_point), history_length)
    _line_search_functor = LineSearchFunctor{F,C,T,N}(
        objective_function, constraint_function!,
        current_point, similar(initial_point), next_step_direction)
    return LBFGSOptimizer{F,G,C,T,N}(
        objective_function,
        gradient_function!,
        constraint_function!,
        iteration_count,
        has_converged,
        current_point,
        current_objective_value,
        current_gradient,
        delta_point,
        delta_gradient,
        last_step_size,
        next_step_direction,
        _alpha_history,
        _rho_history,
        _delta_point_history,
        _delta_gradient_history,
        _line_search_functor)
end

function step!(opt::LBFGSOptimizer{S1,S2,S3,T,N}) where {S1,S2,S3,T,N}

    if !opt.has_converged[]

        point = opt.current_point
        delta_point = opt.delta_point
        gradient = opt.current_gradient
        delta_gradient = opt.delta_gradient
        next_step_direction = opt.next_step_direction

        n = length(point)
        @assert n == length(delta_point)
        @assert n == length(gradient)
        @assert n == length(delta_gradient)
        @assert n == length(next_step_direction)

        alpha = opt._alpha_history
        rho = opt._rho_history
        delta_point_history = opt._delta_point_history
        delta_gradient_history = opt._delta_gradient_history

        m = length(alpha)
        @assert m == length(rho)
        @assert (n, m) == size(delta_point_history)
        @assert (n, m) == size(delta_gradient_history)

        next_step_size, next_obj = quadratic_line_search(
            opt._line_search_functor,
            opt.current_objective_value[],
            one(T))

        if next_obj < opt.current_objective_value[]

            # Accept L-BFGS step.
            opt.current_objective_value[] = next_obj
            opt.last_step_size[] = next_step_size * norm(next_step_direction)
            opt.iteration_count[] += 1

            # Update point and gradient.
            negate!(delta_point, point, n)
            negate!(delta_gradient, gradient, n)
            add!(point, -next_step_size, next_step_direction, n)
            constraint_success = opt.constraint_function!(point)
            @assert constraint_success
            opt.gradient_function!(gradient, point)
            add!(delta_point, point, n)
            add!(delta_gradient, gradient, n)

            # Store delta_point and delta_gradient in history.
            c = Base.srem_int(opt.iteration_count[] - 1, m) + 1 # cyclic index
            @simd ivdep for i = 1 : n
                @inbounds delta_point_history[i,c] = delta_point[i]
            end
            @simd ivdep for i = 1 : n
                @inbounds delta_gradient_history[i,c] = delta_gradient[i]
            end
            delta_overlap = dot(delta_point, delta_gradient, n)
            @inbounds rho[c] = inv(delta_overlap)

            # Compute next step direction, starting from current gradient.
            @simd ivdep for i = 1 : n
                @inbounds next_step_direction[i] = gradient[i]
            end

            # Apply forward L-BFGS correction.
            history_count = max(opt.iteration_count[] - m + 1, 1)
            for iter = opt.iteration_count[] : -1 : history_count;
                c = Base.srem_int(iter - 1, m) + 1
                @inbounds overlap = rho[c] * dot(
                    next_step_direction, delta_point_history, n, c)
                @inbounds alpha[c] = overlap
                add!(next_step_direction,
                     overlap, delta_gradient_history, n, c)
            end

            # Compute natural step size.
            gamma = delta_overlap / norm2(delta_gradient)
            if !isfinite(gamma)
                gamma = sqrt(eps(T)) / max(one(T), norm2(gradient))
            end
            scalar_mul!(next_step_direction, gamma)

            # Apply backward L-BFGS correction.
            for iter = history_count : opt.iteration_count[]
                c = Base.srem_int(iter - 1, m) + 1
                @inbounds overlap = alpha[c] - rho[c] * dot(
                    next_step_direction, delta_gradient_history, n, c)
                add!(next_step_direction,
                     overlap, delta_point_history, n, c)
            end

        else
            opt.has_converged[] = true
        end

    end
    return opt
end

############################################################## OPTIMIZER TESTING

function run_and_test!(opt)

    history = [deepcopy(opt)]
    while !opt.has_converged[]
        step!(opt)
        push!(history, deepcopy(opt))
    end

    # None of the optimizers should have converged, except the last one.
    for i = 1 : length(history) - 1
        @assert !history[i].has_converged[]
    end
    @assert history[end].has_converged[]

    # Verify consistency of opt.iteration_count.
    for i = 1 : length(history) - 1
        @assert history[i].iteration_count[] == i - 1
    end
    history[end-1].iteration_count[] == history[end].iteration_count[]

    # Verify consistency of opt.current_objective_value.
    for opt in history
        @assert opt.objective_function(opt.current_point) ==
                opt.current_objective_value[]
    end

    # Verify consistency of opt.current_gradient.
    for opt in history
        grad = similar(opt.current_gradient)
        opt.gradient_function!(grad, opt.current_point)
        if grad != opt.current_gradient
            println("ERROR: ", norm(grad - opt.current_gradient))
        end
        @assert grad == opt.current_gradient
    end

    # Verify consistency of opt.delta_point.
    for i = 1 : length(history) - 2
        delta_point = history[i+1].current_point - history[i].current_point
        @assert delta_point == history[i+1].delta_point
    end
    @assert history[end-1].current_point == history[end].current_point

    # Verify consistency of opt.delta_gradient.
    for i = 1 : length(history) - 2
        delta_grad = history[i+1].current_gradient - history[i].current_gradient
        @assert delta_grad == history[i+1].delta_gradient
    end
    @assert history[end-1].current_gradient == history[end].current_gradient

    return opt
end

################################################################## MODEL FITTING

function find_saturation_threshold(data::Vector{Tuple{T,T}}) where {T}

    n = length(data)

    @inbounds sx, sy = data[1]
    sxx = sx * sx
    sxy = sx * sy

    y_mean = zero(T)
    @inbounds y_norm = sy^2
    @simd ivdep for i = 2 : n
        @inbounds y_mean += data[i][2]
        @inbounds y_norm += data[i][2]^2
    end
    y_mean /= (n - 1)

    best_loss = typemax(T)
    best_threshold = zero(T)
    best_slope = zero(T)
    best_mean = zero(T)

    for i = 2 : n-1

        @inbounds xi, yi = data[i]
        sx += xi
        sy += yi
        sxx += xi * xi
        sxy += xi * yi

        y_mean = (y_mean * (n - i + 1) - yi) / (n - i)

        @inbounds lower = data[i][1]
        threshold = lower
        begin
            slope = (i * y_mean - sy) / (i * threshold - sx)
            loss = half(T) * (
                y_norm
                + slope^2 * sxx
                - i * (y_mean - slope * threshold)^2
                - (n - i) * y_mean^2
            ) - slope * sxy
            if loss < best_loss
                best_loss = loss
                best_threshold = threshold
                best_slope = slope
                best_mean = y_mean
            end
        end

        @inbounds upper = data[i+1][1]
        threshold = upper
        begin
            slope = (i * y_mean - sy) / (i * threshold - sx)
            loss = half(T) * (
                y_norm
                + slope^2 * sxx
                - i * (y_mean - slope * threshold)^2
                - (n - i) * y_mean^2
            ) - slope * sxy
            if loss < best_loss
                best_loss = loss
                best_threshold = threshold
                best_slope = slope
                best_mean = y_mean
            end
        end

        a = i * sxx - sx * sx
        b = i * sxy - sx * sy
        c = sx * sxy - sxx * sy
        threshold = (c / b) + (a / b) * y_mean
        if lower <= threshold <= upper
            slope = (i * y_mean - sy) / (i * threshold - sx)
            loss = half(T) * (
                y_norm
                + slope^2 * sxx
                - i * (y_mean - slope * threshold)^2
                - (n - i) * y_mean^2
            ) - slope * sxy
            if loss < best_loss
                best_loss = loss
                best_threshold = threshold
                best_slope = slope
                best_mean = y_mean
            end
        end

    end
    return (best_threshold, best_slope, best_mean, best_loss)
end

####################################################################### INCLUDES

include("./ExampleFunctions.jl")

end # module DZOptimization

ext: jl


text:
using JellyFish
using Base.Test

jw = jaro_winkler("jellyfish", "smellyfish")
jw_correct = 0.8962962962962964
@test abs(jw_correct - jw) < 0.0001

jd = jaro_distance("jellyfish", "smellyfish")
jd_correct = 0.8962962962962964
@test abs(jd_correct - jd) < 0.0001

hd = hamming_distance("apple", "pear")
hd_correct = 5
@test abs(hd_correct - hd) < 0.0001

lev = levenshtein_distance("jellyfish", "smellyfish")
lev_correct = 2
@test abs(lev_correct - lev) < 0.0001

dlev = damerau_levenshtein_distance("jellyfish", "jellyfihs")
dlev_correct = 1
@test abs(dlev_correct - dlev) < 0.0001

sx = soundex("jellyfish")
sx_correct = "J412"
@test sx == sx_correct

met = metaphone("jellyfish")
met_correct = "JLFX"
@test met == met_correct

mat = match_rating_codex("jellyfish")
mat_correct = "JLLFSH"
@test mat == mat_correct

mra = match_rating_comparison("jellyfish", "jellyfihs")
mra_correct = true
@test mra == mra_correct

nys = nysiis("jellyfish")
nys_correct = "JALYF"
@test nys == nys_correct


text:
################################################# FILE DESCRIPTION #########################################################

# This file contains tests for graph conversions.

############################################################################################################################

@testset "Eachvertex" begin
   g = completegraph(10)

   # Set vertex properties
   setvprop!(g, :, 1 : 10, :p1)
   setvprop!(g, :, 1 : 10, :p2)
   setvprop!(g, :, 1 : 10, :p3)
   setvprop!(g, :, 1 : 10, :p4)

   arr = getvprop(g, :, :p1)

   # Straightforward fetch prop
   @test @query(g |> eachvertex(v.p1)) == arr

   # Add constant
   @test @query(g |> eachvertex(v.p2 + 3)) == arr .+ 3

   # Add two props
   @test @query(g |> eachvertex(v.p2 + v.p3)) == arr .+ arr

   # Parenthesis
   @test @query(g |> eachvertex((v.p1 + v.p2)/(v.p3 * v.p4))) == (arr .+ arr) ./ (arr .* arr)

   # Special tokens
   setlabel!(g, map(string, 1 : 10))
   @test @query(g |> eachvertex(v.id)) == collect(vertices(g))
   @test @query(g |> eachvertex(v.adj)) == [fadj(g, v) for v in vertices(g)]
   @test @query(g |> eachvertex(v.nbors)) == [g[v] for v in encode(g)]
   @test @query(g |> eachvertex(v.indegree)) == indegree(g)
   @test @query(g |> eachvertex(v.outdegree)) == outdegree(g)
   @test @query(g |> eachvertex(v.label)) == encode(g)
end

@testset "EachEdge" begin
   g = completegraph(10)

   # Set vertex properties
   setvprop!(g, :, 1 : 10, :p1)
   setvprop!(g, :, 1 : 10, :p2)
   setvprop!(g, :, 1 : 10, :p3)
   setvprop!(g, :, 1 : 10, :p4)

   # Set edge properties
   seteprop!(g, :, 1 : 90, :p1)
   seteprop!(g, :, 1 : 90, :p2)
   seteprop!(g, :, 1 : 90, :p3)
   seteprop!(g, :, 1 : 90, :p4)

   arr = geteprop(g, :, :p1)

   # Straightforward fetch prop
   @test @query(g |> eachedge(e.p1)) == arr

   # Add constant
   @test @query(g |> eachedge(e.p2 + 3)) == arr .+ 3

   # Add two props
   @test @query(g |> eachedge(e.p2 + e.p3)) == arr .+ arr

   # Parenthesis
   @test @query(g |> eachedge((e.p1 + e.p2)/(e.p3 * e.p4))) == (arr .+ arr) ./ (arr .* arr)

   # Source == target
   @test @query(g |> eachedge(s.p1 == t.p1)) == falses(90)

   # Source + target
   @test @query(g |> eachedge(s.p1 + t.p2 > 0 )) == trues(90)

   # Special tokens
   setlabel!(g, map(string, 1 : 10))
   eit = edges(g)
   @test @query(g |> eachedge(s.id)) == eit.us
   @test @query(g |> eachedge(t.id)) == eit.vs

   @test @query(g |> eachedge(s.label)) == encode(g, eit.us)
   @test @query(g |> eachedge(t.label)) == encode(g, eit.vs)

   @test @query(g |> eachedge(e.source)) == @query(g |> eachedge(s.label)) == encode(g, eit.us)
   @test @query(g |> eachedge(e.target)) == @query(g |> eachedge(t.label)) == encode(g, eit.vs)

   @test @query(g |> eachedge(s.adj)) == [fadj(g, v) for v in eit.us]
   @test @query(g |> eachedge(t.adj)) == [fadj(g, v) for v in eit.vs]

   @test @query(g |> eachedge(s.nbors)) == [g[v] for v in encode(g, eit.us)]
   @test @query(g |> eachedge(t.nbors)) == [g[v] for v in encode(g, eit.vs)]

   @test @query(g |> eachedge(s.indegree)) == indegree(g, eit.us)
   @test @query(g |> eachedge(s.outdegree)) == outdegree(g, eit.us)

   @test @query(g |> eachedge(t.indegree)) == indegree(g, eit.vs)
   @test @query(g |> eachedge(t.outdegree)) == outdegree(g, eit.vs)
end

@testset "Filter" begin
   g = completegraph(10)
   setlabel!(g, collect(1:10))
   # Set vertex properties
   setvprop!(g, :, 1 : 10, :p1)
   setvprop!(g, :, 1 : 10, :p2)
   setvprop!(g, :, 1 : 10, :p3)
   setvprop!(g, :, 1 : 10, :p4)

   # Set edge properties
   seteprop!(g, :, 1 : 90, :p1)
   seteprop!(g, :, 1 : 90, :p2)
   seteprop!(g, :, 1 : 90, :p3)
   seteprop!(g, :, 1 : 90, :p4)

   # Filter on vertex property
   g1 = @query(g |> filter(v.p1 <= 5))
   @test @query(g1 |> eachvertex(v.p1)) == collect(1:5)
   @test nv(g1) == 5
   @test ne(g1) == 20
   @test encode(g1) == [1,2,3,4,5]

   # Filter on edge property
   g2 = @query(g |> filter(e.p1 <= 45))
   @test @query(g2 |> eachedge(e.p1)) == collect(1:45)
   @test nv(g2) == 10
   @test ne(g2) == 45

   # Multiple vertex properties
   g3 = @query(g |> filter(v.p1 <=10, v.p2 >= 7))
   @test @query(g3 |> eachvertex(v.p1)) == collect(7:10)
   @test nv(g3) == 4
   @test ne(g3) == 12
   @test encode(g3) == [7,8,9,10]

   # Multiple edge properties
   g4 = @query(g |> filter(e.p1 <= 70, e.p2 >= 60))
   @test @query(g4 |> eachedge(e.p1)) == collect(60:70)
   @test nv(g4) == 10
   @test ne(g4) == 11

   # Mixed properties
   g5 = @query(g |> filter(v.p1 <= 5, e.p1 <= 20))
   @test nv(g5) == 5
   @test ne(g5) == 10

   # Source and target properites
   g6 = @query(g |> filter(s.p1 < t.p1))
   @test nv(g6) == 10
   @test ne(g6) == 45
end

@testset "Select" begin
   g = completegraph(10)
   # Set vertex properties
   setvprop!(g, :, 1 : 10, :p1)
   setvprop!(g, :, 1 : 10, :p2)
   setvprop!(g, :, 1 : 10, :p3)
   setvprop!(g, :, 1 : 10, :p4)

   # Set edge properties
   seteprop!(g, :, 1 : 90, :p1)
   seteprop!(g, :, 1 : 90, :p2)
   seteprop!(g, :, 1 : 90, :p3)
   seteprop!(g, :, 1 : 90, :p4)

   # Vertex props
   g1 = @query(g |> select(v.p1, v.p3))
   @test listvprops(g1) == [:p1, :p3]
   @test size(vdata(g1)) == (10,2)

   # Edge props
   g2 = @query(g |> select(e.p2, e.p4))
   @test listeprops(g2) == [:p2, :p4]
   @test size(edata(g2)) == (90,2)

   # Mixed props
   g3 = @query(g |> select(v.p1, v.p4, e.p2, e.p3))
   @test listvprops(g3) == [:p1, :p4]
   @test size(vdata(g3)) == (10,2)
   @test listeprops(g3) == [:p2, :p3]
   @test size(edata(g3)) == (90,2)
end


text:
"""
    Lx,Ly,Lz = makeangularmomentum(c,M,a,basis="Hermite")

Construct a compact representation of angular momentum operators in the given basis using step operators.

`basis`: name of the orthonomal mode space. Currently only `"Hermite"` is implemented.

`c`: cfield array on which `L_j` are to act.

`M`: vector of cutoff quantum numbers, i.e. `Mx,My = M` in 2D.

`a`: vector of oscillator lengths, i.e. `ax,ay = a` in 2D.

"""
function makeangularmomentum(c,M,a,basis="Hermite")
if basis=="Hermite"
  dim=length(size(c))
  if dim==1
    error("L not defined in 1D")
  elseif dim==2
    ax,ay = a
    Mx,My = M
    X,Px = ladderops(Mx,ax,basis)
    Y,Py = ladderops(My,ay,basis)
    Lz = -im*(X*c*Py' - Px*c*Y') #decide if return a function or matrices
    return Lz
  elseif dim==3
    ax,ay,az = a
    Mx,My,Mz = M
    X,Px = ladderops(Mx,ax,basis)
    Y,Py = ladderops(My,ay,basis)
    Z,Pz = ladderops(My,ay,basis)
    Lx = -im*(Y.*Pz - Z.*P̂y)
    Ly = -im*(Ẑ.*P̂x - X̂.*P̂z)
    Lz = -im*(X̂.*P̂y - Ŷ.*P̂x)
    return Lx,Ly,Lz
  end
else error("basis not implemented")
end
end


text:
struct test{R<:Real} # {} declares the condition on the type, it can be omitted, but if omitted, one can only constrain x and y to have the same types, not any other conditions on the type itself!
    x::R
    y::R
    test{R}(;x=4,y=7) where {R<:Real} = new(x,y)
end
#
obj = test{Int64}(x=3)
println("test ...")
println(obj.y)

test(::Type{R}; kwargs...) where {R} = test{R}(;kwargs...)
test(;kwargs...) = test(Int64;kwargs...)
test_(;kwargs...) = test{Int64}(;kwargs...)
# test_(;x::R , y::R) where{R<:Real} = test{R}(;x,y)
#
obj2 = test(Int64;x=2)
obj3 = test(x=3)
obj4 = test_(x=0)

function change_x!(y::Vector{Int64})
    y[1] = -1
end
x = [5,4]
change_x!(x)
println(x)

# function change_x!(x::Int64) #? how to make it mutable?
#     x = 1
# end
# x = 6
# change_x!(x)
# println(x)

# y = copy(x) #? is copy necessary?
# x[1] = 200
# println(y)
# println(x)

mutable struct x_m
    xx
end

function change_x!(x::x_m)
    x.xx = 1
end
x = x_m(10)
change_x!(x)
println(x.xx)
x.xx = 3
println(x.xx)


text:
struct CyclicContainer <: AbstractVector{String}
    c::Vector{String}
    n::Int
end
CyclicContainer(c) = CyclicContainer(c, 0)

Base.length(c::CyclicContainer) = typemax(Int)
Base.size(c::CyclicContainer) = size(c.c)
Base.getindex(c::CyclicContainer, i) = c.c[(i-1)%length(c.c) + 1]
function Base.getindex(c::CyclicContainer)
    c.n += 1
    c[c.n]
end
Base.iterate(c::CyclicContainer, i = 1) = iterate(c.c, i)

COLORSCHEME = [
    "#6D44D0",
    "#2CB3BF",
    "#1B1B1B",
    "#DA5210",
    "#03502A",
    "#866373",
]

COLORS = CyclicContainer(COLORSCHEME)
# LINESTYLES = CyclicContainer(["-", "--", ":", "-."])
LINESTYLES = CyclicContainer(["-", "--", "-."])
MARKERS = CyclicContainer(["o", "v", "s", "p", "X"])

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
using Test
using Dates
using ExtensibleScheduler

@testset "PeriodTrigger" begin
    @testset "TimeTrigger" begin
        td = Dates.Minute(5)
        trigger = Trigger(td)

        n_executed = 0

        function sample_action()
            n_executed += 1
        end
        trigger = Trigger(td)
        dt_start = DateTime(2010, 1, 1, 0, 0)
        dt = dt_start
        clock = SimClock(dt)
        sched = BlockingScheduler(clock=clock, delayfunc=NoSleep)
        action = Action(sample_action)

        add(sched, action, trigger)
        run_pending(sched)
        @test n_executed == 0

        for i in 1:11
            dt += Dates.Minute(1)
            #println("Clock set to $dt")
            set(clock, dt)
            run_pending(sched)
        end

        @test n_executed == 2

    end

    @testset "DateTrigger" begin
        td = Dates.Month(1)

        trigger = Trigger(td)

        n_executed = 0

        function sample_action()
            n_executed += 1
        end
        trigger = Trigger(td)
        dt_start = DateTime(2010, 1, 1, 0, 0)
        dt = dt_start
        clock = SimClock(dt)
        sched = BlockingScheduler(clock=clock, delayfunc=NoSleep)
        action = Action(sample_action)

        add(sched, action, trigger)
        run_pending(sched)
        @test n_executed == 0

        for i in 1:80
            dt += Dates.Day(1)
            #println("Clock set to $dt")
            set(clock, dt)
            run_pending(sched)
        end

        @test n_executed == 2
    end

end


text:
using BiobakeryUtils
using DataFrames
using Random
using Test
using Microbiome

@testset "Biobakery Utilities" begin
    abund = import_abundance_table("metaphlan_test.tsv")

    @test typeof(abund) <: DataFrame
    @test size(abund) == (42, 8)
    spec = taxfilter(abund, keepunidentified=false)
    @test size(spec) == (15, 8)
    phyl = taxfilter(abund, :phylum)
    @test size(phyl) == (2, 8)

    @test !any(occursin.("|", phyl[!, 1]))

    taxfilter!(abund, 2)
    @test abund == phyl

    d = rand(10, 10)
    dm = d + d'
    p = permanova(dm, repeat(["a", "b"], 5))
    @test typeof(p) == DataFrame
    @test size(p) == (3, 6)
end


text:
# This file is a part of SearchModels.jl

using Test

using SearchModels
import SearchModels: combine, mutate, config_type

struct PolyModelSpace <: AbstractSolutionSpace
    degree
end

const PolyModel = Vector{Float64}

config_type(c::PolyModel) = length(c)
Base.eltype(::PolyModelSpace) = PolyModel
Base.rand(space::PolyModelSpace) = randn(rand(space.degree) + 1)

function combine(a::PolyModel, b::PolyModel)
    [(a[i] + b[i]) / 2.0 for i in eachindex(a)]
end

function mutate(space::PolyModelSpace, c::PolyModel, iter)
    step = 1.0 + 1 / (1 + iter)
    [SearchModels.scale(c[i], s=step) for i in eachindex(c)]
end

function poly(coeff, x)::Float64
    s = coeff[1]
    @inbounds @simd for i in 2:length(coeff)
        m = coeff[i]
        s += coeff[i] * x^(i-1)
    end

    s
end

function create_model_1(n)
    X = collect(-10:0.1:10)
    coeff = rand(n + 1)
    X, [poly(coeff, x) for x in X], coeff
end

@testset "SearchModels.jl" begin
    X, y, coeff__ = create_model_1(4)

    space = PolyModelSpace(2:5)
    B = search_models(space, 300, SearchParams(
            bsize=64,
            mutbsize=32,
            crossbsize=32,
            tol=0.0,
            maxiters=100,
            verbose=true
        )
    ) do c
        s = zero(Float64)
        @inbounds @simd for i in eachindex(X)
            m = poly(c, X[i]) - y[i]
            s += m * m
        end

        s / length(X)
    end

    println("===== PolyModel -- poly: $(coeff__)")
    for (i, p) in enumerate(B[1:15])
        println(stderr, i => p)
    end
end


text:
using JD3
using Test

@testset "JD3.jl" begin
    # Write your tests here.
end


text:
include(joinpath(@__DIR__, "..", "src/dynamics", "flamingo", "visuals.jl"))
T = Float64
vis = Visualizer()
# render(vis)
open(vis)

s = get_simulation("flamingo", "flat_2D_lc", "flat")
model = s.model
env = s.env
const ContactImplicitMPC = Main
ref_traj = deepcopy(ContactImplicitMPC.get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/flamingo/gaits/gait_forward_36_4.jld2"),
    load_type = :split_traj_alt))

H = ref_traj.H
h = ref_traj.h
N_sample = 5
H_mpc = 15
h_sim = h / N_sample
H_sim = 2000#35000

# barrier parameter
κ_mpc = 1.0e-4

obj_mpc = quadratic_objective(model, H_mpc,
    q = [Diagonal(0.1 * ones(model.dim.q)) for t = 1:H_mpc+2],
    v = [Diagonal(0.001 * ones(model.dim.q)) for t = 1:H_mpc],
    u = [Diagonal(0.01 * ones(model.dim.u)) for t = 1:H_mpc-1])

p = linearized_mpc_policy(ref_traj, s, obj_mpc,
    H_mpc = H_mpc,
    N_sample = N_sample,
    κ_mpc = κ_mpc,
	ip_type = :interior_point,
	# mode = :configurationforce,
	mode = :configuration,
	newton_mode = :structure,
	# ip_type = :mehrotra,
    n_opts = NewtonOptions(
		# solver = :ldl_solver,
        r_tol = 3e-4,
		# verbose = true,
        max_iter = 5),
    mpc_opts = LinearizedMPCOptions(
        # live_plotting=true,
        # altitude_update = true,
        altitude_impact_threshold = 0.02,
        altitude_verbose = true,
        )
    )

q1_ref = copy(ref_traj.q[2])
q0_ref = copy(ref_traj.q[1])
q1_sim = SVector{model.dim.q}(q1_ref)
q0_sim = SVector{model.dim.q}(copy(q1_sim - (q1_ref - q0_ref) / N_sample))
@assert norm((q1_sim - q0_sim) / h_sim - (q1_ref - q0_ref) / h) < 1.0e-8

sim = simulator(s, q0_sim, q1_sim, h_sim, H_sim,
    p = p,
    ip_opts = InteriorPointOptions(
        r_tol = 1.0e-8,
        κ_init = 1.0e-8,
        κ_tol = 2.0e-8),
    sim_opts = SimulatorOptions(warmstart = true),
	ip_type = :interior_point,
    )

@time status = simulate!(sim)
l = 9
lu = 1
plt = plot(layout=(3,1), legend=false)
plot!(plt[1,1], hcat(Vector.(vcat([fill(ref_traj.q[i], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[1,1], hcat(Vector.([q[l:l] for q in sim.traj.q])...)', color=:blue, linewidth=1.0)
plot!(plt[2,1], hcat(Vector.(vcat([fill(ref_traj.u[i][lu:lu], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[3,1], hcat(Vector.(vcat([fill(ref_traj.γ[i][1:nc], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[2,1], hcat(Vector.([u[lu:lu] for u in sim.traj.u]*N_sample)...)', color=:blue, linewidth=1.0)
# plot!(plt[3,1], hcat(Vector.([γ[1:nc] for γ in sim.traj.γ]*N_sample)...)', color=:blue, linewidth=1.0)
# plot!(plt[3,1], hcat(Vector.([b[1:nb] for b in sim.traj.b]*N_sample)...)', color=:red, linewidth=1.0)

plot_surface!(vis, env, xlims=[-0.5, 1.5], ylims = [-0.5, 0.5])
# anim = visualize_robot!(vis, model, sim.traj, sample=10)
anim = visualize_meshrobot!(vis, model, sim.traj, sample=10)
anim = visualize_force!(vis, model, env, sim.traj, anim=anim, h=h_sim, sample=10)

convert_config(model_sim, sim.traj.q[1])

filename = "flamingo_100_steps"
MeshCat.convert_frames_to_video(
    "/home/simon/Downloads/$filename.tar",
    "/home/simon/Documents/$filename.mp4", overwrite=true)

convert_video_to_gif(
    "/home/simon/Documents/$filename.mp4",
    "/home/simon/Documents/$filename.gif", overwrite=true)

settransform!(vis["/Cameras/default"],
		compose(Translation(0.0, 0.5, -1.0),LinearMap(RotZ(-pi / 2.0))))

# Ghost
flamingo_ghost!(vis, sim, x -> 0.0)

# Animation
flamingo_animation!(vis, sim, x -> 0.0)


text:
struct RecursiveVec{T<:Union{Tuple,AbstractVector}}
    vecs::T

end
function RecursiveVec(arg1::AbstractVector{T}) where T
    if isbitstype(T)
        return RecursiveVec((arg1,))
    else
        return RecursiveVec{typeof(arg1)}(arg1)
    end
end
RecursiveVec(arg1, args...) = RecursiveVec((arg1, args...))

Base.getindex(v::RecursiveVec, i) = v.vecs[i]

Base.iterate(v::RecursiveVec, args...) = iterate(v.vecs, args...)

Base.IteratorEltype(::Type{RecursiveVec{T}}) where T = Base.IteratorEltype(T)
Base.IteratorSize(::Type{RecursiveVec{T}}) where T = Base.IteratorSize(T)

Base.eltype(v::RecursiveVec) = eltype(v.vecs)
Base.size(v::RecursiveVec) = size(v.vecs)
Base.length(v::RecursiveVec) = length(v.vecs)

Base.first(v::RecursiveVec) = first(v.vecs)
Base.last(v::RecursiveVec) = last(v.vecs)

Base.:-(v::RecursiveVec) = RecursiveVec(map(-,v.vecs))
Base.:+(v::RecursiveVec, w::RecursiveVec) = RecursiveVec(map(+,v.vecs, w.vecs))
Base.:-(v::RecursiveVec, w::RecursiveVec) = RecursiveVec(map(-,v.vecs, w.vecs))
Base.:*(v::RecursiveVec, a) = RecursiveVec(map(x->x*a, v.vecs))
Base.:*(a, v::RecursiveVec) = RecursiveVec(map(x->a*x, v.vecs))
Base.:/(v::RecursiveVec, a) = RecursiveVec(map(x->x/a, v.vecs))
Base.:\(a, v::RecursiveVec) = RecursiveVec(map(x->a\x, v.vecs))

function Base.similar(v::RecursiveVec)
    RecursiveVec(similar.(v.vecs))
end

function Base.copy!(w::RecursiveVec, v::RecursiveVec)
    @assert length(w.vecs) == length(v.vecs)
    @inbounds for i = 1:length(w.vecs)
        copyto!(w.vecs[i], v.vecs[i])
    end
    return w
end

function LinearAlgebra.mul!(w::RecursiveVec, a, v::RecursiveVec)
    @assert length(w.vecs) == length(v.vecs)
    @inbounds for i = 1:length(w.vecs)
        mul!(w.vecs[i], a, v.vecs[i])
    end
    return w
end

function LinearAlgebra.mul!(w::RecursiveVec, v::RecursiveVec, a)
    @assert length(w.vecs) == length(v.vecs)
    @inbounds for i = 1:length(w.vecs)
        mul!(w.vecs[i], v.vecs[i], a)
    end
    return w
end

function LinearAlgebra.rmul!(v::RecursiveVec, a)
    for x in v.vecs
        rmul!(x, a)
    end
    return v
end

function LinearAlgebra.axpy!(a, v::RecursiveVec, w::RecursiveVec)
    @assert length(w.vecs) == length(v.vecs)
    @inbounds for i = 1:length(w.vecs)
        axpy!(a, v.vecs[i], w.vecs[i])
    end
    return w
end
function LinearAlgebra.axpby!(a, v::RecursiveVec, b, w::RecursiveVec)
    @assert length(w.vecs) == length(v.vecs)
    @inbounds for i = 1:length(w.vecs)
        axpby!(a, v.vecs[i], b, w.vecs[i])
    end
    return w
end

LinearAlgebra.dot(v::RecursiveVec{T}, w::RecursiveVec{T}) where {T} = sum(dot.(v.vecs, w.vecs))
LinearAlgebra.norm(v::RecursiveVec) = norm(norm.(v.vecs))


text:
import CLIMAParameters
const CP = CLIMAParameters

Base.@kwdef struct ParameterBox{FT}
    molmass_dryair::FT
    gas_constant::FT
    new_parameter::FT
    # Derived parameters
    R_d::FT = gas_constant / molmass_dryair
end

function ParameterBox(toml_dict::CP.AbstractTOMLDict)

    aliases = ["molmass_dryair", "gas_constant", "new_parameter"]

    params = CP.get_parameter_values!(toml_dict, aliases, "ParameterBox")
    # Returns an array of `Pair`s for all given `aliases`

    FT = CP.float_type(toml_dict)
    return ParameterBox{FT}(; params...)
end

@testset "Example use case: parameter box" begin

    # [1.] read from file
    toml_file = joinpath(@__DIR__, "toml", "parambox.toml")
    toml_dict = CP.create_toml_dict(
        Float64;
        override_file = toml_file,
        dict_type = "alias",
    )

    # [2.] build
    param_set = ParameterBox(toml_dict)

    # [3.] log & checks(with warning)
    mktempdir(@__DIR__) do path
        logfilepath = joinpath(path, "logfilepath.toml")
        @test_logs (:warn,) CP.log_parameter_information(toml_dict, logfilepath)
    end

    # [4.] use
    # from default
    @test param_set.molmass_dryair ≈ 0.02897
    # overridden default
    @test param_set.gas_constant ≈ 4.0
    # derived in constructor
    @test param_set.R_d ≈ param_set.gas_constant / param_set.molmass_dryair
    # from toml
    @test param_set.new_parameter ≈ 19.99

end

max_stars_repo_name: climate-machine/CLIMAParameters

max_issues_repo_name: climate-machine/CLIMAParameters

max_forks_repo_name: climate-machine/CLIMAParameters


text:
struct TuplesIterator{G}
    length::Int
    generator::G
end
TuplesIterator(itr) = TuplesIterator(length(itr), itr)

# iteration interface
Base.iterate(t::TuplesIterator, args...) = iterate(t.generator, args...)
Base.length(t::TuplesIterator) = t.length
Base.eltype(::Type{<:TuplesIterator}) = Tuple

function Base.show(io::IO, tuples::TuplesIterator)
    Base.print(io, "$TuplesIterator(")
    Base.show(io, collect(Tuple, tuples.generator))
    Base.print(io, ")")
end

function Base.show(io::IO, ::MIME"text/plain", tuples::TuplesIterator)
    isempty(tuples) && return show(io, tuples)

    print(io, "$TuplesIterator with $(tuples.length) tuples:")
    # isempty(iter) && get(io, :compact, false) && return show(io, iter)
    # summary(io, iter)
    # isempty(iter) && return
    # print(io, ". ", isa(iter,KeySet) ? "Keys" : "Values", ":")
    _show_list(io, tuples)
end

function _show_list(io::IO, iter)
    limit = get(io, :limit, false)::Bool
    if limit
        sz = displaysize(io)
        rows, cols = sz[1] - 3, sz[2]
        rows < 2 && (print(io, " …"); return)
        cols < 4 && (cols = 4)
        cols -= 2 # For prefix "  "
        rows -= 1 # For summary
    else
        rows = cols = typemax(Int)
    end

    for (i, v) in enumerate(iter)
        print(io, "\n  ")
        i == rows < length(iter) && (print(io, "⋮"); break)

        if limit
            str = sprint(show, v, context=io, sizehint=0)
            str = Base._truncate_at_width_or_chars(str, cols, "\r\n")
            print(io, str)
        else
            show(io, v)
        end
    end
end


text:
############################################################################################
# Constants
"RNG for sampling based solutions"
const _RNG = Random.MersenneTwister(1)   # shorthand
Random.seed!(_RNG, 1)

"Tolerance for stochastic solutions"
const _TOL = 1.0e-6

"Number of samples"
N = 10^3

############################################################################################
# Initiate DefaultFlatten structs for strict and flexible (AD) type conversion:
outputtypes = [Float32, Float64]
flattentypes = [FlattenContinuous(), FlattenAll()]
unflattenmethods = [UnflattenStrict(), UnflattenFlexible()]

############################################################################################
#Probabilistic Parameters - Some selected distributions
struct ProbModel <: ModelName end
val_dist = (
    ## Normal distribution
    d_a1=Param(2.0, Distributions.Normal()),
    d_a2=Param(Float32(2.0), Distributions.Normal{Float32}(Float32(1.0), Float32(2.0))),
    d_a3=Param([2.0, 2.0], [Distributions.Normal(), Distributions.Normal(2.0, 3.0)]),
    d_a5=Param(
        Float16.([52.0, 53.0]),
        [
            Distributions.Normal{Float16}(Float16(1.0), Float16(2.0)),
            Distributions.Normal{Float16}(Float16(1.0), Float16(2.0)),
        ],
    ),
    d_a6=Param(
        [[2.0, 2.0], [2.0, 2.0]],
        [
            [Distributions.Normal(), Distributions.Normal(2.0, 3.0)],
            [Distributions.Normal(), Distributions.Normal(2.0, 3.0)],
        ],
    ),
    d_a7=Param(5.0, Distributions.truncated(Distributions.Normal(2.0, 3.0), 0.0, 10.0)),
    d_a8=Param(
        [2.0, 2.0],
        [
            Distributions.truncated(Distributions.Normal(2.0, 3.0), 0.0, 10.0),
            Distributions.truncated(Distributions.Normal(2.0, 3.0), 0.0, 10.0),
        ],
    ),
    ## Multivariate Normal
    d_b1=Param(
        [2.0, 2.0], Distributions.MvNormal([2.0, 3.0], Diagonal(map(abs2, [1.0, 1.0])))
    ),
    d_b2=Param(
        [[2.0, 2.0], [2.0, 2.0]],
        [
            Distributions.MvNormal([2.0, 3.0], Diagonal(map(abs2, [1.0, 1.0]))),
            Distributions.MvNormal(Diagonal(map(abs2, [1.0, 1.0]))),
        ],
    )
    ,
    d_b3=Param(
        [Float32.([54.0, 55.0]), Float32.([56.0, 57.0])],
        [
            Distributions.MvNormal(Float32.([1, 1.0]), Float32.([5.0 0.3; 0.3 5.0])),
            Distributions.MvNormal(Float32.([1, 1.0]), Float32.([5.0 0.3; 0.3 5.0])),
        ],
    ),
    ## Gamma
    d_c1=Param(2.0, Distributions.Gamma(2, 3)),
    d_c2=Param([2.0, 2.0], [Distributions.Gamma(2, 3), Distributions.Gamma(2, 3)]),
    ## Exponential
    d_d1=Param(2.0, Distributions.Exponential(2)),
    d_d2=Param([2.0, 2.0], [Distributions.Exponential(2), Distributions.Exponential(2)]),
    ## Dirichlet Distribution
    d_e1=Param([0.1, 0.9], Distributions.Dirichlet(2, 2.0)),
    d_e2=Param(
        [[0.1, 0.9], [0.8, 0.2]],
        [Distributions.Dirichlet(2, 2.0), Distributions.Dirichlet(2, 2.0)],
    ),
    d_e3=Param([0.2, 0.3, 0.5], Distributions.Dirichlet(3, 3.0)),
    d_e4=Param(
        [[0.2, 0.3, 0.5], [0.2, 0.3, 0.5]],
        [Distributions.Dirichlet(3, 3.0), Distributions.Dirichlet(3, 3.0)],
    ),
    ## Categorical
    d_f1=Param(5, Distributions.Categorical(100)),
    #f2 = Param(5, Distributions.Categorical( Float32.( repeat([1/100], 100) ) ),
    d_f3=Param([10, 10], [Distributions.Categorical(100), Distributions.Categorical(100)]),
    d_f4=Param(
        [[34, 34], [34, 34]],
        [
            [Distributions.Categorical(40), Distributions.Categorical(40)],
            [Distributions.Categorical(40), Distributions.Categorical(40)],
        ],
    ),
    ## Poisson
    d_g1=Param(5, Distributions.Poisson(100)),
    d_g2=Param([10, 10], [Distributions.Poisson(100), Distributions.Poisson(100)]),
    ## Negative Binomial ~ Does not work for Zygote (foreign function call)
    #        d_h1 = Param(5, Distributions.NegativeBinomial(.5, .5)),
    #        d_h2 = Param([10, 10], [Distributions.NegativeBinomial(.5, .5), Distributions.NegativeBinomial(.5, .5)]),
    ## Uniform
    d_i1=Param(1.0, Distributions.Uniform(-5.0, 5.0)),
    #~ Not correctly sampled as Float16 from Distributions package.
    #d_i2 = Param(Float16(5.), Distributions.Uniform(Float16(-10), Float16(10))),
    d_i3=Param(
        [3.0, 3.0], [Distributions.Uniform(-5.0, 5.0), Distributions.Uniform(-5.0, 5.0)]
    ),
    ## InverseWishart ~ Does not work with ReverseDiff due to Diag{} type conversion from Bijector
    #        d_j1 = Param([15. .16 ; .16 18.], Distributions.InverseWishart(10., [1. 0. ; 0. 1.])),
    #        d_j2 = Param([[19. .20 ; .20 22.], [23. .24 ; .24 26.]],[Distributions.InverseWishart(10., [1. 0. ; 0. 1.]), Distributions.InverseWishart(10., [1. 0. ; 0. 1.])]),
    #        d_j3 = Param([15. .16 .16 ; .16 18. .16 ; .16 .16 20.], Distributions.InverseWishart(10., [1. 0. .0 ; 0 1 0 ; 0 0. 1.])),
    #        d_j4 = Param([[15. .16 .16 ; .16 18. .16 ; .16 .16 20.], [15. .16 .16 ; .16 18. .16 ; .16 .16 20.]],[Distributions.InverseWishart(10., [1. 0. .0 ; 0 1 0 ; 0 0. 1.]), Distributions.InverseWishart(10., [1. 0. .0 ; 0 1 0 ; 0 0. 1.])]),
    ## Cauchy
    d_k1=Param(5.0, Distributions.Cauchy(5.0, 10.0)),
    d_k2=Param(
        [2.0, 2.0], [Distributions.Cauchy(5.0, 10.0), Distributions.Cauchy(5.0, 10.0)]
    ),
    ## LKJ
    d_l1=Param([1.0 0.16; 0.16 1.0], Distributions.LKJ(2, 1.0)),
    d_l2=Param(
        [[1.0 0.20; 0.20 1.0], [1.0 0.24; 0.24 1.0]],
        [Distributions.LKJ(2, 1.0), Distributions.LKJ(2, 1.0)],
    ),
    d_l3=Param([1.0 0.16 0.16; 0.16 1.0 0.16; 0.16 0.16 1.0], Distributions.LKJ(3, 1.0)),
    d_l4=Param(
        [
            [1.0 0.16 0.16; 0.16 1.0 0.16; 0.16 0.16 1.0],
            [1.0 0.16 0.16; 0.16 1.0 0.16; 0.16 0.16 1.0],
        ],
        [Distributions.LKJ(3, 1.0), Distributions.LKJ(3, 1.0)],
    ),
    ## Truncated
    d_m1=Param(10.0, Distributions.truncated(Distributions.Normal(5.0, 10.0), 0.0, 20.0)),
    d_m2=Param(10.0, Distributions.truncated(Distributions.Cauchy(5.0, 10.0), 0.0, 20.0)),
    ## Product Distributions -> Check if ReverseDiff works with this structure
    d_n1=Param(
        [[i / 10 for i in Base.OneTo(3)] for _ in Base.OneTo(4)],
        [
            Distributions.product_distribution([Uniform(-1, 1) for _ in Base.OneTo(3)]) for
            _ in Base.OneTo(4)
        ],
    ),
)
#val_dist_length = 58

############################################################################################
# Nested Probabilistic Parameter/Models - Some selected distributions
val_dist_nested = (;
    nd_a1=(a=val_dist.d_a1, b=val_dist.d_f1, c=val_dist.d_i3),
    nd_a2=(
        a=(b=val_dist.d_a1, c=val_dist.d_f1, d=val_dist.d_i3),
        e=val_dist.d_a2,
        f=val_dist.d_e2,
    ),
    nd_a3=(
        a=(b=val_dist.d_a1, c=(g=val_dist.d_f1,), d=val_dist.d_i3),
        e=val_dist.d_a2,
        f=val_dist.d_e2,
    ),
)

############################################################################################
# Non-Probabilistic Parameter (Experimental!):
struct ConstrainedModel <: ModelName end
val_constrained = (
    ## Fixed
    # Floats
    c_a1=Param(1.0, Fixed()),
    c_a2=Param([2.0, 2.0], Fixed()),
    c_a3=Param([14.0 0.14; 0.14 14.0], Fixed()),
    c_b1=Param(Float32(26.0), Fixed()),
    c_b2=Param(Float16.([28.0, 28.0]), Fixed()),
    # Integer - Should not be flattened
    c_c1=Param(34, Fixed()),
    c_c2=Param([34, 34], Fixed()),
    # Fixed buffer parameter
    c_d1=Param(rand(1:100, 100), Fixed()),
    c_d2=Param(randn(100), Fixed()),
    c_d3=Param(rand(1:100, 10, 10), Fixed()),
    c_d4=Param(randn(10, 10), Fixed()),
    # Arbitrary AbstractArrays
    c_e1=Param([[1.0, 2.0], zeros(3, 2, 1)], Fixed()),
    c_e2=Param([[1, 2], [1, 2]], Fixed()),
    c_e3=Param([[1.0, 2.0], [1.0, 2.0]], Fixed()),
    c_e4=Param([(rand(), rand()) for _ in 1:10], Fixed()),
    c_e5=Param([(rand(1:10), rand(1:10)) for _ in 1:10], Fixed()),
    ## Unconstrained Parameter
    # Floats
    u_a1=Param(1.0, Unconstrained()),
    u_a2=Param([2.0, 2.0], Unconstrained()),
    u_a3=Param([14.0 0.14; 0.14 14.0], Unconstrained()),
    u_b1=Param(Float32(26.0), Unconstrained()),
    u_b2=Param(Float16.([28.0, 28.0]), Unconstrained()),
    # Fixed buffer parameter
    u_d1=Param(rand(1:100, 100), Unconstrained()),
    u_d2=Param(randn(100), Unconstrained()),
    u_d3=Param(rand(1:100, 10, 10), Unconstrained()),
    u_d4=Param(randn(10, 10), Unconstrained()),
    # Arbitrary AbstractArrays
    u_e1=Param([[1.0, 2.0], zeros(3, 2, 1)], Unconstrained()),
    u_e2=Param([[1, 2], [1, 2]], Unconstrained()),
    u_e3=Param([[1.0, 2.0], [1.0, 2.0]], Unconstrained()),
    u_e4=Param([(rand(), rand()) for _ in 1:10], Unconstrained()),
    u_e5=Param([(rand(1:10), rand(1:10)) for _ in 1:10], Unconstrained()),
    ## Constrained Parameter -> Scalar only
    con_a1=Param(1.0, Constrained(0.0, 2.0)),
    con_b1=Param(Float32(26.0), Constrained(Float32(20.0), Float32(30.0))),
)

################################################################################
# Parameter for an example Model
_σ = Distributions.truncated(Distributions.Cauchy(10.0, 10), 0.0, 30.0)
_iwish2 = Distributions.InverseWishart(10.0, [1.0 0.0; 0.0 1.0])
_iwish3 = Distributions.InverseWishart(10.0, [1.0 0.0 0.0; 0 1 0; 0 0.0 1.0])
_ρ = [1.0 0.25; 0.25 1.0]
struct ExampleModel <: ModelName end
_val_examplemodel = (
    μ1=Param(1.0, Distributions.Normal()),
    μ2=Param([2.0, 3.0], Distributions.MvNormal(Diagonal(map(abs2, [1.0, 1.0])))),
    μ3=Param([4.0, 5.0], [Distributions.Normal(), Distributions.Normal()]),
    μ4=Param(
        [[6.0, 7.0], [8.0, 9.0]],
        [
            Distributions.MvNormal([2.0, 3.0], Diagonal(map(abs2, [1.0, 1.0]))),
            Distributions.MvNormal([1.0, 2.0], [3.0 0.4; 0.4 5.0]),
        ],
    ),
    σ1=Param(10.0, Distributions.Gamma()),
    σ3=Param([5.0, 5.0], [Distributions.Gamma(), Distributions.Gamma()]),
    σ2=Param([5.0, 5.0], [_σ, _σ]),
    ρ2=Param(copy(_ρ), Distributions.LKJ(2, 1.0)),
    σ4=Param([[5.0, 5.0], [5.0, 5.0]], [[_σ, _σ], [_σ, _σ]]),
    ρ4=Param([copy(_ρ), copy(_ρ)], [Distributions.LKJ(2, 1.0), Distributions.LKJ(2, 1.0)]),
    p=Param([0.2, 0.3, 0.5], Distributions.Dirichlet(3, 3.0)),
)
#_val_examplemodel_length = 23

################################################################################
# Parameter that may work in lower dimensions
struct LowerDims <: ModelName end
_val_lowerdims = (;
    ## Dirichlet
    ldim_e1=Param([0.1, 0.9], Distributions.Dirichlet(2, 2.0)),
    ldim_e2=Param(
        [[0.1, 0.9], [0.8, 0.2]],
        [Distributions.Dirichlet(2, 2.0), Distributions.Dirichlet(2, 2.0)],
    ),
    ldim_e3=Param([0.2, 0.3, 0.5], Distributions.Dirichlet(3, 3.0)),
    ldim_e4=Param(
        [[0.2, 0.3, 0.5], [0.2, 0.3, 0.5]],
        [Distributions.Dirichlet(3, 3.0), Distributions.Dirichlet(3, 3.0)],
    ),
    ## LKJ
    ldim_l1=Param([1.0 0.16; 0.16 1.0], Distributions.LKJ(2, 1.0)),
    ldim_l2=Param(
        [[1.0 0.20; 0.20 1.0], [1.0 0.24; 0.24 1.0]],
        [Distributions.LKJ(2, 1.0), Distributions.LKJ(2, 1.0)],
    ),
    ldim_l3=Param([1.0 0.16 0.16; 0.16 1.0 0.16; 0.16 0.16 1.0], Distributions.LKJ(3, 1.0)),
    ldim_l4=Param(
        [
            [1.0 0.16 0.16; 0.16 1.0 0.16; 0.16 0.16 1.0],
            [1.0 0.16 0.16; 0.16 1.0 0.16; 0.16 0.16 1.0],
        ],
        [Distributions.LKJ(3, 1.0), Distributions.LKJ(3, 1.0)],
    ),
    ## InverseWishart
    #        ldim_j1 = Param([15. .16 ; .16 18.], _iwish2),
    #        ldim_j2 = Param([[19. .20 ; .20 22.], [23. .24 ; .24 26.]],[_iwish2, _iwish2]),
    #        ldim_j3 = Param([15. .16 .16 ; .16 18. .16 ; .16 .16 20.], _iwish3),
    #        ldim_j4 = Param([[15. .16 .16 ; .16 18. .16 ; .16 .16 20.], [15. .16 .16 ; .16 18. .16 ; .16 .16 20.]],[_iwish3, _iwish3]),
)


text:
using LoadflowPkgDemo
using Test

@testset "LoadflowPkgDemo.jl" begin
    # Write your tests here.
end


text:
######
# This file is part of the BivariateCopulas.jl package.
#
#   Definition of bivariate copula type, Joint distribution type and methods
#
#           University of Liverpool, Institute for Risk and Uncertainty
#
#                                           Author: Ander Gray
#                                           Email:  ander.gray@liverpool.ac.uk
######

###
#   Notes: This is a module for performing tests with copulas for plotting, sampling and
#   making bivariate joint distributions
#
#   Would it be possible to represent a copula as an expansion of orthogonal
#   polynomials? Could be used as a sentivitivty analysis for dependancy. Something equivilant
#   to disperive Monte Carlo, where not only the correlatin coefficient is changed but the entire
#   dependancy structure. Could this expansion be used in the arithmetic?
#
#   Or can we do it as sums of the Frechet bounds? (Shuffles of M)
#   See: Nelson page 68.
#
#
###

###
#
#   Known bugs:
#                   ->  Differentiating with the interpolators is really ugly
#
###

#import PyPlot: PyPlot.plot

mutable struct copula <: AbstractCopula

    cdf ::  Array{Float64,2}
    density :: Union{Array{Float64,2},Missing}
    func :: Union{Function,Missing}
    param :: Union{Float64,Missing}

    function copula(cdf = missing; density = missing, func = missing, param = missing)

        if (ismissing(cdf) && ismissing(density)) throw(ArgumentError("Cdf and/or density must be provided")) end
        if (ismissing(cdf)); cdf = calcCdf(density); end
        C = new(cdf,density,func,param)
        if (ismissing(density)); density = calcDensity(C); end

        return new(cdf,density, func, param)

    end
end

function (obj::copula)(X,Y, useInterp = false)             # Evaluating the copula

    #if (!issorted(X) || !issorted(Y)) throw(ArgumentError("cdf evaluation request must be given in assending order"));end

    if !useInterp
        if !ismissing(obj.func)
            if !ismissing(obj.param)
                return obj.func(X,Y,obj.param)
            end
            return obj.func(X,Y)
        end
    end

    cdf = interpolate(obj.cdf, BSpline(Linear()))

    uu = range(0,1,length = n)
    sCdf = Interpolations.scale(cdf, uu, uu);
    return sCdf(X,Y);
end

cdf(C:: copula, X, Y) = C(X,Y)

function density(C :: copula, X, Y)

    den = interpolate(C.density, BSpline(Quadratic(Line(OnCell()))))
    uu = range(0,1,length = n)
    sDen = Interpolations.scale(den, uu, uu);
    return sDen(X, Y);
end

function (obj::copula)(X :: ContinuousUnivariateDistribution,Y :: ContinuousUnivariateDistribution)             #Constructing a joint distribution using the copula
    return Joint(X, Y,obj)
end

function sample(C :: AbstractCopula, N ::Int64 = 1; plot = false, useInterp = false)

    # Set true to use interpolator
    useInterp = useInterp
    if ismissing(C.func) func = 0; useInterp = true else func = C.func end

    if plot return samplePlot(C,N);end
    func == Gau && return CholeskyGaussian(N, C.param)   # Use Cholesky decompostition of the Cov matrix for gaussian copula sampling
    func == Cla && return _claytonsample(N, C.param)

    x = rand(N);    y = rand(N);
    ux = x;         uy = zeros(N);
    m = n;

    if func   == perf return hcat(ux,ux); end
    if func   == opp return hcat(ux,1 .- ux); end
    if func   == indep return hcat(ux,y); end

    # if (!ismissing(C.func) && n < 10^5) m = 10^5; end  # If function is provided more acurate sampling, otherwise must use interpolator
    e = sqrt(eps());      ygrid = range(0,stop = 1,length = m);

    #if (C.func == Gau) useInterp = true ;end           # Should you not want to use Cholesky and interpolator has better performance for gaussian

    if func !== 0; conditional = (C(x .+ e/2,ygrid, useInterp) - C(x .- e/2,ygrid, useInterp))./e
    else

    end
    conditional = (C(x .+ e/2,ygrid, useInterp) - C(x .- e/2,ygrid, useInterp))./e
    for i=1:N
        #conditional = (C(x[i] + e/2,ygrid, true) - C(x[i] - e/2,ygrid, true))./e
        #conditional[1] = 0; conditional[end]=1;       # Small errors sometimes happen
        #if (conditional[end] < y[i]) uy[i] = 1;      # Should linearly interpolate between cond[end] and 1
        #elseif (conditional[1] > y[i]) uy[i] = 0;      # Should linearly interpolate between cond[end] and 1
        #else uy[i] = findlast(conditional .<= y[i]) / m; end        # Finds Nothing??
        #yitp_linear = interpolate(conditional, BSpline(Linear()))
        #uy[i] = yitp_linear(y[i] * m);

        #uy[i] = findlast(conditional .<= y[i])/m; end
        thisCond = conditional[i,:]
        thisGrid = ygrid

        if !issorted(thisCond) thisCond = sort(thisCond) end

        idx = unique(z ->thisCond[z], 1:n)    # Some magic: https://discourse.julialang.org/t/unique-indices-method-similar-to-matlab/34446
                                                    # Gets unique indicies from the vector

        thisCond = thisCond[idx]
        thisGrid = ygrid[idx]

        inpl = LinearInterpolation(thisCond, thisGrid, extrapolation_bc = Flat())
        uy[i] = inpl(y[i])
    end
    return hcat(ux,uy);
end

rand(C :: AbstractCopula, N :: Int64 = 1; plot= false) = sample(C, N, plot=plot)

function CholeskyGaussian(N = 1, correlation = 0)

    Cov = ones(2,2); Cov[1,2] = Cov[2,1] = correlation;
    a = cholesky(Cov).L;
    z = transpose(rand(Normal(),N,2))
    x = a * z;
    u = transpose(cdf.(Normal(),x))

    return hcat(u[:,1],u[:,2])
end

function _claytonsample(n, t)
    r = rand(n, 2)
    t == 0 && return r
    for i in 1:n
        v = r[i, 2]
        u = r[i, 1]
        r[i, 2] = (1 - u^(-t) + (u^(1 + t)*v)^(-(t/(1 + t))))^(-1/t)
    end
    return r
end

function conditional(C :: AbstractCopula, xVal :: Real; plot = false)   # May also work for joint? xVal will be take from invCdf(M1, u1)

    e = 1/n;        ygrid = range(0,stop=1,length = n);
    if !ismissing(C.func) e = (sqrt(eps())); end

    conditional = (C(xVal + e/2,ygrid) - C(xVal - e/2,ygrid))./e;
    if (plot) plot2(C.cdf, conditional, xVal) end
    return conditional

end

function calcCdf(density :: Array{<:Real,2})
    cdf = zeros(size(density))
    for i = 2:n
        for j = 2:n
            cdf[i,j] =  (density[i,j]/(n^2) + cdf[i, j-1] + cdf[i-1, j] - cdf[i-1, j-1]);
        end
    end
    return cdf
end

function calcDensity(C :: AbstractCopula)

    if ismissing(C.func) return calcDensity(C.cdf); end

    e = sqrt(sqrt(eps()));

    x = range(0 , stop = 1 - e    ,length = n);
    y = range(0 , stop = 1 - e    ,length = n);

    der1 = ( C(x .+ e, y)       - C(x, y)     )/e      # df/dx
    der2 = ( C(x .+ e, y .+ e ) - C(x, y .+ e) )/e     # d(df/dx)/dy
    density = (der2-der1)/e

    return density
end

function calcDensity(cdf :: Array{<:Real,2}, xRange = range(0,1,length = n), yRange = range(0,1,length = n))

    xStep = Float64(xRange.step); yStep = Float64(yRange.step);

    density = zeros(size(cdf));
    for i = 1:n-1
        for j = 1:n-1
            density[i,j] = (cdf[i,j] + cdf[i+1,j+1] - cdf[i+1,j] - cdf[i,j+1])/(xStep*yStep);

            #   -> Possible alternative:
            #der1 = ( cdf[i + 1, j]      - cdf[i, j]     )/(xStep^2)             # df/dx
            #der2 = ( cdf[i + 1, j + 1 ] - cdf[i, j + 1] )/(yStep^2)             # d(df/dx)/dy
            #density[i,j] = (der2 - der1)
        end
    end

    # Problem: How do you get the end values?
    #   Two solutions:

    #   -> If copula, it's symetric (exact answer exept at two endpoints)
    #density[end,:] = reverse(density[:,1]);
    #density[:,end] = reverse(density[1,:]);

    #   -> If joint, it may not be symetric
    density[end,:] = density[end-1,:];
    density[:,end] = density[:,end-1];
    return density
end

function M()
    x = y = range(0,stop=1,length = n);
    return copula(perf(x,y), density = Matrix{Float64}(I, n, n), func = perf);
end

function W()
    x = y = range(0,stop=1,length = n);

    den = zeros(n,n)
    for i = 1:n
        for j = 1:n
            if (i-1) == n-(j-1); den[i,j] = 1; end
        end
    end
    return copula(opp(x,y), density = den , func = opp);
end

function Pi()
    x = y = range(0,stop=1,length = n);
    return copula(indep(x,y),density= ones(n,n), func = indep);
end

function Frank(s = 1)                                     #   Frank copula
    x = y = range(0,stop = 1,length = n);                          #   s>0; 0 for perfect, 1 for indep, inf for oposite
    return copula(F(x,y,s), func = F, param = s);
end

function Clayton(t = 0)                                  #   Clayton copula
    x = y = range(0, stop=1, length = n);                          #   t>-1; -1 for opposite, 0 for indep and inf for perfect
    return copula(Cla(x,y,t), func = Cla, param = t);
end

function Gaussian(corr = 0)
    x = y = range(0,stop=1,length = n);
    cdf = Gau(x,y,corr);
    #cdf[:,end] = cdf[end,:] = x;    # This removes the NaNs from Liams mvNormCdf. We know the marginals must be uniform
    return copula(cdf, func = Gau, param = corr);
end

function τCopula( τ = 0 )   # Imprecise copula from kendal tau

    x = y = range(0,stop = 1,length = n);

    cdf = kenLB(x,y,τ)

    return copula(cdf, func = kenLB, param = τ)

end

KendalCopula(τ = 0) = τCopula(τ)

function ρCopula( ρ = 0 ) # Imprecise copula from Spearman rho

    x = y = range(0,stop = 1,length = n);

    cdf = spearLB(x,y,ρ)

    return copula(cdf, func = spearLB, param = ρ)
end

SpearmanCopula(ρ = 0 ) = ρCopula( ρ )

indep(X, Y) = [x*y for x in X, y in Y];
perf(X, Y)  = [min(x,y) for x in X, y in Y];
opp(X, Y)   = [max(x+y-1,0) for x in X, y in Y];
F(X,Y,s = 1)      = [log(1+(s^x-1)*(s^y-1)/(s-1))/log(s) for x in X, y in Y]
Cla(X,Y, t = 0) = [max((x^(-t)+y^(-t)-1)^(-1/t),0) for x in X, y in Y]
Gau(X,Y, corr = 0)  = [bivariate_cdf(quantile.(Normal(),x),quantile.(Normal(),y), corr) for x in X, y in Y];

#kenLB(x, y, τ) = max(0, x+y-1, 0.5*( (x+y)-sqrt( (x-y)^2)+1-τ ) )
kenUB(x, y, τ) = min(x, y, 0.5*( (x+y-1) + sqrt( (x+y-1)^2 +1+τ )) )
kenLB(x, y, τ) = [x - kenUB(x, 1 - y, - τ) for x in x, y in y]

#spearϕ(a, b) = 1/6 * ( (max(0, 9*b + 3*sqrt( 9*b^2 - 3*a^6)))^(1/3) + ( max(0,9*b - 3*sqrt(9*b^2 - 3*a^6)))^(1/3) )

function spearϕ(a, b)
    A = 9*b
    B = max(9*b^2 - 3*a^6, 0)
    C = (max(0, A + 3*sqrt(B)))^(1/3)
    D = (max(0, A - 3*sqrt(B)))^(1/3)
    return 1/6 * (C + D)
end

#spearLB(x, y, ρ) = max(0, x + y - 1, (x + y)/2 - spearϕ(x - y, 1 - ρ))
spearUB(x, y, ρ) = max( opp(x,y)[1], min(x, y, (x + y -1)/2 + spearϕ(x + y - 1, 1 + ρ)))
spearLB(x, y, τ) = [x - spearUB(x, 1 - y, - τ) for x in x, y in y]

function GauCopulaDen(corr = 0, X = range(bOt, tOp, length=n), Y = range(bOt, tOp, length=n))

    COR = [1. corr; corr 1.];
    I = [1. 0; 0 1.]
    den = zeros(length(X),length(Y))
    for i = 1:length(X)
        for j = 1:length(Y)
            Vec = [invCdf(Normal(),X[i]); invCdf(Normal(),Y[j])]
            den[i,j] = 1/sqrt(det(COR)) * exp(-1/2 * Vec' * (inv(COR) - I) *Vec)
        end
    end
    return den
end

Arch1(X,Y,theta = 0) = [1 - ((1-x)^theta + (1-y)^theta - (1-x)^theta * (1-y)^theta )^theta for x in X, y in Y];

mvNormCdf(X,Y,cor) = [bivariate_cdf(x,y,cor) for x in X, y in Y];

###
#   Copula rotations
###
function rotate90(C :: copula)

    C_func = C.func;
    C_param = C.param;
    if !ismissing(C_func)
        if !ismissing(C_param)
            Rotate90(xs, ys, C_param) = [y - C_func(1 - x, y, C_param)[1] for x in xs, y in ys]
            u = v = range(0, stop=1,length = n);
            cdf = Rotate90(u, v, C_param)
            cdf[1,:] .= 0
            cdf[:,1] .= 0
            return copula(cdf, func = Rotate90, param = C_param)
        else
            Rotate90(xs, ys) = [y - C_func(1 - x, y)[1] for x in xs, y in ys]
            u = v = range(0, stop=1, length = n);
            cdf = Rotate90(u, v)
            return copula(cdf, func = Rotate90)
        end
    end

    ###
    #   Here we need to figure out how to rotate the cdf matrix
    ###

    throw(ArgumentError("Currently can only rotate copulas with functions defined"))

end

function rotate180(C :: copula)

    C_func = C.func;
    C_param = C.param;
    if !ismissing(C_func)
        if !ismissing(C_param)
            Rotate180(xs, ys, C_param) = [x + y - 1 + C_func(1 - x, 1 - y, C_param)[1] for x in xs, y in ys]
            u = v = range(0, stop=1,length = n);
            cdf = Rotate180(u, v, C_param)
            cdf[1,:] .= 0
            cdf[:,1] .= 0
            return copula(cdf, func = Rotate180, param = C_param)
        else
            Rotate180(xs, ys) = [x + y - 1 + C_func(1 - x, 1 - y)[1] for x in xs, y in ys]
            u = v = range(0, stop=1, length = n);
            cdf = Rotate180(u, v)
            return copula(cdf, func = Rotate180)
        end
    end

    ###
    #   Here we need to figure out how to rotate the cdf matrix
    ###

    throw(ArgumentError("Currently can only rotate copulas with functions defined"))

end

function rotate270(C :: copula)

    C_func = C.func;
    C_param = C.param;
    if !ismissing(C_func)
        if !ismissing(C_param)
            Rotate270(xs, ys, C_param) = [x - C_func(x, 1 - y, C_param)[1] for x in xs, y in ys]
            u = v = range(0, stop=1,length = n);
            cdf = Rotate270(u, v, C_param)
            cdf[1,:] .= 0
            cdf[:,1] .= 0
            return copula(cdf, func = Rotate270, param = C_param)
        else
            Rotate270(xs, ys) = [x - C_func(x, 1 - y)[1] for x in xs, y in ys]
            u = v = range(0, stop=1, length = n);
            cdf = Rotate270(u, v)
            return copula(cdf, func = Rotate270)
        end
    end

    ###
    #   Here we need to figure out how to rotate the cdf matrix
    ###

    throw(ArgumentError("Currently can only rotate copulas with functions defined"))

end

mutable struct Joint <: AbstractJoint

    marginal1 :: ContinuousUnivariateDistribution
    marginal2 :: ContinuousUnivariateDistribution
    copula  :: AbstractCopula

    xRange
    yRange

    cdf ::  Union{Array{Float64,2}, Missing}     # Or should the density/ cdf be calculated on the fly from the copula?
    density :: Union{Array{Float64,2}, Missing}

    function Joint(marginal1, marginal2, copula)

        xRange = range(invCdf(marginal1,bOt), stop=invCdf(marginal1,tOp), length = n)
        yRange = range(invCdf(marginal2,bOt), stop=invCdf(marginal2,tOp), length = n)

        j = new(marginal1, marginal2, copula, xRange, yRange, missing)
        cdf = j(xRange,yRange)
        j = new(marginal1, marginal2, copula, xRange, yRange, cdf)
        density = calcDensity(j)                                    # Looks really ugly when using the interpolator for a marginal

        return new(marginal1, marginal2, copula, xRange, yRange, cdf, density)

    end
end

function (obj :: Joint)(X,Y)    # Evaluating the joint cdf
                                # BiLinear interpolator if cdf is provided

    if !ismissing(obj.copula.func)
        return obj.copula(cdf.(obj.marginal1,X), cdf.(obj.marginal2,Y))  # This is the case where we have the functions
    end
    if ismissing(obj.cdf)
        return obj.copula(cdf.(obj.marginal1,X), cdf.(obj.marginal2,Y))
    end

    cdf1 = interpolate(obj.cdf, BSpline(Quadratic(Line(OnCell()))))
    sCdf1 = Interpolations.scale(cdf1, obj.xRange, obj.yRange);

    return sCdf1(X, Y);

    # Should use new scaling from interpolatorrs
end

cdf(J::Joint, X, Y) = J(X, Y)

function conditional(J :: Joint, xVal :: Real; plot = false)   # May also work for joint? xVal will be take from invCdf(M1, u1)

    e = 1/n;        ygrid = J.yRange;
    if !ismissing(J.copula.func) e = (sqrt(eps())); end

    conditional = (J(xVal + e/2,ygrid) - J(xVal - e/2,ygrid))./e;

    conditional = conditional/pdf(J.marginal1,xVal)

    if (plot) plot2(J.cdf, conditional, xVal) end
    return conditional

end

function density2(J :: Joint, X = J.xRange, Y = J.yRange)

    den = interpolate(J.density, BSpline(Quadratic(Line(OnCell()))))
    sDen = Interpolations.scale(den, X, Y);
    return sDen(X,Y);

end

function density(J :: Joint,X = J.xRange, Y = J.yRange)        # An alternative way to construct the Joint. fxy / fx*fy = c[Fx,Fy]. Where c is the copula Density
                                                                # More error found against Multivariate gaussian pdf. Due to interpolation error in the copula density.
    denCop = density(J.copula, cdf.(J.marginal1,X), cdf.(J.marginal2,Y))
    return denCop .* [pdf(J.marginal1, x) * pdf(J.marginal2, y) for x in X, y in Y]

end

function sample(J :: AbstractJoint, N :: Int64 =1)

    copulaSamples = sample(J.copula, N)

    x = invCdf.(J.marginal1,copulaSamples[:,1]);
    y = invCdf.(J.marginal2,copulaSamples[:,2]);
    return hcat(x,y)

end

rand(J :: AbstractJoint, N :: Int64 = 1) = sample(J :: AbstractJoint, N)

function invCdf(X :: Sampleable{Univariate}, u)
    return quantile.(X,u)                           # We will also need the interpolartor here
end

function calcDensity(J :: AbstractJoint)

    #return calcDensity(J.cdf, J.xRange, J.yRange)

    if ismissing(J.copula.func) return calcDensity(J.cdf,J.xRange,J.yRange) end

    e = sqrt(sqrt(eps()));

    #x = J.xRange - e;f
    #y = J.yRange - e;

    x = range(J.xRange[1], stop = J.xRange[end]-e, length = length(J.xRange))
    y = range(J.xRange[1], stop = J.yRange[end]-e, length = length(J.yRange))

    der1 = ( J(x .+ e, y)       - J(x, y)     )/e     # df/dx
    der2 = ( J(x .+ e, y .+ e ) - J(x, y .+ e) )/e     # d(df/dx)/dy
    density = (der2-der1)/e                            #d^2f/dxdy

    return density

end

M(M1::ContinuousUnivariateDistribution, M2::ContinuousUnivariateDistribution) = M()(M1,M2)
W(M1::ContinuousUnivariateDistribution, M2::ContinuousUnivariateDistribution) = W()(M1,M2)
Pi(M1::ContinuousUnivariateDistribution, M2::ContinuousUnivariateDistribution) = Pi()(M1,M2)
Frank(M1::ContinuousUnivariateDistribution, M2::ContinuousUnivariateDistribution, s = 1) = Frank(s)(M1,M2)
Clayton(M1::ContinuousUnivariateDistribution, M2::ContinuousUnivariateDistribution, t = 0) = Clayton(t)(M1,M2)
Gaussian(M1::ContinuousUnivariateDistribution, M2::ContinuousUnivariateDistribution, corr = 0) = Gaussian(corr)(M1,M2)

struct marginal <: AbstractMarginal

    cdf :: Array{Float64,1}
    range
    cdfInp
    invCdfInp

    # Could have density and a calcDensity function for pdf

    function marginal(cdf :: Array{Float64,1}, range)

        if (length(cdf) !== length(range)) throw(ArgumentError("Cdf and range must be the same length")) end

        inpl = interpolate(cdf, BSpline(Quadratic(Line(OnCell()))));
        cdfInp = Interpolations.scale(inpl, range);

        # Error here and also in differention
        invCdfInp = LinearInterpolation(cdf, range, extrapolation_bc = Flat())

        new(cdf, range, cdfInp, invCdfInp)
    end
end

function (obj:: marginal)(x)
    return cdf(obj, x)
end

function cdf(D :: AbstractMarginal, x :: Real)

    if x < D.range[1] return 0; end
    if x > D.range[end] return 1; end

    return D.cdfInp(x)

end

function invCdf(D :: AbstractMarginal, u :: Real)

    return D.invCdfInp(u)

end

quantile(D :: marginal, u) = invCdf(D, u);

function sample(D :: marginal, N = 1)
    return invCdf(D, rand(N))
end

rand(D::marginal, N=1) = sample(D, N)

function Base.show(io::IO, z::copula)

    statement1 = "Arbitrary"
    statement2 = ""

    if (!ismissing(z.func))
        func = z.func
        if (func == indep); func ="π";end
        if (func == perf); func = "M";end
        if (func == opp); func = "W";end
        statement1 = "$(func)"
    end

    if (!ismissing(z.param))
        func = z.func
        parName = "par"
        if (func == Gau) parName = "r";end
        if (func == F) parName = "s";end
        if (func == Cla) parName = "t";end
        statement2 = "$parName=$(z.param)"
    end

    print(io, "Copula ~ $statement1($statement2)");
end

function Base.show(io::IO, z::Joint)

    statement1 = "Arbitrary"
    statement2 = ""

    if (!ismissing(z.copula.func))
        func = z.copula.func
        if (func == indep); func ="π";end
        if (func == perf); func = "M";end
        if (func == opp); func = "W";end
        statement1 = "$(func)"
    end

    if (!ismissing(z.copula.param))
        func = z.copula.func
        parName = "par"
        if (func == Gau) parName = "r";end
        if (func == F) parName = "s";end
        if (func == Cla) parName = "t";end
        if (func == spearLB) parName = "ρ";end
        if (func == kenLB) parName = "τ";end
        statement2 = "$parName=$(z.copula.param), "
    end

    print(io, "Joint ~ $statement1( $statement2$(z.marginal1), $(z.marginal2) )");
end

function Base.show(io::IO, z::marginal)

    print(io, "Distribution ~ ( range=[$(z.range[1]),$(z.range[end])])");
end

function linInterp(A, x)

    ind0 = findlast(A .<= x);
    x0 = A[ind0];   x1 = A[ind0+1]
    y0 = ind0/n;     y1 = (ind0+1)/n

    return y0 + (x-x0) * (y1-y0)/(x1-x0)
end


text:
import Pkg

# tries to load a Python dependency; on failure, adds the dependency via Conda
function check_add_dep(pkg; channel="")
    try
        @info "Checking dependency: $pkg"
        pyimport(pkg)
    catch
        @info "Installing $pkg..."
        Conda.add(pkg, channel=channel)
        pyimport(pkg)
        @info "$pkg install verified."
    end
end

@info "Setting up Python environment..."
ENV["PYTHON"]=""
Pkg.add("PyCall")
Pkg.build("PyCall")
using PyCall
pyimport("sys")
@info "PyCall verified."

# check for Conda; if not found, install it
try
    @info "Checking Conda."
    using Conda
    @info "Conda verified."
catch
    @info "Installing Conda..."
    Pkg.add("Conda")
    using Conda
    @info "Conda verified."
end

# check for Xtals; if not found, add it
try
    @info "Checking Xtals."
    using Xtals
    @info "Xtals verified."
catch # if not, install it
    @info "Installing Xtals..."
    Pkg.add("Xtals")
    Pkg.build("Xtals")
    using Xtals
    @info "Xtals install verified."
end

# check for XtalsPyTools; if not found, add it
try
    @info "Checking XtalsPyTools."
    using XtalsPyTools
    @info "XtalsPyTools verified."
catch # if not, install it
    @info "Installing XtalsPyTools..."
    Pkg.add("XtalsPyTools")
    Pkg.build("XtalsPyTools")
    using XtalsPyTools
    @info "XtalsPyTools install verified."
end

# check the deps, add if missing
check_add_dep("scipy")
check_add_dep("pymatgen", channel="conda-forge")

@info "Setup complete!"

max_stars_repo_name: SimonEnsemble/XtalsPyTools

max_issues_repo_name: SimonEnsemble/XtalsPyTools