text:
"""
    NestedSamplers.Proposals

This module contains the different algorithms for proposing new points within a bounding volume in unit space.

The available implementations are
* [`Proposals.Uniform`](@ref) - samples uniformly within the bounding volume
* [`Proposals.RWalk`](@ref) - random walks to a new point given an existing one
* [`Proposals.RStagger`](@ref) - random staggering away to a new point given an existing one
* [`Proposals.Slice`](@ref) - slicing away to a new point given an existing one
* [`Proposals.RSlice`](@ref) - random slicing away to a new point given an existing one
"""
module Proposals

using ..Bounds

using Random
using LinearAlgebra
using Parameters

export AbstractProposal

"""
    NestedSamplers.AbstractProposal

The abstract type for live point proposal algorithms.

# Interface

Each `AbstractProposal` must have this function,
```julia
(::AbstractProposal)(::AbstractRNG, point, loglstar, bounds, loglikelihood, prior_transform)
```
which, given the input `point` with loglikelihood `loglstar` inside a `bounds`, returns a new point in unit space, prior space, the loglikelihood, and the number of function calls.
"""
abstract type AbstractProposal end

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

# Helper for checking unit-space bounds
unitcheck(us) = all(u -> 0 < u < 1, us)

"""
    Proposals.Uniform()

Propose a new live point by uniformly sampling within the bounding volume.
"""
struct Uniform <: AbstractProposal end

function (::Uniform)(rng::AbstractRNG,
    point::AbstractVector,
    logl_star,
    bounds::AbstractBoundingSpace,
    loglike,
    prior_transform)

    ncall = 0
    while true
        u = rand(rng, bounds)
        unitcheck(u) || continue
        v = prior_transform(u)
        logl = loglike(v)
        ncall += 1
        logl ≥ logl_star && return u, v, logl, ncall
    end
end

Base.show(io::IO, p::Uniform) = print(io, "NestedSamplers.Proposals.Uniform")

"""
    Proposals.RWalk(;ratio=0.5, walks=25, scale=1)

Propose a new live point by random walking away from an existing live point.

## Parameters
- `ratio` is the target acceptance ratio
- `walks` is the minimum number of steps to take
- `scale` is the proposal distribution scale, which will update _between_ proposals.
"""
@with_kw mutable struct RWalk <: AbstractProposal
    ratio = 0.5
    walks = 25
    scale = 1.0

    @assert 1 / walks ≤ ratio ≤ 1 "Target acceptance ratio must be between 1/`walks` and 1"
    @assert walks > 1 "Number of steps must be greater than 1"
    @assert scale ≥ 0 "Proposal scale must be non-negative"
end

function (prop::RWalk)(rng::AbstractRNG,
                       point::AbstractVector,
                       logl_star,
                       bounds::AbstractBoundingSpace,
                       loglike,
                       prior_transform;
                       kwargs...)
    # setup
    n = length(point)
    scale_init = prop.scale
    accept = reject = fail = nfail = nc = ncall = 0
    local du, u_prop, logl_prop, u, v, logl

    while nc < prop.walks || iszero(accept)
        # get proposed point
        while true
            # check scale factor to avoid over-shrinking
            prop.scale < 1e-5 * scale_init && error("Random walk sampling appears to be stuck.")
            # transform to proposal distribution
            du = randoffset(rng, bounds)
            u_prop = @. point + prop.scale * du
            # inside unit-cube
            unitcheck(u_prop) && break

            fail += 1
            nfail += 1
            # check if stuck generating bad numbers
            if fail > 100 * prop.walks
                @warn "Random number generation appears extremely inefficient. Adjusting the scale-factor accordingly"
                fail = 0
                prop.scale *= exp(-1/n)
            end
        end
        # check proposed point
        v_prop = prior_transform(u_prop)
        logl_prop = loglike(v_prop)
        if logl_prop ≥ logl_star
            u = u_prop
            v = v_prop
            logl = logl_prop
            accept += 1
        else
            reject += 1
        end
        nc += 1
        ncall += 1

        # check if stuck generating bad points
        if nc > 50 * prop.walks
            @warn "Random walk proposals appear to be extremely inefficient. Adjusting the scale-factor accordingly"
            prop.scale *= exp(-1/n)
            nc = accept = reject = 0
        end
    end

    # update proposal scale using acceptance ratio
    update_scale!(prop, accept, reject, n)

    return u, v, logl, ncall
end

# update proposal scale using target acceptance ratio
function update_scale!(prop, accept, reject, n)
    ratio = accept / (accept + reject)
    norm = max(prop.ratio, 1 - prop.ratio) * n
    scale = prop.scale * exp((ratio - prop.ratio) / norm)
    prop.scale = min(scale, sqrt(n))
    return prop
end

"""
    Proposals.RStagger(;ratio=0.5, walks=25, scale=1)

Propose a new live point by random staggering away from an existing live point.
This differs from the random walk proposal in that the step size here is exponentially adjusted
to reach a target acceptance rate _during_ each proposal, in addition to _between_
proposals.

## Parameters
- `ratio` is the target acceptance ratio
- `walks` is the minimum number of steps to take
- `scale` is the proposal distribution scale, which will update _between_ proposals.
"""
@with_kw mutable struct RStagger <: AbstractProposal
    ratio = 0.5
    walks = 25
    scale = 1.0

    @assert 1 / walks ≤ ratio ≤ 1 "Target acceptance ratio must be between 1/`walks` and 1"
    @assert walks > 1 "Number of steps must be greater than 1"
    @assert scale ≥ 0 "Proposal scale must be non-negative"
end

function (prop::RStagger)(rng::AbstractRNG,
                          point::AbstractVector,
                          logl_star,
                          bounds::AbstractBoundingSpace,
                          loglike,
                          prior_transform;
                          kwargs...)
    #setup
    n = length(point)
    scale_init = prop.scale
    accept = reject = fail = nfail = nc = ncall = 0
    stagger = 1
    local du, u_prop, logl_prop, u, v, logl

    while nc < prop.walks || iszero(accept)
        # get proposed point
        while true
            # check scale factor to avoid over-shrinking
            prop.scale < 1e-5 * scale_init && error("Random walk sampling appears to be stuck.")
            # transform to proposal distribution
            du = randoffset(rng, bounds)
            u_prop = @. point + prop.scale * stagger * du
            # inside unit-cube
            unitcheck(u_prop) && break

            fail += 1
            nfail += 1
            # check if stuck generating bad numbers
            if fail > 100 * prop.walks
                @warn "Random number generation appears extremely inefficient. Adjusting the scale-factor accordingly"
                fail = 0
                prop.scale *= exp(-1/n)
            end
        end
        # check proposed point
        v_prop = prior_transform(u_prop)
        logl_prop = loglike(v_prop)
        if logl_prop ≥ logl_star
            u = u_prop
            v = v_prop
            logl = logl_prop
            accept += 1
        else
            reject += 1
        end
        nc += 1
        ncall += 1

        # adjust _stagger_ to target an acceptance ratio of `prop.ratio`
        ratio = accept / (accept + reject)
        if ratio > prop.ratio
            stagger *= exp(1 / accept)
        elseif ratio < prop.ratio
            stagger /= exp(1 / reject)
        end

        # check if stuck generating bad points
        if nc > 50 * prop.walks
            @warn "Random walk proposals appear to be extremely inefficient. Adjusting the scale-factor accordingly"
            prop.scale *= exp(-1 / n)
            nc = accept = reject = 0
        end
    end

    # update proposal scale using acceptance ratio
    update_scale!(prop, accept, reject, n)

    return u, v, logl, ncall
end

"""
    Proposals.Slice(;slices=5, scale=1)

Propose a new live point by a series of random slices away from an existing live point.
This is a standard _Gibbs-like_ implementation where a single multivariate slice is a combination of `slices` univariate slices through each axis.

## Parameters
- `slices` is the minimum number of slices
- `scale` is the proposal distribution scale, which will update _between_ proposals.
"""
@with_kw mutable struct Slice <: AbstractProposal
    slices = 5
    scale = 1.0

    @assert slices ≥ 1 "Number of slices must be greater than or equal to 1"
    @assert scale ≥ 0 "Proposal scale must be non-negative"
end

function (prop::Slice)(rng::AbstractRNG,
                       point::AbstractVector,
                       logl_star,
                       bounds::AbstractBoundingSpace,
                       loglike,
                       prior_transform;
                       kwargs...)
    # setup
    n = length(point)
    nc = nexpand = ncontract = 0
    local u, v, logl

    # modifying axes and computing lengths
    axes = Bounds.axes(bounds)
    axes = prop.scale .* axes'
    # slice sampling loop
    for it in 1:prop.slices

        # shuffle axis update order
        idxs = shuffle!(rng, collect(Base.axes(axes, 1)))

        # slice sample along a random direction
        for idx in idxs

            # select axis
            axis = axes[idx, :]

            u, v, logl, nc, nexpand, ncontract = sample_slice(rng, axis, point, logl_star, loglike,
                                                              prior_transform, nc, nexpand, ncontract)
        end # end of slice sample along a random direction
    end # end of slice sampling loop

    # update slice proposal scale based on the relative size of the slices compared to the initial guess
    prop.scale = prop.scale * nexpand / (2.0 * ncontract)

    return u, v, logl, nc
end   # end of function Slice

"""
    Proposals.RSlice(;slices=5, scale=1)
Propose a new live point by a series of random slices away from an existing live point.
This is a standard _random_ implementation where each slice is along a random direction based on the provided axes.
## Parameters
- `slices` is the minimum number of slices
- `scale` is the proposal distribution scale, which will update _between_ proposals.
"""
@with_kw mutable struct RSlice <: AbstractProposal
    slices = 5
    scale = 1.0

    @assert slices ≥ 1 "Number of slices must be greater than or equal to 1"
    @assert scale ≥ 0 "Proposal scale must be non-negative"
end

function (prop::RSlice)(rng::AbstractRNG,
                       point::AbstractVector,
                       logl_star,
                       bounds::AbstractBoundingSpace,
                       loglike,
                       prior_transform;
                       kwargs...)
    # setup
    n = length(point)
    nc = nexpand = ncontract = 0
    local u, v, logl
    # random slice sampling loop
    for it in 1:prop.slices

        # propose a direction on the unit n-sphere
        drhat = randn(rng, n)
        drhat /= norm(drhat)

        # transform and scale into parameter space
        axis = prop.scale .* (Bounds.axes(bounds) * drhat)
        u, v, logl, nc, nexpand, ncontract = sample_slice(rng, axis, point, logl_star, loglike,
                                                          prior_transform, nc, nexpand, ncontract)
    end # end of random slice sampling loop

    # update random slice proposal scale based on the relative size of the slices compared to the initial guess
    prop.scale = prop.scale * nexpand / (2.0 * ncontract)

    return u, v, logl, nc
end # end of function RSlice

# Method for slice sampling
function sample_slice(rng, axis, u, logl_star, loglike, prior_transform, nc, nexpand, ncontract)
    # define starting window
    r = rand(rng)  # initial scale/offset
    u_l = @. u - r * axis  # left bound
    if unitcheck(u_l)
        v_l = prior_transform(u_l)
        logl_l = loglike(v_l)
    else
        logl_l = -Inf
    end
    nc += 1
    nexpand += 1

    u_r = u_l .+ axis # right bound
    if unitcheck(u_r)
        v_r = prior_transform(u_r)
        logl_r = loglike(v_r)
    else
        logl_r = -Inf
    end
    nc += 1
    nexpand += 1

    # stepping out left and right bounds
    while logl_l ≥ logl_star
        u_l .-= axis
        if unitcheck(u_l)
            v_l = prior_transform(u_l)
            logl_l = loglike(v_l)
        else
            logl_l = -Inf
        end
        nc += 1
        nexpand += 1
    end

    while logl_r ≥ logl_star
        u_r .+= axis
        if unitcheck(u_r)
            v_r = prior_transform(u_r)
            logl_r = loglike(v_r)
        else
            logl_r = -Inf
        end
        nc += 1
        nexpand += 1
    end

    # sample within limits. If the sample is not valid, shrink the limits until the `logl_star` bound is hit
    window_init = norm(u_r - u_l)  # initial window size
    while true
        # define slice and window
        u_hat = u_r - u_l
        window = norm(u_hat)

        # check if the slice has shrunk to be ridiculously small
        window < 1e-5 * window_init && error("Slice sampling appears to be stuck.")

        # propose a new position
        r = rand(rng)
        u_prop = @. u_l + r * u_hat
        if unitcheck(u_prop)
            v_prop = prior_transform(u_prop)
            logl_prop = loglike(v_prop)
        else
            logl_prop = -Inf
        end
        nc += 1
        ncontract += 1

        # if success, then move to the new position
        if logl_prop ≥ logl_star
            return u_prop, v_prop, logl_prop, nc, nexpand, ncontract
        # if fail, then check if the new point is to the left/right of the original point along the proposal axis and update the bounds accordingly
        else
            s = dot(u_prop - u, u_hat) # check sign (+/-)
            if s < 0   # left
                u_l = u_prop
            elseif s > 0  # right
                u_r = u_prop
            else # if `s = 0` something has gone wrong
                error("Slice sampler has failed to find a valid point.")
            end
        end
    end # end of sample within limits while
end

end # module Proposals


text:
# Show different wasy to index associative arrays.

# Read CSV file into associative array.
E = ReadCSV(joinpath(Base.source_dir(),"Edge.csv"))
printFull(E)

# Get orange edges.
Eo = E[(E[:,"Color,"] == "Orange" ).row,:]
printFull(Eo)

# Get orange and green edges.
Eog = E[ StartsWith("O,G,") ,:]
printFull(Eog)

ext: jl


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# test meta-expressions that annotate blocks of code

const inlining_on = Base.JLOptions().can_inline != 0

function f(x)
    y = x+5
    z = y*y
    q = z/y
    m = q-3
end

@inline function f_inlined(x)
    y = x+5
    z = y*y
    q = z/y
    m = q-3
end

g(x) = f(2x)
g_inlined(x) = f_inlined(2x)

@test g(3) == g_inlined(3)
@test f(3) == f_inlined(3)

f() = backtrace()
@inline g_inlined() = f()
@noinline g_noinlined() = f()
h_inlined() = g_inlined()
h_noinlined() = g_noinlined()

function foundfunc(bt, funcname)
    for b in bt
        for lkup in StackTraces.lookup(b)
            if lkup.func == funcname
                return true
            end
        end
    end
    false
end
@test foundfunc(h_inlined(), :g_inlined)
@test foundfunc(h_noinlined(), :g_noinlined)

using Base: pushmeta!, popmeta!

macro attach(val, ex)
    esc(_attach(val, ex))
end

_attach(val, ex) = pushmeta!(ex, :test, val)

@attach 42 function dummy()
    false
end

asts = code_lowered(dummy, Tuple{})
@test length(asts) == 1
ast = asts[1]

body = Expr(:block)
body.args = ast.code

@test popmeta!(body, :test) == (true, [42])
@test popmeta!(body, :nonexistent) == (false, [])

# Simple popmeta!() tests
ex1 = quote
    $(Expr(:meta, :foo))
    x*x+1
end
@test popmeta!(ex1, :foo)[1]
@test !popmeta!(ex1, :foo)[1]
@test !popmeta!(ex1, :bar)[1]
@test !(popmeta!(:(x*x+1), :foo)[1])

# Find and pop meta information from general ast locations
multi_meta = quote
    $(Expr(:meta, :foo1))
    y = x
    $(Expr(:meta, :foo2, :foo3))
    begin
        $(Expr(:meta, :foo4, Expr(:foo5, 1, 2)))
    end
    x*x+1
end
@test popmeta!(deepcopy(multi_meta), :foo1) == (true, [])
@test popmeta!(deepcopy(multi_meta), :foo2) == (true, [])
@test popmeta!(deepcopy(multi_meta), :foo3) == (true, [])
@test popmeta!(deepcopy(multi_meta), :foo4) == (true, [])
@test popmeta!(deepcopy(multi_meta), :foo5) == (true, [1,2])
@test popmeta!(deepcopy(multi_meta), :bar)  == (false, [])

# Test that popmeta!() removes meta blocks entirely when they become empty.
for m in [:foo1, :foo2, :foo3, :foo4, :foo5]
    @test popmeta!(multi_meta, m)[1]
end
@test Base.findmeta(multi_meta.args)[1] == 0

# Test that pushmeta! can push across other macros,
# in the case multiple pushmeta!-based macros are combined

@attach 40 @attach 41 @attach 42 dummy_multi() = return nothing

asts = code_lowered(dummy_multi, Tuple{})
@test length(asts) == 1
ast = asts[1]

body = Expr(:block)
body.args = ast.code

@test popmeta!(body, :test) == (true, [40])
@test popmeta!(body, :test) == (true, [41])
@test popmeta!(body, :test) == (true, [42])
@test popmeta!(body, :nonexistent) == (false, [])

# tests to fully cover functions in base/meta.jl
using Base.Meta

@test isexpr(:(1+1),Set([:call]))
@test isexpr(:(1+1),Vector([:call]))
@test isexpr(:(1+1),(:call,))
@test isexpr(1,:call)==false
@test isexpr(:(1+1),:call,3)
ioB = IOBuffer()
show_sexpr(ioB,:(1+1))

show_sexpr(ioB,QuoteNode(1),1)

# test base/expr.jl
baremodule B
    eval = 0
    x = 1
    module M; x = 2; end
    import Base
    @Base.eval x = 3
    @Base.eval M x = 4
end
@test B.x == 3
@test B.M.x == 4

# specialization annotations

function _nospec_some_args(@nospecialize(x), y, @nospecialize z::Int)
end
@test first(methods(_nospec_some_args)).nospecialize == 5
@test first(methods(_nospec_some_args)).sig == Tuple{typeof(_nospec_some_args),Any,Any,Int}
function _nospec_some_args2(x, y, z)
    @nospecialize x y
    return 0
end
@test first(methods(_nospec_some_args2)).nospecialize == 3
function _nospec_with_default(@nospecialize x = 1)
    2x
end
@test collect(methods(_nospec_with_default))[2].nospecialize == 1
@test _nospec_with_default() == 2
@test _nospec_with_default(10) == 20

let oldout = STDOUT
    local rdout, wrout, out
    try
        rdout, wrout = redirect_stdout()
        out = @async read(rdout, String)

        @test eval(:(@dump x + y)) === nothing

        redirect_stdout(oldout)
        close(wrout)

        @test wait(out) == """
            Expr
              head: Symbol call
              args: Array{Any}((3,))
                1: Symbol +
                2: Symbol x
                3: Symbol y
              typ: Any
            """
    finally
        redirect_stdout(oldout)
    end
end

macro is_dollar_expr(ex)
    return Meta.isexpr(ex, :$)
end

module TestExpandModule
macro is_in_def_module()
    return __module__ === @__MODULE__
end
end

let a = 1
    @test @is_dollar_expr $a
    @test !TestExpandModule.@is_in_def_module
    @test @eval TestExpandModule @is_in_def_module

    @test Meta.lower(@__MODULE__, :($a)) === 1
    @test !Meta.lower(@__MODULE__, :(@is_dollar_expr $a))
    @test Meta.@lower @is_dollar_expr $a
    @test Meta.@lower @__MODULE__() @is_dollar_expr $a
    @test !Meta.@lower TestExpandModule.@is_in_def_module
    @test Meta.@lower TestExpandModule @is_in_def_module

    @test macroexpand(@__MODULE__, :($a)) === 1
    @test !macroexpand(@__MODULE__, :(@is_dollar_expr $a))
    @test @macroexpand @is_dollar_expr $a
end

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
import ..UncertainValues:
    AbstractScalarPopulation,
    UncertainScalarPopulation,
    UncertainValue

import Distributions
UVAL_TYPES = Union{T1, T2} where {T1 <: AbstractUncertainValue, T2 <: Distributions.Distribution}

"""
    truncate(population::AbstractScalarPopulation, constraint::NoConstraint)

Get the elements and the associated sampling weights of the `population` members
satisfying the sampling `constraint`.

- If `constraint` is a `NoConstraint` instance, then all members and weights are
    returned unmodified.
"""
function Base.truncate(pop::AbstractScalarPopulation, constraint::NoConstraint, n::Int = 30000)
    ConstrainedUncertainScalarPopulation(pop.values, pop.probs)
end

############################################################
# Populations whose members are strictly real-valued scalars
############################################################
function Base.truncate(pop::UncertainScalarPopulation{T, PW}, constraint::TruncateRange, n::Int = 30000) where {T <: Number, PW}
    inds = findall(constraint.min .<= pop.values .<= constraint.max)

    if length(inds) == 0
        #throw(ArgumentError("$pop could not be truncated. No values left after truncation."))
        return nothing
    end

    ConstrainedUncertainScalarPopulation(pop.values[inds], pop.probs[inds])
end

function Base.truncate(pop::UncertainScalarPopulation{T, PW}, constraint::TruncateMaximum, n::Int = 30000) where {T <: Number, PW}
    inds = findall(pop.values .<= constraint.max)

    if length(inds) == 0
        #throw(ArgumentError("$pop could not be truncated. No values left after truncation."))
        return nothing
    end

    ConstrainedUncertainScalarPopulation(pop.values[inds], pop.probs[inds])
end

function Base.truncate(pop::UncertainScalarPopulation{T, PW}, constraint::TruncateMinimum, n::Int = 30000) where {T <: Number, PW}
    inds = findall(constraint.min .<= pop.values)

    if length(inds) == 0
        #throw(ArgumentError("$pop could not be truncated. No values left after truncation."))
        return nothing
    end

    ConstrainedUncertainScalarPopulation(pop.values[inds], pop.probs[inds])
end

function Base.truncate(pop::UncertainScalarPopulation{T, PW}, constraint::TruncateLowerQuantile, n::Int = 30000) where {T <: Number, PW}
    lower_bound = quantile(pop, constraint.lower_quantile)

    inds = findall(lower_bound .<= pop.values)

    if length(inds) == 0
        #throw(ArgumentError("$pop could not be truncated. No values left after truncation."))
        return nothing
    end

    ConstrainedUncertainScalarPopulation(pop.values[inds], pop.probs[inds])
end

function Base.truncate(pop::UncertainScalarPopulation{T, PW}, constraint::TruncateUpperQuantile, n::Int = 30000) where {T <: Number, PW}
    upper_bound = quantile(pop, constraint.upper_quantile)

    inds = findall(pop.values .<= upper_bound)

    if length(inds) == 0
        #throw(ArgumentError("$pop could not be truncated. No values left after truncation."))
        return nothing
    end

    ConstrainedUncertainScalarPopulation(pop.values[inds], pop.probs[inds])
end

function Base.truncate(pop::UncertainScalarPopulation{T, PW}, constraint::TruncateQuantiles, n::Int = 30000) where {T <: Number, PW}
    lower_bound = quantile(pop, constraint.lower_quantile)
    upper_bound = quantile(pop, constraint.upper_quantile)

    inds = findall(lower_bound .<= pop.values .<= upper_bound)

    if length(inds) == 0
        #throw(ArgumentError("$pop could not be truncated. No values left after truncation."))
        return nothing
    end

    ConstrainedUncertainScalarPopulation(pop.values[inds], pop.probs[inds])
end

function Base.truncate(pop::UncertainScalarPopulation{T, PW}, constraint::TruncateStd, n::Int = 30000) where {T <: Number, PW}
    # Draw a sample of size `n` of the member of `p` according to their weights.
    s = rand(pop, n)

    # Compute mean and standard deviation
    p_mean = mean(s)
    p_stdev = std(s)
    nσ  = constraint.nσ
    inds = findall(p_mean - p_stdev*nσ .<= pop.values .<= p_mean + p_stdev*nσ)

    if length(inds) == 0
        throw(ArgumentError("$pop could not be truncated. No values left after truncation."))
        return nothing
    end

    ConstrainedUncertainScalarPopulation(pop.values[inds], pop.probs[inds])
end

############################################################
# Populations whose members are some sort of uncertain value
############################################################
TRUNCVAL_TYPES = Union{T1, T2} where {
    T1 <: AbstractUncertainValue,
    T2 <: Distributions.Distribution}

export TRUNCVAL_TYPES

function Base.truncate(pop::UncertainScalarPopulation{T, PW}, constraint::TruncateMinimum) where {
        T <: AbstractUncertainValue, PW}
    maxs = [maximum(uv) for uv in pop]

    # Find all distributions whose supports start *above or at* the minimum value imposed by the constraint.
    inds = findall(constraint.min .<= maxs)

    if length(inds) == 0
        throw(ArgumentError("$pop could not be truncated. No values left after truncation."))
    end

    # Constrain those distributions and match them with their respective probabilities
    truncated_vals = Vector{TRUNCVAL_TYPES}(undef, length(inds))

    for (i, val) in enumerate(pop[inds])
        truncated_vals[i] = constrain(val, constraint)
    end
    truncated_vals = Vector{TRUNCVAL_TYPES}(undef, 0)
    inds = Vector{Int}(undef, 0)

    for (i, val) in enumerate(pop.values)
        if maximum(val) >= constraint.min
            c = constrain(val, constraint)

            if !(c isa Nothing)
                push!(inds, i)
                push!(truncated_vals, c)
            end
        end
    end

    if length(inds) > 0
        ConstrainedUncertainScalarPopulation(truncated_vals, pop.probs[inds])
    else
        nothing
    end
end

function Base.truncate(pop::UncertainScalarPopulation{T, PW}, constraint::TruncateMinimum,
        n::Int) where {T <: AbstractUncertainValue, PW}
    Base.truncate(pop, constraint)
end

function Base.truncate(pop::UncertainScalarPopulation{T, PW}, constraint::TruncateMaximum) where {
        T <: AbstractUncertainValue, PW}

    mins = [minimum(uv) for uv in pop]

    # Find all distributions whose supports start *below or at* the maximum value imposed by the constraint.
    inds = findall(mins .<= constraint.max)

    if length(inds) == 0
        #throw(ArgumentError("$pop could not be truncated. No values left after truncation."))
        return nothing
    end

    # Constrain those distributions and match them with their respective probabilities
    truncated_vals = Vector{TRUNCVAL_TYPES}(undef, 0)
    inds = Vector{Int}(undef, 0)

    for (i, val) in enumerate(pop.values)
        if minimum(val) < constraint.max
            c = constrain(val, constraint)

            if !(c isa Nothing)
                push!(inds, i)
                push!(truncated_vals, c)
            end
        end
    end

    if length(inds) > 0
        return ConstrainedUncertainScalarPopulation(truncated_vals, pop.probs[inds])
    else
        return nothing
    end
end

function Base.truncate(pop::UncertainScalarPopulation{T, PW}, constraint::TruncateMaximum,
        n::Int) where {T <: AbstractUncertainValue, PW}
    Base.truncate(pop, constraint)
end

function Base.truncate(pop::UncertainScalarPopulation{T, PW}, constraint::TruncateRange) where {
        T <: AbstractUncertainValue, PW}

    mins = [minimum(uv) for uv in pop]
    maxs = [maximum(uv) for uv in pop]

    # Find all distributions whose supports start *above or at* the minimum value imposed by the constraint.
    satisfies_minconstraint = constraint.min .<= maxs

    # Find all distributions whose supports start *below or at* the maximum value imposed by the constraint.
    satisfies_maxconstraint = mins .<= constraint.max

    # Find all distributions that satisfy both the lower constraint and the upper constraint
    inds = findall(satisfies_minconstraint .& satisfies_maxconstraint)

    if length(inds) == 0
        #throw(ArgumentError("$pop could not be truncated. No values left after truncation."))
        return nothing
    end

    # Constrain those distributions and match them with their respective probabilities
    truncated_vals = Vector{TRUNCVAL_TYPES}(undef, length(inds))

    for (i, val) in enumerate(pop[inds])
        truncated_vals[i] = constrain(val, constraint)
    end

    ConstrainedUncertainScalarPopulation(truncated_vals, pop.probs[inds])
end

function Base.truncate(pop::UncertainScalarPopulation{T, PW}, constraint::TruncateRange,
        n::Int) where {T <: AbstractUncertainValue, PW}
    Base.truncate(pop, constraint)
end

function Base.truncate(pop::UncertainScalarPopulation{T, PW},
        constraint::TruncateLowerQuantile,
        n::Int = 30000) where {T <: AbstractUncertainValue, PW}

    # We want to truncate each of the population members so that their furnishing distributions
    # cannot yield numbers smaller than the overall quantile of the population. First
    # find the overall lower quantile
    population_lower_quantile = quantile(pop, constraint.lower_quantile, n)

    # Now, truncate each of the population members below at the overall population lower
    # quantile. Probabilities are kept the same.
    truncated_vals = Vector{TRUNCVAL_TYPES}(undef, 0)
    inds = Vector{Int}(undef, 0)

    for (i, val) in enumerate(pop.values)
        if maximum(val) > population_lower_quantile
            push!(inds, i)
            push!(truncated_vals, truncate(val, TruncateMinimum(population_lower_quantile)))
        end
    end

    #println("There were $(length(inds)) population members left after truncation")

    if length(inds) > 0
        return ConstrainedUncertainScalarPopulation(truncated_vals, pop.probs[inds])
    else
        return nothing
    end
end

function Base.truncate(pop::UncertainScalarPopulation{T, PW},
        constraint::TruncateUpperQuantile,
        n::Int = 30000) where {T <: AbstractUncertainValue, PW}

    # We want to truncate each of the population members so that their furnishing distributions
    # cannot yield numbers larger than the overall quantile of the population. First
    # find the overall lower quantile
    population_upper_quantile = quantile(pop, constraint.upper_quantile, n)
    #@show "Overall quantile", population_upper_quantile

    # Now, truncate each of the population members above at the overall population upper
    # quantile. Probabilities are kept the same. We initialise an empty array, because
    # some values of the population may be dropped during the truncation process.
    truncated_vals = Vector{TRUNCVAL_TYPES}(undef, 0)
    inds = Vector{Int}(undef, 0)

    for (i, val) in enumerate(pop.values)
        if minimum(val) < population_upper_quantile
            push!(inds, i)
            push!(truncated_vals, truncate(val, TruncateMaximum(population_upper_quantile)))
        end
    end
    #println("There were $(length(inds)) population members left after truncation")

    if length(inds) > 0
        return ConstrainedUncertainScalarPopulation(truncated_vals, pop.probs[inds])
    else
        return nothing
    end
end

function Base.truncate(pop::UncertainScalarPopulation{T, PW},
        constraint::TruncateQuantiles,
        n::Int = 30000) where {T <: AbstractUncertainValue, PW}

    # We want to truncate each of the population members so that their furnishing distributions
    # cannot yield numbers larger than the overall upper quantile of the population, nor
    # numbers smaller than the overall lower quantile of the population. Find these.
    population_upper_quantile = quantile(pop, constraint.upper_quantile, n)
    population_lower_quantile = quantile(pop, constraint.lower_quantile, n)

    # Now, truncate each of the population members at the range given by the overall quantiles
    # Probabilities are kept the same.
    truncated_vals = Vector{TRUNCVAL_TYPES}(undef, 0)
    inds = Vector{Int}(undef, 0)

    for (i, val) in enumerate(pop.values)
        if maximum(val) > population_lower_quantile && minimum(val) < population_upper_quantile
            push!(inds, i)
            push!(truncated_vals, truncate(val, TruncateRange(population_lower_quantile, population_upper_quantile)))
        end
    end

    if length(inds) > 0
        return ConstrainedUncertainScalarPopulation(truncated_vals, pop.probs[inds])
    else
        return nothing
    end
end

function Base.truncate(pop::UncertainScalarPopulation{T, PW},
        constraint::TruncateStd{TN},
        n::Int = 30000) where {T <: AbstractUncertainValue, PW, TN <: Number}

    # We want to truncate each of the population members so that their furnishing distributions
    # cannot yield numbers outside +- constraint.nσ*population_std.
    s = rand(pop, n)
    population_std = std(s)
    population_mean = mean(s)

    upper_bound = population_mean + (constraint.nσ*population_std)
    lower_bound = population_mean - (constraint.nσ*population_std)

    # Now, truncate each of the population members at the range given by the overall quantiles
    # Probabilities are kept the same.

    truncated_vals = Vector{TRUNCVAL_TYPES}(undef, 0)
    inds = Vector{Int}(undef, 0)

    for (i, val) in enumerate(pop.values)
        if maximum(val) > lower_bound && minimum(val) < upper_bound
            push!(inds, i)
            push!(truncated_vals, truncate(val, TruncateRange(lower_bound, upper_bound)))
        end
    end

    if length(inds) > 0
        return ConstrainedUncertainScalarPopulation(truncated_vals, pop.probs[inds])
    else
        return nothing
    end

end

export truncate


text:
module ColorEdits

Base.Experimental.@optlevel 1

using Redux
using CImGui
import ..ReduxImGui: get_label, get_color

# actions
abstract type AbstractColorEditAction <: AbstractSyncAction end

"""
    Rename(label, new_label)
Change label to `new_label`.
Note that renaming may also change the identifier, please refer to `help?> CImGui.PushID`.
"""
struct Rename <: AbstractColorEditAction
    label::String
    new_label::String
end

"""
    EnableAlpha(label)
Enable editing alpha. This also changes widget appearance.
"""
struct EnableAlpha <: AbstractColorEditAction
    label::String
end

"""
    DisableAlpha(label)
Disable editing alpha. This also changes widget appearance.
"""
struct DisableAlpha <: AbstractColorEditAction
    label::String
end

"""
    ChangeFlags(label, flags)
Change widget's flags to `flags`.
"""
struct ChangeFlags <: AbstractColorEditAction
    label::String
    flags::CImGui.ImGuiColorEditFlags_
end

# state
struct State <: AbstractImmutableState
    label::String
    rgba::Vector{Cfloat}
    flags::CImGui.ImGuiColorEditFlags_
    no_alpha::Bool
end
State(label::AbstractString, rgba::Vector{Cfloat}, flags = 0) =
    State(label, rgba, CImGui.ImGuiColorEditFlags_(flags), true)

# reducers
color_edit(state::AbstractState, action::AbstractAction) = state
color_edit(state::Vector{<:AbstractState}, action::AbstractAction) = state
color_edit(state::Dict{String,<:AbstractState}, action::AbstractAction) = state

function color_edit(state::Dict{String,State}, action::AbstractColorEditAction)
    s = Dict{String,State}()
    for (k, v) in state
        s[k] = get_label(v) == action.label ? color_edit(v, action) : v
    end
    return s
end

color_edit(s::State, a::Rename) = State(a.new_label, s.rgba, s.flags, s.no_alpha)
color_edit(s::State, a::EnableAlpha) = State(s.label, s.rgba, s.flags, false)
color_edit(s::State, a::DisableAlpha) = State(s.label, s.rgba, s.flags, true)
color_edit(s::State, a::ChangeFlags) = State(s.label, s.rgba, a.flags, s.no_alpha)

# helper
"""
    ColorEdit(store::AbstractStore, get_state=Redux.get_state) -> Bool
Return `true` when triggered.
"""
function ColorEdit(store::AbstractStore, get_state = Redux.get_state)
    s = get_state(store)
    if s.no_alpha
        is_triggered = CImGui.ColorEdit3(get_label(s), s.rgba, s.flags)
    else
        is_triggered = CImGui.ColorEdit4(get_label(s), s.rgba, s.flags)
    end
    return is_triggered
end

get_label(s::State) = s.label
get_color(s::State) = s.rgba

end # module


text:
using NumericalMethodsforEngineers, DataFrames, Plots
pyplot(size=(700,700))

ProjDir = dirname(@__FILE__)
cd(ProjDir) #do

  x = [1.0, 3.0, 6.0, 5.0]
  y = [1.0, 5.0, 10.0, 9.0]
  xi = [2.0, 4.5]
  (dfin, dfxi) = lagrangianpolynomial(length(x), x, y, xi)

  xint = 1:0.1:5
  (dfin, dfxint) = lagrangianpolynomial(length(x), x, y, collect(xint))

  dfin |> display
  println()

  dfxi |> display
  println()

  dfxint |> display
  println()

  p = plot(dfxint[:xi], dfxint[:yi], line=(:path, 1), label="interpolated curve")
  scatter!(p, dfin[:x], dfin[:y], marker=(:circle, 4), label="input points", color=:blue)
  scatter!(p, dfxi[:xi], dfxi[:yi], marker=(:star, 8), color=:red, label="interpolated points")

  plot(p)
  savefig("ex.5.1.png")
  gui()

#end


text:
# https://diffeqflux.sciml.ai/dev/examples/LV-Univ/

using DiffEqFlux, Flux, Optim, OrdinaryDiffEq, Plots
using DataInterpolations

u0 = 1.1f0
tspan = (0.0f0, 25.0f0)
tsteps = 0.0f0:1.0:25.0f0

#
noise_itp = LinearInterpolation( 0.001*randn(Float32, length(tsteps)), tsteps )

model_univ = FastChain(FastDense(2, 16, tanh),
                       FastDense(16, 16, tanh),
                       FastDense(16, 1))

# The model weights are destructured into a vector of parameters
p_model = initial_params(model_univ)
n_weights = length(p_model)

# Parameters of the second equation (linear dynamics)
p_system = Float32[0.5, -0.5]

p_all = [p_model; p_system]

function dudt_univ!(du, u, p, t)
    # Destructure the parameters
    model_weights = p[1:n_weights]
    α = p[end - 1]
    β = p[end]

    # The neural network outputs a control taken by the system
    # The system then produces an output
    model_control, system_output = u

    # adding some noise to the system output
    system_output += noise_itp(t)

    # Dynamics of the control and system
    dmodel_control = model_univ(u, model_weights)[1]
    dsystem_output = α*system_output + β*model_control

    # Update in place
    du[1] = dmodel_control
    du[2] = dsystem_output
end

prob_univ = ODEProblem(dudt_univ!, u0, tspan, p_all)
sol_univ = concrete_solve(prob_univ, Tsit5(),[0f0, u0], p_all,
                          abstol = 1e-8, reltol = 1e-6)

function predict_univ(θ)
  return Array(concrete_solve(prob_univ, Tsit5(), [0f0, θ[1]], θ[2:end],
                              saveat = tsteps))
end

loss_univ(θ) = sum(abs2, predict_univ(θ)[2,:] .- 1)

θ = Float32[u0; p_all]
l = loss_univ(θ)

list_plots = []
iter = 0
callback = function (θ, l)
  global list_plots, iter

  if iter == 0
    list_plots = []
  end
  iter += 1

  println(l)

  plt = Plots.plot(predict_univ(θ), ylim = (0, 6))
  push!(list_plots, plt)
  display(plt)
  return false
end

result_univ = DiffEqFlux.sciml_train(loss_univ, θ,
                                     BFGS(initial_stepnorm = 0.01),
                                     cb = callback)
#

#


text:
# CUDA CuContext

immutable CuContext
    handle::Ptr{Void}
end

const CTX_SCHED_AUTO  = 0x00
const CTX_SCHED_SPIN  = 0x01
const CTX_SCHED_YIELD = 0x02
const CTX_SCHED_BLOCKING_SYNC = 0x04
const CTX_MAP_HOST = 0x08
const CTX_LMEM_RESIZE_TO_MAX = 0x10

function create_context(dev::CuDevice, flags::Integer)
    a = Array(Ptr{Void}, 1)
    @cucall(cuCtxCreate, (Ptr{Ptr{Void}}, Cuint, Cint), a, flags, dev.handle)
    return CuContext(a[1])
end

create_context(dev::CuDevice) = create_context(dev, 0)

function destroy(ctx::CuContext)
    @cucall(cuCtxDestroy, (Ptr{Void},), ctx.handle)
end

function push(ctx::CuContext)
    @cucall(cuCtxPushCurrent, (Ptr{Void},), ctx.handle)
end

function pop(ctx::CuContext)
    a = Array(Ptr{Void}, 1)
    @cucall(cuCtxPopCurrent, (Ptr{Ptr{Void}},), a)
    return CuContext(a[1])
end


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ bef0f870-eb91-11ea-02b6-dd0f1a40947c
md"# Introduction"

# ╔═╡ e6959750-eb91-11ea-3aea-519bdec8f9f6
md"""## Who are we?

Lecturer: LCL IMM Ben Lauwens / D30.20 / [ben.lauwens@mil.be]()

Assistant: CPN Bart De Clerck / D30.20 / [bart.declerck@mil.be]()"""

# ╔═╡ 0f6a2ab2-eb92-11ea-2a06-0d9b69dd1358
md"""## Why Modelling and Simulation

- What is modelling?

- What is simulation?

Reality is often too complex to calculate ..."""

# ╔═╡ 29d92270-eb92-11ea-374f-45b11865dd9a
md"""## Documentations

All slides can be found on github: [https://github.com/BenLauwens/ES313.jl.git]()."""

# ╔═╡ ee8e57c0-eb92-11ea-193a-5bfcf7800d99
md"""## Schedule

### Theory

- 31/08: Cellular Automaton + Game of Life
- 01/09: Year Coord
- 14/09: Physical Modelling + Self-Organization
- 21/09: Optimisation Techniques
- 22/09: Linear Programming I
- 28/09: Linear Programming II + Applications I
- 29/09: Applications of Linear Programming II
- 19/10: Introduction to Discrete Event Simulation
- 20/10: Process Driven DES: SimJulia I
- 26/10: Process Driven DES II + Applications I
- 27/10: Applications with SimJulia II

### Practice

- 07/09: Visualisation
- 08/09: Cellular Automaton
- 15/09: Game of Life
- 05/10: Physical Modelling + Self-Organization
- 06/10: Optimisation Techniques I
- 12/10: Optimisation Techniques II + Linear Programming I
- 13/10: Linear Programming II
- 09/11: Introduction to Discrete Event Simulation
- 10/11: Process Driven DES: SimJulia I
- 16/11: Process Driven DES II + Applications with SimJulia I
- 17/11: Applications with SimJulia II
- 24/11: Performance

### Project

- 09/11: List of projects available
- we are available during contact hours
- 23/11: obligatory meeting: understanding of the problem
- 07/12: obligatory meeting: progress"""

# ╔═╡ 0b3e8cf2-eb93-11ea-140d-35c9f1074ea9
md"""## Evaluation

Test: Oct - 2Hr
- Cellular Automaton + Game of Life
- Physical Modelling + Self-Organization

Examen: Project with Oral Defense
- Visualisation
- Optimisation Techniques
- Linear Programming
- Discrete Event Simulation"""

# ╔═╡ 5fe3a710-ebbe-11ea-3987-3fec2c7f3057
md"""## Julia

- Install Julia on CDN laptop
- Install Notepad++ on CDN laptop
- Start Julia and install Pluto:

```julia
using Pkg
pkg"add Pluto"
pkg"add PlutoUI"
```

- Start Pluto:

```julia
using Pluto
Pluto.run()
```
"""

# ╔═╡ Cell order:
# ╟─bef0f870-eb91-11ea-02b6-dd0f1a40947c
# ╟─e6959750-eb91-11ea-3aea-519bdec8f9f6
# ╟─0f6a2ab2-eb92-11ea-2a06-0d9b69dd1358
# ╟─29d92270-eb92-11ea-374f-45b11865dd9a
# ╟─ee8e57c0-eb92-11ea-193a-5bfcf7800d99
# ╟─0b3e8cf2-eb93-11ea-140d-35c9f1074ea9
# ╟─5fe3a710-ebbe-11ea-3987-3fec2c7f3057


text:
# using Distributions, Bijectors

# Define a deterministic distribution
struct Determin{T<:Real} <: ContinuousUnivariateDistribution
  val::T
end

Distributions.rand(rng::AbstractRNG, d::Determin) = d.val
Distributions.logpdf(d::Determin, x::T) where T<:Real = zero(x)

Bijectors.bijector(d::Determin) = Identity{0}() # <= `0` indicates 0-dimensional, i.e. univariate

# # useful functions
function compute_stats(chain::Chains, name::Symbol)
    # x = Array(x)
    # x = DataFrame(x)
    x = Array(group(chain, name))
    mean_x = mean(x, dims=1)
    xlb = mean_x - mapslices(x -> quantile(x, 0.025), x, dims=1)
    xub = mean_x - mapslices(x -> quantile(x, 0.975), x, dims=1)
    return (μ = mean_x, lb = xlb, ub = xub)
end

function compute_stats(x::Array)
    mean_x = mean(x, dims=1)
    xlb = mean_x - mapslices(x -> quantile(x, 0.025), x, dims=1)
    xub = mean_x - mapslices(x -> quantile(x, 0.975), x, dims=1)
    return (μ = mean_x, lb = xlb, ub = xub)
end

function get_omega(chain::Chains)
    w1 = Array(group(chain, :ω1))
    w2 = Array(group(chain, :ω2))
    w = w1 .+ w2
    return (w = w, w1 = w1, w2 = w2)
end

function summarize_omega(chain::Chains)
    w1_array = Array(group(chain, :ω1))
    w2_array = Array(group(chain, :ω2))
    w_array = Array(group(chain, :ω))

    w = compute_stats(w_array)
    w1 = compute_stats(w1_array)
    w2 = compute_stats(w2_array)
    return (w = w, w1 = w1, w2 = w2)
end

invlogit(x) = 1/(1 + exp(-x))

function beta_params_mode(mode, k)
    #= http://doingbayesiandataanalysis.blogspot.com/2012/06/beta-distribution-parameterized-by-mode.html
    mode = (α - 1)/(α + β - 2)
    k = a + b # sample size
    =#
    @assert 0 < mode < 1
    @assert k > 2

    α = mode * (k-2) + 1
    β = (1-mode) * (k-2) + 1
    return α, β
end

function beta_params(μ, σ)
    @assert 0 < μ < 1
    @assert σ > 0

    α = ((1-μ)/σ^2 - (1/μ)) * μ^2
    β = α * (1 / μ - 1)
    return α, β
end


text:
module NamedIndexing

export NamedAxisArray
export axisnames

struct NamedAxisArray{T, N, A<:AbstractArray{T, N}, Names} <: AbstractArray{T, N}
    data::A
end

function NamedAxisArray(data::A, names::NTuple{N, Symbol}) where {T, N, A<:AbstractArray{T, N}}
    NamedAxisArray{T, N, A, names}(data)
end

Base.size(AA::NamedAxisArray) = size(AA.data)

function Base.getindex(AA::NamedAxisArray{T, N, A, Names}, axinds::IndAx) where {T, N, A, Names, IndAx<:NamedTuple}
    fullinds = merge(
        NamedTuple{Names, NTuple{N, Colon}}(ntuple((_) -> (:), Val{N}())),
        axinds
    )

    fullkeys = keys(fullinds)
    if fullkeys != Names
        throw(ArgumentError("Unexpected named indexes $(setdiff(fullkeys, Names))"))
    end

    newdata = getindex(AA.data, fullinds...)

    # falls apart if you do APL-style reshaping with indexes
    NamedAxisArray{T, N, A, Names}(newdata)
end

Base.getindex(AA::NamedAxisArray; kwargs...) = getindex(AA, kwargs.data)

axisnames(AA::NamedAxisArray{T, N, A, Names}) where {T, N, A, Names} = Names

end # module


text:
using BattleshipSimulator

gr()

# Create a random 10x10 board with the classic ship lengths
ship_lengths = [5,4,3,3,2]
M = layoutBoard(10,10; ships=ship_lengths)

# Create the distribution of all possible ship placements
D = distributionBoard(zeros(Int64,size(M)); ships=ship_lengths)

# Visualize the board and the distribution of possible ships
drawBoardComp(M,D; name1 = "True Hidden Map", name2="Distribution of Possible Ships")
savefig("./Random_Board_Layout_and_Ship_Distribution.png")

# Play the board once with a random hunt strategy
hits_rand = randomHunter(M)

# Animate the random game
anim = animateGame(M,hits_rand)
gif(anim,"./anim_random.gif", fps = 2) # save animation to gif

# Play with the search and destroy strategy
hits_target = targetHunter(M)

# Animate the random game
anim = animateGame(M,hits_target)
gif(anim,"./anim_target.gif", fps = 2) # save animation to gif

# Play with the distribution hunter
hits_dist = distHunter(M)

# Animate the random game
anim = animateGame(M,hits_dist)
gif(anim,"./anim_dist.gif", fps = 2) # save animation to gif

# Play the same board 100_000 times with the random hunt strategy (about 5s on my machine)
@time hits_many_rand = playManySingleGames(M, 100_000; hunt=randomHunter)

# Check out the distribution of game lengths
histogram(hits_many_rand, normed=true, xlabel="Hits Needed to Win",ylabel="Probability", legend=(:none),xlim=(0,length(M)))
savefig("./dist_hits_many_rand.png")

# Play the same board 10_000 times with the random hunt strategy (about 5s on my machine)
@time hits_many_target = playManySingleGames(M, 100_000; hunt=targetHunter)

# Check out the distribution of game lengths
histogram(hits_many_target, normed=true, xlabel="Hits Needed to Win",ylabel="Probability", legend=(:none),xlim=(0,length(M)))
savefig("./dist_hits_many_target.png")

max_stars_repo_name: mstobb/BattleshipSimulator

max_issues_repo_name: mstobb/BattleshipSimulator

max_forks_repo_name: mstobb/BattleshipSimulator


text:
"""
    bipartite_map(g) -> Vector{UInt8}

For a bipartite graph `g`, return a vector `c` of size ``|V|`` containing
the assignment of each vertex to one of the two sets (``c_i == 1`` or ``c_i == 2``).
If `g` is not bipartite, return an empty vector.

### Implementation Notes
Note that an empty vector does not necessarily indicate non-bipartiteness.
An empty graph will return an empty vector but is bipartite.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = SimpleGraph(3);

julia> bipartite_map(g)
3-element Array{UInt8,1}:
 0x01
 0x01
 0x01

julia> add_vertices!(g, 3);

julia> add_edge!(g, 1, 2);

julia> add_edge!(g, 2, 3);

julia> bipartite_map(g)
3-element Array{UInt8,1}:
 0x01
 0x02
 0x01
```
"""
function bipartite_map(g::AbstractGraph{T}) where T
    nvg = nv(g)
    if !is_directed(g)
        ccs = filter(x -> length(x) >= 2, connected_components(g))
    else
        ccs = filter(x -> length(x) >= 2, weakly_connected_components(g))
    end
    seen = zeros(Bool, nvg)
    colors = zeros(Bool, nvg)
    for cc in ccs
        Q = Vector{T}()
        s = cc[1]
        push!(Q, s)
        bipartitemap = zeros(UInt8, nvg)
        while !isempty(Q)
            u = popfirst!(Q)
            for v in outneighbors(g, u)
                if !seen[v]
                    colors[v] = !colors[u]
                    push!(Q, v)
                    seen[v] = true
                elseif colors[v] == colors[u]
                    return Vector{UInt8}()
                end
            end
        end
    end
    return UInt8.(colors).+(one(UInt8))
end

"""
    is_bipartite(g)

Return `true` if graph `g` is [bipartite](https://en.wikipedia.org/wiki/Bipartite_graph).

# Examples
```jldoctest
julia> using LightGraphs

julia> g = SimpleGraph(3);

julia> add_edge!(g, 1, 2);

julia> add_edge!(g, 2, 3);

julia> is_bipartite(g)
true

julia> add_edge!(g, 1, 3);

julia> is_bipartite(g)
false
```
"""
is_bipartite(g::AbstractGraph) = length(bipartite_map(g)) == nv(g)


text:
using LightGraphs
using DataStructures

function readinput(strs)

    function parse_contents(s)::Dict{String,Int}
        contents = Dict{String,Int}()
        bs = split(s, ",")
        if bs[1] == "no other bags"
            return contents
        else
            for b in bs
                m = match(r".*(\d+)\s(.+) bag.*", b)
                contents[m[2]] = parse(Int, m[1])
            end
        end
        contents
    end

    function parse_rules!(rules, vertex, s)
        m = match(r"(.+) bags contain (.+).", s)
        rules[m[1]] = (vertex = vertex, contents = parse_contents(m[2]))
        m[1]
    end

    rules = Dict{String,@NamedTuple{vertex::Int, contents::Dict{String,Int64}}}()
    vtxclr = Vector{String}(undef, length(strs))

    for (i, s) in enumerate(strs)
        vtxclr[i] = parse_rules!(rules, i, s)
    end

    rulesgrf = SimpleDiGraph(length(strs))

    for (_, r) in rules
        rv = r.vertex
        cs = r.contents
        for (color, _) in pairs(r.contents)
            add_edge!(rulesgrf, rv, rules[color].vertex)
        end
    end

    return rules, vtxclr, rulesgrf
end

# First Half

function howmany_contains(rules, vtxclr, rulesgrf, color)

    function add_predecessors!(v, preds)

        inns = inneighbors(rulesgrf, v)

        if isempty(inns)
            return
        end
        push!(preds, inns...)

        for nv in inns
            add_predecessors!(nv, preds)
        end
        return
    end

    v = rules[color].vertex
    containers = Set()
    add_predecessors!(v, containers)
    length(containers)
end

puzzleinput = readlines("inputs/day07-test.txt")
(rules, vtxclr, rulesgrf) = readinput(puzzleinput)

howmany_contains(rules, vtxclr, rulesgrf, "shiny gold")

# Second Half

function howmany_inside(rules, color)

    function howmany_inside(color)

        contents = rules[color].contents

        if isempty(contents)
            return 0
        end
        count = sum(counter(contents))

        for (c, q) in contents
            count += q * howmany_inside(c)
        end
        return count
    end

    howmany_inside(color)
end

puzzleinput = readlines("inputs/day07.txt")
(rules, vtxclr, rulesgrf) = readinput(puzzleinput)

howmany_inside(rules, "shiny gold")


text:
import Symbolics
using Test

@test !(@isdefined Num)
Symbolics.@variables t
Symbolics.@register fff(t)
@test isequal(fff(t), Symbolics.Num(Symbolics.Term{Real}(fff, [Symbolics.value(t)])))


text:
export NKpLandscape

# Construct an NKp landscape, see:
# Barnett, L. (1998). Ruggedness and Neutrality - The NKp Family of Fitness
# Landscapes. In Artificial Life VI (pp. 18–27).
type NKpLandscape <: Landscape
  n::Int64
  k::Int64
  p::Float64
  a::Int64
  links::AlleleLinks
  contribs::AlleleContribs
end

function NKpLandscape(n::Int64, k::Int64, p::Float64; a::Int64=2, near::Bool=false)
  if a != 2
    # TODO:  implement a values other than 2
    error("only a == 2 is implemented at this time")
  end
  links = makelinks(n, k, near)
  contribs = newAlleleContribs(n)
  return NKpLandscape(n, k, p, a, links, contribs)
end

show(io::Base.IO, l::NKpLandscape) = print(io, "NKpLandscape($(l.n), $(l.k), $(l.p))")

# Compute an estimate of the probability that a mutation within an NKp
# landscape is neutral.
# Barnett, L. (1998). Ruggedness and Neutrality - The NKp Family of Fitness
# Landscapes. In Artificial Life VI (pp. 18–27).
function prob_neutral(l::NKpLandscape)
  l.p ^ 2 * (1 - l.k / (l.n - 1) * (1 - l.p ^ 2)) ^ (l.n - 1)
end


text:
using Documenter
using DayCounts

makedocs(
    pages = [
        "index.md",
        "api.md",
        ],
    sitename = "DayCounts",
    format = Documenter.HTML(),
    modules = [DayCounts]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/JuliaFinance/DayCounts.jl.git",
    deploy_config = Documenter.GitHubActions()
)


text:
import Base: +, length, size, ==, isempty
import ViewerGL
GL = ViewerGL

#-------------------------------------------------------------------------------
#   MODEL STRUCT DEFINITION
#-------------------------------------------------------------------------------

mutable struct Model
    G::Lar.Points
    T::Array{Lar.ChainOp, 1}

    function Model(V::Lar.Points, T::Array{Lar.ChainOp, 1})
        dim, npts = size(V)
        dim > 0 ||
            throw(ArgumentError("At least one point is needed."))
        length(T) == dim ||
            throw(ArgumentError("Topology is not coherent with Geometry."))

        isempty(T[1]) || size(T[1], 2) == npts ||
            throw(ArgumentError("Topology not coherent with Geometry."))
        for i = 2 : dim
            isempty(T[i-1]) || isempty(T[i]) ||
                size(T[i], 2) == size(T[i-1],1) ||
                throw(ArgumentError("Topology not coherent with Topology."))
        end

        new(V, T)
    end

    function Model(V::Lar.Points)
        T = convert(Array{Lar.ChainOp,1},
            [SparseArrays.spzeros(Int8, 0, 0) for i = 1 : size(V, 1)]
        )
        T[1] = convert(Lar.ChainOp, SparseArrays.spzeros(Int8, 0, size(V,2)))
        Lar.Model(V, T)
    end

    function Model()
        nothing
    end
end

#-------------------------------------------------------------------------------
#   BASIC PROPERTIES
#-------------------------------------------------------------------------------

length(m::Lar.Model)               = size(m.G, 1)
size(m::Lar.Model)		           = size(m.G)
size(m::Lar.Model, i::Int)         = i == 0 ? size(m.G   ) : size(m.T[i]   )
size(m::Lar.Model, i::Int, j::Int) = i == 0 ? size(m.G, j) : size(m.T[i], j)
==(m1::Lar.Model, m2::Lar.Model)   = m1.G == m2.G && m1.T == m2.T
isempty(m::Lar.Model, d::Int)      = isempty(m.T[d])

#-------------------------------------------------------------------------------
#   OTHER BASE REIMPLEMENTATIONS
#-------------------------------------------------------------------------------

Base.copy(m::Lar.Model)     = Lar.Model(m.G, m.T)
Base.deepcopy(m::Lar.Model) = Lar.Model(Base.deepcopy(m.G), Base.deepcopy(m.T))

#-------------------------------------------------------------------------------
#   GEOMETRY MANIPULATION
#-------------------------------------------------------------------------------

function addModelVertices!(m::Lar.Model, V::Lar.Points)::Nothing
    length(m) == size(V, 1) ||
        throw(ArgumentError("Point dimension mismatch."))
    m.G = [m.G V];

    m.T[1] = SparseArrays.sparse(
        findnz(m.T[1])...,
        size(m, 1, 1),
        size(m, 1, 2) + size(V, 2)
    )
    return
end

function addModelVertex!(m::Lar.Model, v::Array{Float64, 1})::Nothing
    return addModelVertices!(m, v[:, :])
end

function deleteModelVertex!(m::Lar.Model, v::Int)::Nothing
    deleteModelVertices!(m, [v])
end

function deleteModelVertices!(m::Lar.Model, vs::Array{Int,1})::Nothing
    tokeepV = setdiff(collect(1 : size(m, 0, 2)), vs);
    m.G     = m.G[:, tokeepV]

    if !isempty(m.T[1])
        # From EV must be deleted all the columns linked to a deleted vertex
        #  and all rows related to a dangling edge
        todelE = abs.(m.T[1][:, vs])
        todelE = [l for l = 1 : todelE.m if sum(todelE[l, :]) != 0]
        m.T[1] = m.T[1][:, tokeepV]
        if !isempty(todelE) Lar.deleteModelCells!(m, 1, todelE) end
    end

    return
end

#-------------------------------------------------------------------------------
#   TOPOLOGY MANIPULATION
#-------------------------------------------------------------------------------

function addModelCells!(m::Lar.Model, deg::Int, cs::Lar.ChainOp)::Nothing
    deg > 0 || throw(ArgumentError("Degree must be a non negative value"))
    deg ≤ length(m) || throw(ArgumentError("The model do not have degree $deg"))
    (cs.n == size(m, deg - 1, 1) && deg > 1) ||
        (cs.n == size(m, deg - 1, 2) && deg == 1) ||
        throw(ArgumentError("Incoherent Chain Operator size"))

    # Adding input cells to the topological structure
    m.T[deg] = [m.T[deg]; cs]

    # Adding slot for the new cells to the higher order cells
    if deg < length(m)
        m.T[deg + 1] = [
            m.T[deg + 1] SparseArrays.spzeros(Int8, size(m, deg + 1, 1), cs.m)
        ]
    end

    return
end
function addModelCell!(m::Lar.Model, deg::Int, c::Lar.Cell)::Nothing
    Lar.addModelCells!(m, deg, convert(Lar.ChainOp, c))
end
function addModelCells!(m::Lar.Model, deg::Int, cs::Lar.Cells)::Nothing
    I = Array{Int,1}()
    J = Array{Int,1}()
    K = Array{Int8,1}()
    for i = 1 : length(cs)
        for j = 1 : length(cs[i])
            push!(I, i)
            push!(J, cs[i][j])
            push!(K, 1)
        end
    end

    scs = SparseArrays.sparse(I, J, K, length(cs), size(m, deg, 2));

    return addModelCells!(m, deg, scs);
end
function addModelCell!(m::Lar.Model, deg::Int, c::Array{Int64,1})::Nothing
    Lar.addModelCells!(m, deg, [c])
end

function deleteModelCells!(m::Lar.Model, deg::Int, cs::Array{Int, 1})::Nothing
    !isempty(cs) || return
    deg > 0 || throw(ArgumentError("Degree must be a non negative value"))
    deg ≤ length(m) || throw(ArgumentError("The model do not have degree $deg"))
    !isempty(m.T[deg]) ||
        throw(ArgumentError("There are no cells of degree $deg"))
    max(cs...) ≤ size(m, deg, 1) ||
        throw(ArgumentError("There are not $(max(cs...)) cells of degree $deg"))

    # Removing `cs` rows from `m.T[deg]`
    tokeep     = setdiff(collect(1 : size(m, deg, 1)), cs)
    m.T[deg]   = m.T[deg][tokeep, :]

    if deg == length(m) return end

    # Removing `cs` cols from `m.T[deg+1]` by checking if some `deg+1` cell
    #  has to be removed too.
    # An higher ord cell must be removed if it contains a removed lower ord cell
    todelHo    = m.T[deg + 1][:, cs]
    m.T[deg+1] = m.T[deg + 1][:, tokeep]
    if isempty(m.T[deg+1]) return end
    todelHo    = [l for l = 1 : todel.m if sum(todelE[l, :]) != 0]
    if !isempty(todelHo) Lar.deleteModelCells!(m, deg+1, todelHo) end
    return
end

function deleteModelCell!(m::Lar.Model, deg::Int, c::Int)::Nothing
    Lar.deleteModelCells!(m, deg, [c])
end

function getModelLoCell(m::Lar.Model, deg::Int, c::Int)::Array{Int, 1}
    deg > 0 || throw(ArgumentError("Degree must be a non negative value"))
    deg ≤ length(m) || throw(ArgumentError("The model do not have degree $deg"))
    return m.T[deg][c, :].nzind
end

function getModelCellVertices(m::Lar.Model, deg::Int, c::Int, ret_idx=false)
    deg > 0 || throw(ArgumentError("Degree must be a non negative value"))
    deg ≤ length(m) || throw(ArgumentError("The model do not have degree $deg"))
    set = [c]
    for d = deg : -1 : 1
        set = ∪([m.T[d][el, :].nzind for el in set]...)
    end

    if ret_idx == true
        return (map(i -> m.G[:, set[i]], 1:length(set)), set)
    end
    return map(i -> m.G[:, set[i]], 1:length(set))
end

function getModelCellGeometry(m::Lar.Model, deg::Int, c::Int, ret_idx=false)
    deg > 0 || throw(ArgumentError("Degree must be a non negative value"))
    deg ≤ length(m) || throw(ArgumentError("The model do not have degree $deg"))
    set = [c]
    for d = deg : -1 : 1
        set = ∪([m.T[d][el, :].nzind for el in set]...)
    end

    if ret_idx == true
        return (m.G[:, set], set)
    end
    return m.G[:, set]
end

#-------------------------------------------------------------------------------
#   MODEL MANIPULATION
#-------------------------------------------------------------------------------

"""
    modelPurge!(m::Lar.Model, [depth::Int = 0])::Nothing

Purges the model Geometry from unused elements.

Purges each Geometry element that is not related to a topological Cell.
If the optional argument `depth` is specified than it also purges all the
topological cells of degree lower or equal to depth that are not related
to an higher order topological cell.
"""
function modelPurge!(m::Lar.Model, depth::Int = 0)::Nothing
    depth < length(m) ||
        throw(ArgumentError("The model do not have degree $(depth+1)"))

    for d = depth : -1 : 1
        todel = [i for i = 1 : size(m, d+1, 2) if sum(abs.(m.T[d+1][:,i])) == 0]
        if !isempty(todel) Lar.deleteModelCells!(m, d, todel) end
    end
    todel = [i for i = 1 : size(m, 1, 2) if sum(abs.(m.T[1][:, i])) == 0]
    if !isempty(todel) Lar.deleteModelVertices!(m, todel) end
    return
end

#-------------------------------------------------------------------------------
#   MODEL CONVERSION TOPOLOGICAL ↔ GEOMETRICAL
#-------------------------------------------------------------------------------

#=
function getGeometricalModel(tm::Lar.Model)::Lar.Model
    tm.topological || throw(ArgumentError("Model must be Topological"))
    gm = deepcopy(tm)
    gm.topological = false
    for i = 2 : length(gm)
        for j = i-1 : -1 : 1
            gm.T[i] = abs.(gm.T[i])*gm.T[j]
        end
    end
    return gm
end

function getTopologicalModel(gm::Lar.Model)::Lar.Model
    !gm.topological || throw(ArgumentError("Model must be Geometrical"))
    tm = deepcopy(gm)
    tm.topological = true
    for i = 2 : length(tm)
        for j = i-1 : -1 : 1
            #...
            throw(Error("Not Coded"))
        end
    end
    return tm
end
=#

#TODO check if removing sparse zeros is necessary

#-------------------------------------------------------------------------------
#   MODEL INTERACTION
#-------------------------------------------------------------------------------

function uniteModels!(model::Lar.Model, m2::Lar.Model)::Nothing

    length(model) == length(m2) ||
        throw(ArgumentError("ERROR: Inconsistent models dimension!"))

    Lar.addModelVertices!(model, m2.G)
    for d = 1 : length(model)
        shift = size(model, d, 2) - size(m2, d, 2)
        I, J, K = SparseArrays.findnz(m2.T[d])
        J = J .+ shift
        Lar.addModelCells!(model, d, SparseArrays.sparse(
            I, J, K, size(m2, d, 1), size(m2, d, 2) + shift
        ))
    end

    return
end

function uniteModels(m1::Lar.Model, m2::Lar.Model)::Lar.Model
    model = Lar.deepcopy(m1)
    Lar.uniteModels!(model, m2)
    return model
end

function mergeMultipleModels(models::Array{Lar.Model,1}; err=1e-6)::Lar.Model

    model = deepcopy(models[1]);

    for i = 2 : length(models)
        Lar.uniteModels!(model, models[i])
    end

    return Lar.mergeModelVertices(model, err=err)
end

function mergeModelVertices(model::Lar.Model; err=1e-4)
    V = convert(Lar.Points, model.G')
    EV = model.T[1]
    FE = model.T[2]
    vertsnum = size(V, 1)
    edgenum = size(EV, 1)
    facenum = size(FE, 1)
    newverts = zeros(Int, vertsnum)
    # KDTree constructor needs an explicit array of Float64
    V = Array{Float64,2}(V)
    W = convert(Lar.Points, LinearAlgebra.transpose(V))
    kdtree = KDTree(W)
	# remove vertices congruent to a single representative
    todelete = []
    i = 1
    for vi in 1:vertsnum
        if !(vi in todelete)
            nearvs = Lar.inrange(kdtree, V[vi, :], err)
            newverts[nearvs] .= i
            nearvs = setdiff(nearvs, vi)
            todelete = union(todelete, nearvs)
            i = i + 1
        end
    end
    nV = V[setdiff(collect(1:vertsnum), todelete), :]

    # translate edges to take congruence into account
    edges = Array{Tuple{Int, Int}, 1}(undef, edgenum)
    oedges = Array{Tuple{Int, Int}, 1}(undef, edgenum)
    for ei in 1:edgenum
        v1, v2 = EV[ei, :].nzind
        edges[ei] = Tuple{Int, Int}(sort([newverts[v1], newverts[v2]]))
        oedges[ei] = Tuple{Int, Int}(sort([v1, v2]))
    end
    nedges = union(edges)
    # remove edges of zero length
    nedges = filter(t->t[1]!=t[2], nedges)
    nedgenum = length(nedges)
    nEV = spzeros(Int8, nedgenum, size(nV, 1))

    etuple2idx = Dict{Tuple{Int, Int}, Int}()
    for ei in 1:nedgenum
    	begin
        	nEV[ei, collect(nedges[ei])] .= 1
        	nEV
        end
        etuple2idx[nedges[ei]] = ei
    end
    for e in 1:nedgenum
    	v1,v2 = findnz(nEV[e,:])[1]
    	nEV[e,v1] = -1; nEV[e,v2] = 1
    end

    # compute new faces to take congruence into account
    faces = [[
        map(x->newverts[x], FE[fi, ei] > 0 ? oedges[ei] : reverse(oedges[ei]))
        for ei in FE[fi, :].nzind
    ] for fi in 1:facenum]

    visited = []
    function filter_fn(face)

        verts = []
        map(e->verts = union(verts, collect(e)), face)
        verts = Set(verts)

        if !(verts in visited)
            push!(visited, verts)
            return true
        end
        return false
    end

    nfaces = filter(filter_fn, faces)

    nfacenum = length(nfaces)
    nFE = spzeros(Int8, nfacenum, size(nEV, 1))

    for fi in 1:nfacenum
        for edge in nfaces[fi]
            ei = etuple2idx[Tuple{Int, Int}(sort(collect(edge)))]
            nFE[fi, ei] = sign(edge[2] - edge[1])
        end
    end

    model = Lar.Model(convert(Lar.Points, nV'))
    Lar.addModelCells!(model, 1, nEV)
    Lar.addModelCells!(model, 2, nFE)
    return model
end

#=
function mergeModelVertices(model::Lar.Model; err=1e-6)
    dims = length(model)
    Gfloat = convert(Array{Float64,2}, model.G)
    kdtree = NearestNeighbors.KDTree(Gfloat)
    todel  = []
    tokeep = []
    I      = [Array{Int, 1}() for i = 1 : dims]
    J      = [Array{Int, 1}() for i = 1 : dims]
    K      = [Array{Int8,1}() for i = 1 : dims]

    # Collide Vertices and update 1-Cochain
    let nidx = 1;  for vidx = 1 : size(model, 0, 2)
        if !(vidx in todel)
            nearvs = NearestNeighbors.inrange(kdtree, Gfloat[:, vidx], err)
            # for d = 1 : dims
            if !isempty(model, 1)
                # @assert sum(sum(abs.(model.T[1][:, nearvs]), dims=2).>1) == 0 "IMPRECISION"
                newVL = sum(model.T[1][:, nearvs], dims=2)
                for i = 1 : length(newVL)  if abs(newVL[i]) == 1
                    push!(I[1], i)
                    push!(J[1], nidx)
                    push!(K[1], newVL[i])
                end  end
            end
            # end
            nidx = nidx + 1
            push!(todel, setdiff(nearvs, vidx)...)
            push!(tokeep, vidx)
        end
    end  end

    # Copy Upper Degree Topology
    for d = 2 : dims
        (I[d], J[d], K[d]) = SparseArrays.findnz(model.T[d])
    end

    model = Lar.Model(model.G[:, tokeep])

    # Collide Topology                                                          #TODO check
    let lolength = length(tokeep);  for d = 1 : dims  if !isempty(I[d])
        # Define compact indices for `d` dimension

        # Recostructiong Topology of degree `d`
        Td = SparseArrays.sparse(I[d], J[d], K[d], max(I[d]...), lolength)
        # Convert Topology to Array of Arrays
        aaT = [SparseArrays.findnz(Td[row, :])[1] for row = 1 : size(Td, 1)]
        nEl = length(aaT)
        # Find new Topology as array of arrays
        newT = sort(unique(aaT))
        # get the indices of the ordered topology
        ordT = sortperm([aaT; newT])
        # building the new Topology as SparseMatrix
        In = Array{Int , 1}()
        Jn = Array{Int , 1}()
        Kn = Array{Int8, 1}()
        let newidx = 1;  for eidx = 1 : length(ordT)  if newidx == 1 || ordT[eidx] < ordT[eidx-1]
            #                            v---co_edge------
            Je, Ke = SparseArrays.findnz(Td[ordT[eidx], :])
            push!(In, (ones(Int, length(Je)) * newidx)...)
            push!(Jn, Je...)
            push!(Kn, Ke...)
            newidx = newidx+1
        end  end  end
        Tn = SparseArrays.sparse(In, Jn, Kn, length(newT), lolength)
        # building the map originalEdge -> (collidedEdge, orientation)
        perm = [(0, 0) for i = 1 : nEl]
        let newidx = 1;  for e in ordT
            e <= nEl ? perm[e] = (newidx, Td[e, :] == Tn[newidx, :] ? 1 : -1) : newidx = newidx + 1
        end  end

        # updating HO topology indices
        for i = 1 : length(J[d+1])
            J[d+1][i], K[d+1][i] = perm[J[d+1][i]] .* (1, K[d+1][i])
        end

        # setting the new topology in the model
        Lar.addModelCells!(model, d, Tn);

        # getting the next dimension topology correct size
        lolength = size(model, d, 1)
    end  end  end

    # Last checks
    for d = 1 : length(model)
        todel = Array{Int, 1}()
        for i = 1 : size(model, d, 1)
            if sum(abs.(model.T[d][i, :])) <= d
                push!(todel, i)
            end
        end
        Lar.deleteModelCells!(model, d, todel)
    end

    return model
end
=#

#-------------------------------------------------------------------------------
#   MODEL REPRESENTATION
#-------------------------------------------------------------------------------

function view2DModel(model::Lar.Model, exp::Float64 = 1.)
    # visualization of numbered arrangement
    #VV = [[k] for k = 1 : size(model, 0, 2)]
    #EV = LAR.cop2lar(model.T[1]);
    #FE = LAR.cop2lar(model.T[2]);
    #FV = cat([[EV[e] for e in face] for face in FE])

    # final solid visualization
    FVs = Lar.triangulate2D(model)
    GL.VIEW(GL.GLExplode(model.G,FVs,exp,exp,exp,99,1));

    # polygonal face boundaries
    EVs = Lar.FV2EVs(model)
    GL.VIEW(GL.GLExplode(model.G,EVs,exp,exp,exp,3,1));
end

function viewModel(model::Lar.Model, depth::Int64 = -1, exp::Float64 = 1.)
    if depth == -1  depth = length(model)  end

    if depth == 0
        GL.VIEW( [GL.GLPoints(convert(Lar.Points, model.G'))] )
    end

    if depth == 1
        if length(model) < 3 || !isempty(model.T[3])
            copCE = abs.(model.T[3]) * abs.(model.T[2]) .!= 0
            comp = [(copCE[c, :]).nzind for c = 1 : size(copCE, 1)]
        elseif !isempty(model.T[2])
            comp = [(model.T[2][c, :]).nzind for c = 1 : size(model.T[2], 1)]
        else
            comp = Lar.Arrangement.biconnected_components(model.T[1])
        end

        mesh = [
            GL.GLGrid(
                model.G,
                LAR.cop2lar(model.T[1][comp[i], :]),
                GL.COLORS[(i-1)%12+1],
                1.
            )
            for i = 1 : length(comp)
        ]
        GL.VIEW( mesh )
    end

    if depth == 2
        model = deepcopy(model)
        modelPurge!(model, 1)
        FVs = Lar.triangulate2D(model)
        GL.VIEW(GL.GLExplode(model.G,FVs,exp,exp,exp,99,1));
    end

    if depth == 3
        GL.VIEW(  )
    end
end

#-------------------------------------------------------------------------------
#   MODEL EXPORT
#-------------------------------------------------------------------------------

function jlexportModel(
        model::Lar.Model,
        filename::String,
        modelname::String = "model"
    )::Nothing

    open(filename, "w") do f
        # Write imports
        write(f, "# using Lar\n")
        write(f, "using SparseArrays\n\n")

        # Write Geometry
        write(f, "#=== GEOMETRY ===#\n")
        write(f, "$modelname = Lar.Model([")
        for d = 1 : length(model)
            write(f, "\n\t")
            for p in model.G[d, :]  write(f, "$p ")  end
        end
        write(f, "\n]);\n\n")

        # Write Topology
        for d = 1 : length(model)
            write(f, "#=== TOPOLOGY : $d-CELLS ===#\n")
            I, J, K = findnz(model.T[d])
            write(f, "I = Array{Int64}([")
            for i in I  write(f, "$i, ")  end
            write(f, "]);\nJ = Array{Int64}([")
            for j in J  write(f, "$j, ")  end
            write(f, "]);\nK = Array{Int8}([")
            for k in K  write(f, "$k, ")  end
            write(f, "]);\n\n")
            write(f, "Lar.addModelCells!($modelname, $d, SparseArrays.sparse(")
            write(f, "I, J, K, $(size(model, d, 1)), $(size(model, d, 2))")
            write(f, "))\n\n\n")
        end
    end

    return
end


text:
using CUDA
using LinearAlgebra
using BenchmarkTools

function ms!(X::CuArray, s::Number)
    function kernel(X, s)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        @inbounds X[i] *= s
        return
    end
    @cuda blocks=length(X) kernel(X, s)
    X
end

function ms1!(X::CuArray, a, b, c)
    k = f(a)
    function kernel(X, a, b, c)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        k(X, i)
        k(X, i)
        k(X, i)
        #@inbounds X[i] *= a
        #@inbounds X[i] *= b
        #@inbounds X[i] *= c
        return
    end
    @cuda blocks=length(X)÷256 threads=256 kernel(X, a, b, c)
    X
end

@inline function f(s)
    @inline function kernel(X, i)
        X[i]*=s
    end
end
function ms2!(X::CuArray, s::Number)
    function kernel(X, s)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        @inbounds X[i] *= s
        return
    end
    @cuda blocks=length(X)÷256 threads=256 kernel(X, s)
    X
end

a = randn(1<<10)
cua = cu(a)
@benchmark ms1!($cua, 1.001, 1.001, 1.001)
@benchmark ms2!(ms2!(ms2!($cua, 1.001), 1.001), 1.001)
bss = 1:length(cua)
@benchmark f(1.001).(cua, bss)

@benchmark rmul!($a, 1.01)
@benchmark ms!($cua, 1.0)


text:
struct ZeroField{N} <: AbstractField{Nothing, Nothing, Nothing, Nothing, Nothing, N} end

ZeroField() = ZeroField{3}() # default

@inline Base.getindex(::ZeroField, ind...) = 0


text:
# Numerically stable implementation of `w -> log(sum(exp, w))`.
# https://discourse.julialang.org/t/fast-logsumexp/22827/9
function logsumexp!(w)
    N = length(w)
    offset, maxind = findmax(w)
    w .= exp.(w .- offset)
    Σ = _sum_all_but(w, maxind)
    log1p(Σ) + offset
end

# Add all elements of vector `w` except for index `i`.
# The element at index `i` is assumed to have value 1
function _sum_all_but(w, i)
    w[i] -= 1
    s = sum(w)
    w[i] += 1
    s
end

"""
    fdot(f, u, v) -> Number

A generic, allocation-free implementation of `dot(u, f.(v))`. It may be faster
to provide a specialized method to dispatch to BLAS or so forth.
"""
function fdot(f, u, v)
    T = promote_type(eltype(u), eltype(v))
    s = zero(T)
    @inbounds for i in eachindex(u, v)
        s += conj(u[i]) * f(v[i])
    end
    s
end

"""
    cost_matrix([C,] a, b) -> Matrix

Precompute the cost matrix given a cost function `C`. If no function `C` is provided, the default is `C(x,y) = norm(x-y)`.
"""
cost_matrix(C, a, b) = [C(x, y) for x in a.set, y in b.set]
cost_matrix(a, b) = cost_matrix((x, y) -> norm(x - y), a, b)


text:
#=
# The BUGS Book
# Example 2.1.2, page 17:
# For
#   Y ~ Binomial(0.5, 8)
# we want to know Pr(Y <= 2)
#

  See ~/jags/bugs_book_2_1_2.jags
       Mean     SD Naive SE Time-series SE
  P2 0.1442 0.3513 0.001434       0.001434
  Y  4.0025 1.4177 0.005788       0.005839

  ~/webppl/bugs_book_2_1_2.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function bugs_book_2_1_2()

    y ~ Binomial(8,0.5)
    p2 ~ y <= 2 ? flip(1.0) : flip(0.0);

end

model = bugs_book_2_1_2()

num_chns = 4
# chns = sample(model, Prior(), MCMCThreads(), 10_000, num_chns)

# chns = sample(model, MH(), MCMCThreads(), 40_000, num_chns)

# chns = sample(model, PG(15), MCMCThreads(), 10_000, num_chns)

# chns = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chns)

chns = sample(model, IS(), MCMCThreads(), 10_000, num_chns)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns,:y)
show_var_dist_pct(chns,:p2)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
run(`gcc -c -Wall -fpic uv.c`)
run(`gcc -shared -o libshimuv.so uv.o`)


text:
"""Start from rest or load initial conditions from netCDF file."""
function initial_conditions(starti=-1)
    # initialise the state matrices u,v,η and set their initial conditions

    if initial_cond == "rest"
        # sub domain sub_dom is the entire domain for a single process (no domain decomposition)
        # u = zeros(Numtype,sub_dom["nux"],sub_dom["nuy"])
        # v = zeros(Numtype,sub_dom["nvx"],sub_dom["nvy"])
        # η = zeros(Numtype,sub_dom["nx"],sub_dom["ny"])
        u = zeros(Numtype,nux,nuy)
        v = zeros(Numtype,nvx,nvy)
        η = zeros(Numtype,nx,ny)
        sst = sst_initial()

    elseif initial_cond == "ncfile"

        #TODO for domain decomposition: slice ncfile for different processors
        inirunpath = initpath*"run"*@sprintf("%04d",init_run_id)*"/"

        # take starti time step from existing netcdf files
        ncu = NetCDF.open(inirunpath*"u.nc")

        if starti == -1
            # replace -1 with length of time dimension
            starti = size(ncu.vars["t"])[end]
        end

        u = ncu.vars["u"][:,:,starti]
        NetCDF.close(ncu)

        ncv = NetCDF.open(inirunpath*"v.nc")
        v = ncv.vars["v"][:,:,starti]
        NetCDF.close(ncv)

        ncη = NetCDF.open(inirunpath*"eta.nc")
        η = ncη.vars["eta"][:,:,starti]
        NetCDF.close(ncη)

        if sstrestart
            ncsst = NetCDF.open(inirunpath*"sst.nc")
            sst = ncsst.vars["sst"][:,:,starti]
            NetCDF.close(ncsst)

            sst = Numtpe.(reshape(sst,size(sst)[1:2]))
        else
            sst = sst_initial()
        end

        # remove singleton time dimension
        # and convert from Float32 to Numtype
        u = Numtype.(reshape(u,size(u)[1:2]))
        v = Numtype.(reshape(v,size(v)[1:2]))
        η = Numtype.(reshape(η,size(η)[1:2]))
    end

    return u,v,η,sst
end

"""Zonally constant hyperbolic tangent initial conditions for the tracer determined by SSTmax, SSTmin, SSTw, SSTϕ."""
function sst_south()
    xx_T,yy_T = meshgrid(x_T,y_T)
    sst = (SSTmin+SSTmax)/2 .+ tanh.(2π*(Ly/(4*SSTw))*(yy_T/Ly .- SSTϕ))*(SSTmin-SSTmax)/2
    return Numtype.(sst)
end

"""Meriodionally constant hyperbolic tangent initial conditions for the tracer determined by SSTmax, SSTmin, SSTw, SSTϕ"""
function sst_west()
    xx_T,yy_T = meshgrid(x_T,y_T)
    sst = (SSTmin+SSTmax)/2 .+ tanh.(2π*(Lx/(4*SSTw))*(xx_T/Lx .- SSTϕ))*(SSTmin-SSTmax)/2
    return Numtype.(sst)
end

"""Initial conditions as a rectangle with SSTmax inside, SSTmin outside."""
function sst_rect()

    x0,x1 = 0.05,0.45     # left right border in [0,1]
    y0,y1 = 0.,1.0     # north south border in [0,1]

    xx_T,yy_T = meshgrid(x_T,y_T)

    sst = fill(SSTmin,size(xx_T))
    inside = (xx_T/Lx .> x0) .* (xx_T/Lx .< x1) .* (yy_T/Ly .> y0) .* (yy_T/Ly .< y1)
    sst[inside] .= SSTmax
    return Numtype.(sst)
end

if injection_region == "south"
    sst_initial = sst_south
elseif injection_region == "west"
    sst_initial = sst_west
elseif injection_region == "rect"
    sst_initial = sst_rect
end

max_stars_repo_name: eschnett/juls

max_issues_repo_name: eschnett/juls

max_forks_repo_name: eschnett/juls


text:
using Test
using FluxRM

if !haskey(ENV, "FLUX_URI")
    flux_available = try
        success(`flux env`)
    catch
        false
    end

    if !flux_available
        @error "flux-core is not installed, skipping further tests"
        exit()
    end

    @info "relaunching under Flux"
    current_file = @__FILE__ # bug in 1.5 can't be directly interpolated
    jlcmd = `$(Base.julia_cmd()) $(current_file)`
    cmd = `flux start -o,-Slog-forward-level=7 -- $jlcmd`
    @test success(pipeline(cmd, stdout=stdout, stderr=stderr))
    exit()
end


text:
using ApproxFun, Base.Test

## Dual Numbers

f=Fun(exp,Interval(DualNumbers.dual(1.0,1),DualNumbers.dual(2.0)),20)
@test_approx_eq Fun(h->Fun(exp,Interval(1.0+h,2.0)).coefficients[1],[0.,.1])'(0.) DualNumbers.epsilon(f.coefficients[1])

## Eig test #336

d = Interval(0.,π)
A=Derivative(d)^2
λ=eigvals([dirichlet(d);A],100)
@test_approx_eq sort(λ)[end-5:end] -(-6:-1).^2

F = x->x.^8
d = Interval(0.0,1.0)
f = Fun(F,d)
ginf = Fun(x->exp(-x),d)
gp = ginf'
Af = Fun(x->x+f(x),d)
transport_ = Fun(x-> x - 1,d)
damping = Fun(x-> 1 - f(x),d)
A = transport_*Derivative(d) + damping
P = -DefiniteIntegral(Chebyshev(d))[LowRankFun((x,y)->gp(x)*(y+f(y)),d^2)];
λ,V = ApproxFun.eigs([A],100)
@test norm(sort(real(filter(x->isreal(x),λ)))[1:5]-(0:4)) ≤ 100000eps()

λ,V = ApproxFun.eigs([A+P],100)
@test_approx_eq_eps sort(real(filter(x->isreal(x),λ)))[5] 3.93759261234502 1E-3


text:
track(m, tape = ReverseDiff.InstructionTape()) = fmap(x -> x isa AbstractArray ? ReverseDiff.track(x, tape) : x, m)
untrack(m) = fmap(ReverseDiff.value, m)

function Flux.Optimise.update!(opt, x, x̄)
    Δ = -Flux.Optimise.apply!(opt, ReverseDiff.value(x), ReverseDiff.value(x̄))
    x_data = ReverseDiff.value(x)
    x_deriv = ReverseDiff.deriv(x)
    x_data .+= ReverseDiff.value(Δ)
    x_deriv .= 0
    return x
end


text:
export pred_ared_computation

"""
This function computes the actual and predicted reduction for one dimensionnal
TR or ARC algorithm. It also returns the ratio bewtween the two.
Inputs:
  - gₖ, the value of the current derivative of the function
  - fₖ, the value of the current value of the function
  - dersec, the value of the current second derivative of the function
  - d, the value of the descent direction
  - ftestTR, the candiate fonction value (h(x + t * d))
  - gtestTR, the candiate derivative value (h'(x + t * d))
  - (opt) seuil, numerical tolerance for the actual reduction. Default value
    is -1e-10.
"""
function pred_ared_computation(gₖ :: Float64, fₖ :: Float64,  dersec :: Float64,
                               d        :: Float64, ftestTR :: Float64,
                               gtestTR  :: Float64; seuil :: Float64 = -1e-10)

  pred = gₖ * d + 0.5 * dersec * d^2
  if pred > seuil
    ared = (gₖ + gtestTR)*d/2
  else
    ared = ftestTR-fₖ
  end

  ratio=ared/pred
  return (pred,ared,ratio)
end


text:
## Matrix helpers
#===========================================================================#

function wake_panels(panels, bound, num)
    lastx, lasty = (p2 ∘ last)(panels)
    bounds = range(lastx, bound, length = num)
    @. WakePanel2D(SVector(bounds[1:end-1], lasty), SVector(bounds[2:end], lasty))
end

ψp(str1, str2, panel_i :: AbstractPanel2D, point) = panel_scalar(vortex_stream_plus, str1, panel_i, first(point), last(point)) + panel_scalar(vortex_stream_minus, str2, panel_i, first(point), last(point))

ψm(str1, str2, panel_i :: AbstractPanel2D, point) = panel_scalar(vortex_stream_plus, str1, panel_i, first(point), last(point)) - panel_scalar(vortex_stream_minus, str2, panel_i, first(point), last(point))

trailing_edge(str1, str2, panel_i :: AbstractPanel2D, point) = panel_scalar(source_stream, str1, panel_i, first(point), last(point)) + panel_scalar(vortex_stream_plus, str2, panel_i, first(point), last(point))

source(str, panel_i, point) = panel_scalar(source_stream, str, panel_i, first(point), last(point))

trailing_edge_source(str, panel_i, point) = panel_scalar(source_stream_trailing_edge, str, panel_i, first(point), last(point))

## Matrix construction
#===========================================================================#

# Need to check singularity conditions
ψp_matrix(panels :: Vector{<: AbstractPanel2D}, points) = [ ψp(1., 1., panel_j, point_i) for (point_i, panel_j) in product(points, panels) ]

ψm_matrix(panels :: Vector{<: AbstractPanel2D}, points) = [ ψm(1., 1., panel_j, point_i) for (point_i, panel_j) in product(points, panels) ]

kutta_condition(points) = [ 1 zeros(length(points) - 2)' 1 ]

streamline(f_pan, l_pan) = (p2(l_pan) - p1(l_pan) + p1(f_pan) - p2(f_pan)) / 2

trailing_edge_influence(s, t, te_panel, pt) = ifelse(p1(te_panel) == pt, 0.0, trailing_edge_source(1., te_panel, pt) * norm(t × s) + ψp(1., 1., te_panel, pt) * (abs ∘ dot)(t, s))

function trailing_edge(panels :: Vector{<: AbstractPanel2D}, pts)
    te_panel  = trailing_edge_panel(panels)
    s         = p2(te_panel) - p1(te_panel)
    t         = streamline(panels[1], panels[end])
    te_vector = trailing_edge_influence.(Ref(s), Ref(t), Ref(te_panel), pts)

    te_panel, te_vector
end

function influence_matrix(panels :: Vector{<: AbstractPanel2D})
    rev_panels  = reverse_panel.(panels)[end:-1:1]
    pts         = panel_points(rev_panels)
    num_pts     = length(pts)

    ψM      = [ ψm_matrix(rev_panels, pts) zeros(num_pts) ]
    ψP      = [ zeros(num_pts) ψp_matrix(rev_panels, pts) ]
    ψM, ψP

    ψ_matrix = ψM .+ ψP

    # Trailing edge setup
    te_panel, te_vector = trailing_edge(rev_panels, pts)
    te_matrix           = [ te_vector zeros(num_pts, num_pts - 2) -te_vector ]
    r_te                = panel_length(te_panel)

    inf_matrix          = ψ_matrix + te_matrix

    ψ0_vector           = fill(-1, length(pts))

    [     inf_matrix        ψ0_vector  ;
      kutta_condition(pts)      0      ]
end

boundary(colpoints, uni) = - [ [ stream(uni, x, y) for (x, y) in colpoints ]; 0 ]

# function aerodynamic_coefficients(panels :: Vector{<: AbstractPanel2D}, γs, u)
#       speed = norm(u)

#       vels  = sum(γs)
#       Δrs   = midpair_map(distance, panels)
#       cps   = @. pressure_coefficient(speed, vels)
#       cls   = @. lift_coefficient(cps, Δrs / 2, panel_angle(panels))

#       cps, cls
# end

# source_matrix(panels_1 :: Vector{<: AbstractPanel2D}, panels_2 :: Vector{<: AbstractPanel2D}) = [ source_influence(panel_j, panel_i) for panel_i ∈ panels_1, panel_j ∈ panels_2 ]

# source_strengths(panels :: Vector{<: AbstractPanel2D}, u) = dot.(Ref(u), panel_normal.(panels))

# boundary_vector(panels :: Vector{<: AbstractPanel2D}, u) = [ - source_matrix(panels, panels) * source_strengths(panels, u); 0 ]

# function boundary_vector(panels :: Vector{<: AbstractPanel2D}, wakes :: Vector{<: AbstractPanel2D}, u)
#       source_panels = [ panels; wakes ]
#       [ - source_matrix(panels, source_panels) * source_strengths(source_panels, u); 0 ]
# end


text:
# This file can be used to create the validation file for deterministic scc values.
# It will create a validation file containing all possibilities of parameter values
# defined in the specs dictionary below produce the same results.

using MimiDICE2010
using DataFrames
using Query
using CSVFiles
using Test

specs = Dict([
    :year => [2015, 2055],
    :eta => [0, 1.5],
    :prtp => [0.015, 0.03],
    :last_year => [2295, 2595],
    :pulse_size => [1e3, 1e7, 1e10]
])

results = DataFrame(year = [], eta = [], prtp = [], last_year = [], pulse_size = [], SC = [])

for year in specs[:year]
    for eta in specs[:eta]
        for prtp in specs[:prtp]
            for last_year in specs[:last_year]
                for pulse_size in specs[:pulse_size]
                    sc = MimiDICE2010.compute_scc(year=Int(year), eta=eta, prtp=prtp, last_year=Int(last_year), pulse_size=pulse_size)
                    push!(results, (year, eta, prtp, last_year, pulse_size, sc))
                end
            end
        end
    end
end

path = joinpath(@__DIR__, "deterministic_sc_values_v1-0-1.csv")
save(path, results)