text:
"""
    TypeCore

Defines basic Types and DataStructures.
Methods binded with Types are separately defined in other modules.
It includes:
1. Stochastic Tools
    1. AbstractMarkovChain <: Any
    2. OrdinaryMarkovChain <: AbstractMarkovChain
2.
"""
module TypeCore
    import LinearAlgebra

    export OrdinaryMarkovChain, AR1
# ==============================

# ------------------
"""
    AbstractMarkovChain

The abstract type of Markov Chains.
"""
abstract type AbstractMarkovChain <: Any end
# ------------------
"""
    OrdinaryMarkovChain(P::Matrix{Float64} ; S::Vector{T} where T <: Real = 1:size(P)[1])

Discrete time, finite number of states,
states in number, constant transition matrix,
Markov Chain.
Receives: `P` - one-step transition matrix `P[i,j]` (prob from state i to state j);
`S` - vector of state values, using 1,2,... in default.

Requires `P`, the transition probability matrix must be square.
An `OrdinaryMarkovChain` must have a stationary distribution.
"""
struct OrdinaryMarkovChain <: AbstractMarkovChain
    N::Int  # the number of states
    S::Vector{T} where T <: Real  # states (in numbers)
    P::Matrix{Float64}  # (one-step) transition matrix P_{ij}, prob from state i to state j
    function OrdinaryMarkovChain(P::Matrix{Float64} ; S::Vector{T} where T <: Real = [x for x in 1:size(P)[1]])
        local Pdim::Tuple = size(P)
        local Pnew::Matrix = zeros(Float64, Pdim)
        @assert(Pdim[1] == Pdim[2], "requires a square matrix")
        @assert(Pdim[1] == length(S), "dimension mismatch between state vector and transition matrix")
        for z in 1:Pdim[1] # rescale each row to make row-sum is exact one
            Pnew[z,:] = LinearAlgebra.normalize(P[z,:], 1) # the 2nd par indicates linearly normalization
        end # for
        return new(Pdim[1], S, Pnew)
    end # constructor
end # OrdinaryMarkovChain

# ---------------------
"""
    AR1

AR(1) process ``x_{t+1}=c+\\rho x_{t}+\\varepsilon_{t}``
where ``\\varepsilon_{t}\\sim N(0,\\sigma)``.
"""
struct AR1 <: Any
    c::Real
    rho::Real
    sigma::Real
end # AR1

# ==============================
end # TypeCore

max_stars_repo_name: Clpr/EconToolBox

max_issues_repo_name: Clpr/EconToolBox

max_forks_repo_name: Clpr/EconToolBox


text:
using CUDA

function jitter!(model::BinomialGridModel, width)
    prob1 = Float32(1/width)
    prob2 = Float32(1/(width-1))

    jitter!(model.Nind, length(model.Nrng), prob1, prob2)
    jitter!(model.pind, length(model.prng), prob1, prob2)
    jitter!(model.qind, length(model.qrng), prob1, prob2)
    jitter!(model.σind, length(model.σrng), prob1, prob2)
    jitter!(model.τind, length(model.τrng), prob1, prob2)

    refresh!(model)
    return model
end

function jitter!(indices, maxindex, prob1, prob2)
    function kernel(indices, range, prob1, prob2, randstates)
        i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
        @inbounds if i <= length(indices)
            r   = GPUArrays.gpu_rand(Float32, CUDA.CuKernelContext(), randstates)
            idx = indices[i]
            if idx == 1
                r < prob2 && (idx += 1)
            elseif idx == maxindex
                r < prob2 && (idx -= 1)
            elseif r < prob1
                idx += 1
            elseif r > 1-prob1
                idx -= 1
            end
            indices[i] = idx
        end#if
        return nothing
    end

    rng = GPUArrays.default_rng(CuArray)

    kernel  = @cuda launch=false kernel(indices, maxindex, prob1, prob2, rng.state)
    config  = launch_configuration(kernel.fun)
    threads = Base.min(length(indices), config.threads, 256)
    blocks  = cld(length(indices), threads)
    kernel(indices, maxindex, prob1, prob2, rng.state; threads=threads, blocks=blocks)
    return indices
end

ext: jl


text:
struct MRT{
    Force,
    T <: Real,
    VT <: AbstractVector{T},
    HST, #<: AbstractVector{AbstractArray{T,1}},
    AST, #<: AbstractVector{AbstractArray{T,1}}
} <: CollisionModel
    τs::VT

    # We will be using Hermite polynomials to compute the corresponding coefficients
    H0::T
    Hs::HST
    # Keep higher order coefficients allocated
    As::AST

    force::Force
end

function MRT(q::Quadrature, τ::T, force = nothing) where {T <: Real}
    N = round(Int, order(q) / 2)
    MRT(q, fill(τ, N))
end

function MRT(q::Quadrature, τ_s::T, τ_a::T, force = nothing) where {T <: Real}
    N = round(Int, order(q) / 2)
    MRT(q, repeat([τ_s, τ_a], outer = N)[1:N])
end

function MRT(q::Quadrature, τs::VT, force = nothing) where {VT <: AbstractVector{<:Real}}
    N = round(Int, order(q) / 2)
    Hs = [[hermite(Val{n}, q.abscissae[:, i], q) for i in 1:length(q.weights)] for n in 1:N]

    MRT(τs, 1.0, Hs, copy(Hs), force)
end

function CollisionModel(cm::Type{<:MRT}, q::Quadrature, problem::FluidFlowProblem)
    τ = q.speed_of_sound_squared * lattice_viscosity(problem) + 0.5
    τs = fill(τ, order(q))

    if has_external_force(problem)
        force_field = (x_idx, y_idx, t) -> lattice_force(problem, x_idx, y_idx, t)

        return MRT(q, τs, force_field)
    end

    return MRT(q, τs)
end

function collide!(collision_model::MRT{Force}, q, f_in, f_out; time = 0.0) where {Force}
    collide_mrt!(collision_model, q, f_in, f_out, time = time)
end

# NOTE: we can do some clever optimizations where we check the value of τ_n
# If it is equal to 1, then we don't need to compute the nth hermite coefficient
# of f
function collide_mrt!(
    collision_model::MRT{Force},
    q,
    f_in::Populations,
    f_out::Populations;
    time = 0.0,
) where {Force, T <: Real, Populations <: AbstractArray{T, 3}}
    cs = q.speed_of_sound_squared

    τs = collision_model.τs

    D = dimension(q)
    N = div(LatticeBoltzmann.order(q), 2)

    # Compute (get!?) the hermite polynomials for this quadrature
    Hs = collision_model.Hs

    nx, ny, nf = size(f_in)

    if !(Force <: Nothing)
        F = zeros(dimension(q))
    end

    f = Array{T}(undef, nf)
    u = zeros(dimension(q))
    @inbounds for x_idx in 1:nx, y_idx in 1:ny
        @inbounds for f_idx in 1:nf
            f[f_idx] = f_in[x_idx, y_idx, f_idx]
        end

        ρ = density(q, f)
        velocity!(q, f, ρ, u)
        # temperature = temperature(q, f, ρ, a_f[1], a_f[2])
        temperature = 1.0

        if !(Force <: Nothing)
            F .= collision_model.force(x_idx, y_idx, time)

            u += τs[2] * F
        end

        # NOTE: we don't need to compute the 0th and 1st coefficient as these are equal
        # to a_f[0] and a_f[1]
        a_eq = [equilibrium_coefficient(Val{n}, q, ρ, u, temperature) for n in 1:N]

        # NOTE: we could optimize this by only computing upto n = 2 when τ = 1.0
        a_f = [sum([f[idx] * Hs[n][idx] for idx in 1:length(q.weights)]) for n in 1:N]

        a_coll = [(1 - 1 / τs[n]) .* a_f[n] .+ (1 / τs[n]) .* a_eq[n] for n in 1:N]

        @inbounds for f_idx in 1:nf
            f_out[x_idx, y_idx, f_idx] =
                q.weights[f_idx] * (
                    ρ +
                    cs * ρ * sum(u .* Hs[1][f_idx]) +
                    sum([
                        cs^n * dot(a_coll[n], Hs[n][f_idx]) / (factorial(n)) for n in 2:N
                    ])
                )
        end
    end
    return
end


text:
"""
Computations related to grid graphs
"""
module Grid
using SparseArrays
import LinearAlgebra
import LinearAlgebra: norm
import ..GraphIdx: enumerate_edges, num_nodes, num_edges, Graph

"""
Position in a grid graph
"""
struct Pixel
    x::Int
    y::Int
end

Base.convert(::Type{Pixel}, p::Tuple{Int,Int}) = Pixel(p[1], p[2])

Base.convert(::Type{Tuple{Int,Int}}, p::Pixel) = (p.x, p.y)

@inline Base.getindex(p::Pixel, i::Int) =
    if i == 1
        p.x
    elseif i == 2
        p.y
    else
        throw(ArgumentErrow("i = $i"))
    end

@inline Base.iterate(p::Pixel, i::Int = 1) =
    if i == 1
        (p.x, 2)
    elseif i == 2
        (p.y, 3)
    else
        nothing
    end

LinearAlgebra.norm2(p::Pixel) = p.x^2 + p.y^2

LinearAlgebra.norm(x::Pixel) = sqrt(LinearAlgebra.norm2(x))

"""
    compute_dirs(dn)

Setting the neighborhood degree (`dn`) to `i` will set all numbers
less or equal `i` in the following diagram as neighbors of `[*]`.
```
[…] […] […] […] […] […] […]
[…] […] [4] [3] [4] […] […]
[…] [4] [2] [1] [2] [4] […]
[…] [3] [1] [*] [1] [3] […]
[…] [4] [2] [1] [2] [4] […]
[…] […] [4] [3] [4] […] […]
[…] […] […] […] […] […] […]
```

# Example

```jldoctest
julia> import GraphIdx.Grid: compute_dirs, Pixel

julia> compute_dirs(2)
2-element Array{Pixel,1}:
 Pixel(1, 0)
 Pixel(1, 1)
```
"""
function compute_dirs(dn::Int)::Vector{Pixel}
    if dn == 1
        Pixel[(1,0)]
    elseif dn == 2
        Pixel[(1,0), (1,1)]
    elseif dn == 3
        Pixel[(1,0), (1,1), (2,0)]
    elseif dn == 4
        Pixel[(1,0), (1,1), (2, 0), (2,1), (1,2)]
    else
        throw(ArgumentError("dn >= 4 not supported"))
    end
end

"""
Capture all we need to know about a grid graph.
"""
struct GridGraph <: Graph
    """Number of rows"""
    n1::Int

    """Number of columns"""
    n2::Int

    """Neighbors (in positive direction) of an inner node"""
    dirs::Vector{Pixel}

    """Include weights in `enumerate_edges`?"""
    weighted::Bool
end

GridGraph(n1, n2, dn::Int = Int(1), weighted::Bool = false) =
    GridGraph(n1, n2, compute_dirs(dn), weighted)

num_edges(g::GridGraph) =
    num_edges(g.n1, g.n2, g.dirs)

num_nodes(g::GridGraph) = g.n1 * g.n2

Base.size(g::GridGraph) = (g.n1, g.n2)

Base.adjoint(g::GridGraph) = GridGraph(g.n2, g.n1, g.dirs, g.weighted)

"""
    iter_edges_pixel(proc, n1, n2, dirs)

Call `proc(i1, j1, i2, j2, len)` for every grid edge
`(i1, j1) -- (i2, j2)` having `len`gth.
"""
@inline function iter_edges_pixel(
    proc::F, n1::Int, n2::Int, dirs::Vector{Pixel}
) where {F<:Function}
    for d::Pixel in dirs
        local len::Float64 = 1/norm(d)
        for i = 1:(n1-d.x)
            for j = 1:(n2-d.y)
                proc(i, j, i+d.x, j+d.y, len)
            end
        end
        for i = 1+d.y:n1
            for j = 1:(n2-d.x)
                proc(i, j, i-d.y, j+d.x, len)
            end
        end
    end
end

@inline function enumerate_edges(
    proc::F, n1::Int, n2::Int, dirs::Vector{Pixel}; weighted::Bool = true
) where {F<:Function}
    @inline pix2ind(i, j) = i + (j-1)*n1
    local no::Int = 0
    if weighted
        for d::Pixel in dirs
            local len::Float64 = 1/norm(d)
            for i = 1:(n1-d.x)
                for j = 1:(n2-d.y)
                    no += 1
                    proc(no, pix2ind(i, j), pix2ind(i+d.x, j+d.y), len)
                end
            end
            for i = 1+d.y:n1
                for j = 1:(n2-d.x)
                    no += 1
                    proc(no, pix2ind(i, j), pix2ind(i-d.y, j+d.x), len)
                end
            end
        end
    else
        for d::Pixel in dirs
            for i = 1:(n1-d.x)
                for j = 1:(n2-d.y)
                    no += 1
                    proc(no, pix2ind(i, j), pix2ind(i+d.x, j+d.y))
                end
            end
            for i = 1+d.y:n1
                for j = 1:(n2-d.x)
                    no += 1
                    proc(no, pix2ind(i, j), pix2ind(i-d.y, j+d.x))
                end
            end
        end
    end
end

@inline enumerate_edges(
    f::F, g::GridGraph; weighted::Bool = false
) where {F<:Function} =
    enumerate_edges(f, g.n1, g.n2, g.dirs, weighted=g.weighted || weighted)

"""
    iter_edges_pixel(proc, g::GridGraph)

Iterate over the edges on a grid graph.
Same as `iter_edges_pixel(::Function, ::Int, ::Int, ::Vector{Pixel}).

# Example

```jldoctest
julia> import GraphIdx.Grid: GridGraph, iter_edges_pixel

julia> iter_edges_pixel(GridGraph(2, 3)) do i1, j1, i2, j2, len
           println("(\$i1,\$j1) -- (\$i2,\$j2): \$len")
       end
(1,1) -- (2,1): 1.0
(1,2) -- (2,2): 1.0
(1,3) -- (2,3): 1.0
(1,1) -- (1,2): 1.0
(1,2) -- (1,3): 1.0
(2,1) -- (2,2): 1.0
(2,2) -- (2,3): 1.0

```
"""
iter_edges_pixel(f::F, g::GridGraph) where {F<:Function} =
    iter_edges_pixel(f, g.n1, g.n2, g.dirs)

"""
    iter_edges(proc, n1, n2, dirs)

Call `proc(u::Int, v::Int, len::Float64)` for every grid
edge `u -- v` having `len`gth.
Hereby `u` and `v` are the index of the corresponding grid-nodes.
"""
function iter_edges(proc::F, n1::Int, n2::Int, dirs::Vector{Pixel}) where {F<:Function}
    """Fortran index of the matrix entry `(i,j)`"""
    pix2ind(i, j) = i + (j-1)*n1

    iter_edges_pixel(n1, n2, dirs) do i1::Int, j1::Int, i2::Int, j2::Int, len::Float64
        proc(pix2ind(i1, j1), pix2ind(i2, j2), len)
    end
end

iter_edges(f::F, g::GridGraph) where {F<:Function} =
    iter_edges(f, g.n1, g.n2, g.dirs)

"""
    num_edges(n1, n2, dirs)
    num_edges(n1, n2, dn::Int)

Number of edges in a grid graph `n1×n2` along directions `dirs`.
If called with last argument a number, first generate `dirs`.
"""
function num_edges(n1::Int, n2::Int, dirs::Vector{Pixel})::Int
    local m::Int = 0
    for d::Pixel in dirs
        m += (n1-d.x)*(n2-d.y) + (n1-d.y)*(n2-d.x)
    end
    return m
end

num_edges(n1::Int, n2::Int, dn::Int = 1)::Int =
    num_edges(n1, n2, compute_dirs(dn))

function incmat(n1::Int, n2::Int, dn::Int = 1)
    local dirs::Vector{Pixel} = compute_dirs(dn)
    local n::Int = n1 * n2
    local m::Int = num_edges(n1, n2, dirs)

    local I::Vector{Int} = Vector{Int}(undef, 2m)
    local J::Vector{Int} = Vector{Int}(undef, 2m)
    local W::Vector{Float64} = zeros(Float64, 2m)

    enumerate_edges(n1, n2, dirs, weighted=true) do k::Int, u::Int, v::Int, l::Float64
        I[2k-1] = k
        J[2k-1] = u
        W[2k-1] = +l
        I[2k-0] = k
        J[2k-0] = v
        W[2k-0] = -l
    end
    D = sparse(I, J, W, m, n)
end

function adjlist(n1::Int, n2::Int, dirs::Vector{Pixel})
    local m::Int = num_edges(n1, n2, dirs)
    local head::Vector{Int} = Vector{Int}(undef, m)
    local tail::Vector{Int} = Vector{Int}(undef, m)
    local lam::Vector{Float64} =  Vector{Float64}(undef, m)

    enumerate_edges(n1, n2, dirs, weighted=true) do k::Int, u::Int, v::Int, l::Float64
        head[k] = u
        tail[k] = v
        lam[k] = l
    end
    head, tail, lam
end

adjlist(n1::Int, n2::Int, dn::Int = 1) =
    adjlist(n1, n2, compute_dirs(dn))

adjlist(g::GridGraph) =
    ajdlist(g.n1, g.n2, g.dirs)

"""
    line_D(n)

Return the incidence matrix for a line graph of length `n`.
Is equivalent to `incmat(1, n, 1)`.

# Example

```jldoctest
julia> import GraphIdx.Grid: line_D

julia> line_D(3)
2×3 SparseArrays.SparseMatrixCSC{Float64,Int64} with 4 stored entries:
  [1, 1]  =  1.0
  [1, 2]  =  -1.0
  [2, 2]  =  1.0
  [2, 3]  =  -1.0

```
"""
function line_D(n::Integer)
    D = spzeros(n-1, n)
    for i = 1:n-1
        D[i, i+0] = +1
        D[i, i+1] = -1
    end
    return D
end

"""
    lipschitz(n1, n2, [dn | dirs])

Compute an upper bound for the Lipschitz-constant for ...
TODO: Be more precise
"""
function lipschitz(n1, n2, dirs::Vector{Pixel})
    lmax = maximum(lam)
    Lip = lmax * 2 * 4 * sum([l for (d,l) in dirs])
    error("not implemented, yet")
end

lipschitz(n1::Int, n2::Int, dn::Int)::Float64 =
    lipschitz(n1, n2, compute_dirs(dn))

function neighboridx_lambda(g::Grid.GridGraph)
    local n::Int = num_nodes(g)
    local m::Int = num_edges(g)
    local idx::Vector{Int} = zeros(Int, n+1)
    iter_edges(g) do u::Int, v::Int, lam_i::Float64
        idx[u] += 1
        idx[v] += 1
    end
    acc = 1                        # accumulate degrgees ==> positions
    deg_i = 0
    deg_ii = idx[1]
    for i = 1:n
        idx[i] = acc
        acc += deg_i
        deg_i, deg_ii = deg_ii, idx[i+1]
    end
    idx[n+1] = acc
    @assert(idx[end] + deg_i == 2m + 1,
            "idx[$(length(idx))]: $(idx[end] + deg_i) != $(2m + 1)")
    local pi::Vector{Tuple{Int,Int}} = Vector{Tuple{Int,Int}}(undef, 2m)
    local lam::Vector{Float64} = Vector{Float64}(undef, m)
    enumerate_edges(g, weighted=true) do i::Int, u::Int, v::Int, lam_i::Float64
        lam[i] = lam_i
        pi[idx[u+1]] = (v, i)
        idx[u+1] += 1
        pi[idx[v+1]] = (u, i)
        idx[v+1] += 1
    end
    @assert(idx[end] == 2m + 1, "$(idx[end]) vs $(2m + 1)")
    return idx, pi, lam
end

"""How is a node within the grid?"""
@enum BorderType begin
    Middle      = 0b000
    North       = 0b100
    NorthWest   = 0b101
    NorthEast
    SouthWest
    SouthEast
    South
    West
    East
end

@inline coordinates(g::Grid.GridGraph, v::Int)::Tuple{Int,Int} =
    coordinates(g.n1, v)

@inline function coordinates(n1::Int, v::Int)::Tuple{Int,Int}
    i, j = divrem(v-1, n1)
    i += 1
    j += 1
    (i, j)
end

function BorderType(g::Grid.GridGraph, v::Int)::BorderType
    i, j = coordinates(g, v)
    if (i, j) == (1, 1)
        NorthWest
    elseif (i, j) == (1, g.n2)
        NorthEast
    elseif (i, j) == (g.n1, 1)
        SouthWest
    elseif (i, j) == (g.n1, g.n2)
        SouthEast
    elseif i == 1
        North
    elseif i == g.n1
        South
    elseif j == 1
        West
    elseif j == g.n2
        East
    else
        Middle
    end
end

"""
Idea: Do not store the graph explicitly but compute the neighbors as needed.
For avoiding too many allocations, return just a view to the buffer
"""
struct ImplicitGridGraph <: Graph
    grid::GridGraph
    _buffer::Vector{Tuple{Int,Int}}

    ImplicitGridGraph(n1::Int, n2::Int, dn::Int = Int(1)) =
        new(GridGraph(n1, n2, dn), [])
end

struct GridNeighbors
    p::Complex{Int}
    bt::BorderType
    graph::GridGraph
end

Base.Complex(p::Pixel)::Complex{Int} =
    Complex{Int}(p.x, p.y)

Pixel(p::Complex{Int}) =
    Pixel(real(p), imag(p))

num_edges(bt::BorderType) =
    if bt == Middle
        4
    elseif bt == South || bt == North || bt == West || bt == East
        2
    else
        1
    end

Base.length(gn::GridNeighbors) =
    num_edges(gn.bt) * length(gn.graph.dirs)

function Base.iterate(gn::GridNeighbors, state::Tuple{Int,Int} = (1, 1))
    pix2ind(p::Pixel) = pix2ind(p...)
    pix2ind(p::Complex{Int}) = pix2ind(Pixel(p))
    pix2ind(px::Int, py::Int) = px + (py-1) * gn.graph.n1
    let (di, ii) = state
        if di > length(gn.graph.dirs)
            return nothing
        end
        local dir::Complex{Int} = gn.graph.dirs[di]
        if gn.bt == Middle
            if ii > 4
                return iterate(gn, (di+1, 1))
            else
                return (pix2ind(gn.p + (1im^ii) .* Complex(dir)), (di, ii+1))
            end
        elseif gn.bt == North
            if ii > 3
                return iterate(gn, (di+1, 1))
            else
                return (pix2ind(gn.p + (1im^ii) .* Complex(dir)), (di, ii+1))
            end
            notimplemented("todo")
        end
    end
end

function Base.getindex(graph::GridGraph, v::Int = 2)
    bt = BorderType(graph, v)
    return GridNeighbors(Complex(coordinates(graph, v)...), bt, graph)
end

"""
    collect_edges(::GridGraph)

Return a `Vector{Tuple{Int,Int}}` of edges and `Vector{Float64}` of weights.
"""
function collect_edges(graph::GridGraph)
    local m = num_edges(graph)
    local n = num_nodes(graph)
    edges = Vector{Tuple{Int,Int}}(undef, m)
    lam = Vector{Float64}(undef, m)
    enumerate_edges(graph, weighted=true) do ei::Int, u::Int, v::Int, l::Float64
        edges[ei] = (u, v)
        lam[ei] = l
    end
    return edges, lam
end

end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
#--------------------------------------------------------------------
# Spacetime Discretization methods in Julia
# Soham 01-2018
#--------------------------------------------------------------------

module ScalarWave
using Luxor, Colors

export Patch, Boundary

export cheb, chebx, chebd, chebw, chebgrid, chebweights, vandermonde, pseudovandermonde,
	   delta, coordtransL2G, coordtransG2L, jacobian, shapeH2L, shapeL2H,
       LInfnorm, L1norm, L2norm, array2dict, dict2array, savegrid, loadgrid,
       derivOP, boundaryOP, RHS, modal2nodal, nodal2modal,
       getPatchIC, getPatchBnd, calcPatch, extractPatchCoeffs, interpolatePatch,
       restrictmodes!, prolongatemodes, restrictOP, prolongateOP, restrictPatch, prolongatePatch,
       projectonPatchBndbyRestriction, projectonPatchbyRestriction,
       distribute, sconv, showconv, drawmultipatch,
       fdistribute, fgetPatchBnd, fgetPatchIC, fRHS, fcalcPatch

include("Types.jl")
include("SpecCalculus.jl")
include("Utilities.jl")
include("Physics.jl")
include("Patch.jl")
include("Grid.jl")
include("Projection.jl")
include("Dispatch.jl")
include("Convergence.jl")
include("Visualization.jl")
include("Futures.jl")

end


text:
using SpreadsheetFunctions
using Test

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


text:
# calculateposterior.jl

"""
    calculatedependentfaberschauderfunctions(higestlevel::Int64)

Internal function, not exported.

Calculates which Faber-Schauder functions (ψ_i,ψ_j) have essential overlapping
support, i ≤ j. All other combinations (ψ_i, ψ_j), i ≤ j, have essential
nonoverlapping support and hence int_0^T ψ_i(X_t)ψ_j(X_t)dt=0, which we use in
calculategirsanovmatrix.

Returns a triple (lengthvectors, rowindices, columnindices) where
```julia
lengthvectors == length(rowindices) == length(columnindices),
```
and for i in 1:lengthvectors (ψ_rowindices[i],ψ_columnindices[i]) have
essentially overlapping support. For all i, rowindices[i] ≤ columnindices[i].
"""
function calculatedependentfaberschauderfunctions(higestlevel::Int64)
    lengthvectors = (higestlevel+1)*2^(higestlevel+1)+1 # Correct.
    rowindices = Vector{Int}(undef, lengthvectors)
    columnindices = Vector{Int}(undef, lengthvectors)
    # Ψ_1 is dependent with all basis functions, including itself.
    rowindices[1:2^(higestlevel+1)] .= 1
    columnindices[1:2^(higestlevel+1)] = 1:2^(higestlevel+1)
    index = 2^(higestlevel+1) + 1 # point where we are in constructing rowindices and columnindices
    # psi_{j,k}=psi_{2^j+k} is dependent with
    # psi_{j+d,(k-1)2^d+1},...,psi_{j+d, k2^d}, d\ge0 (which includes itself).
    for jone in 0:higestlevel
        twotothepowerjone = 2^jone
        for kone in 1:2^jone
            ione = twotothepowerjone+kone # other index system
            for  jtwo in jone:higestlevel
                d = jtwo - jone
                twotothepowerjtwo = 2^jtwo
                for ktwo in (kone-1)*2^d+1:kone*2^d
                    rowindices[index] = ione
                    columnindices[index] = twotothepowerjtwo+ktwo
                    index += 1
                end
            end
        end
    end
    @assert lengthvectors + 1 == index
    return (lengthvectors, rowindices, columnindices)
end

# calculategirsanovmatrixelement methods:

"""
    calculategirsanovmatrixelement(ψ1Xt, ψ2Xt, σ, Δt)
    calculategirsanovmatrixelement(samplevalueindices1, samplevalueindices2, ψ1Xt, ψ2Xt, σ, Δt)

Internal function, not exported!

Where ψ1Xt, ψ2Xt are arrays, and σ and Δt are either arrays or numbers.

Let ψ_a and ψ_b denote two basis elements.
calculategirsanovmatrixelement calculates the (a,b) element of the Girsanov
matrix defined by int_0^T ψ_a(X_t)ψ_b(X_t)/(σ^2(X_t)) dt, where ψ_a is given
by ψ1Xt (already evaluated in ψ_a) and ψ_b by ψ2Xt (already evaluated in ψ_b).

samplevalueindices1[i] == false should correspond to Ψ1Xt[i] == 0.0.
Similar for samplevalueindices2 and Ψ2Xt.
"""
calculategirsanovmatrixelement(ψ1Xt, ψ2Xt, σXt, Δt) = sum(Δt .* ψ1Xt .* ψ2Xt ./ (σXt.^2))
#test written
calculategirsanovmatrixelement(ψ1Xt, ψ2Xt, σXt, Δt::Number) = Δt * sum(ψ1Xt .* ψ2Xt ./ (σXt.^2))
#test written
calculategirsanovmatrixelement(ψ1Xt, ψ2Xt, σXt::Number, Δt) = sum(Δt .* ψ1Xt .* ψ2Xt) / (σXt^2)
#test written
calculategirsanovmatrixelement(ψ1Xt, ψ2Xt, σXt::Number, Δt::Number) = Δt * sum(ψ1Xt .* ψ2Xt) / (σXt^2)
#test written
function calculategirsanovmatrixelement(samplevalueindices1, samplevalueindices2, ψ1Xt, ψ2Xt, σXt, Δt)
    samplevalueindices = samplevalueindices1 .& samplevalueindices2
    filteredψ1Xt = ψ1Xt[samplevalueindices]
    filteredψ2Xt = ψ2Xt[samplevalueindices]
    filteredσXt = σXt[samplevalueindices]
    filteredΔt = Δt[samplevalueindices]
    return sum(filteredΔt .* filteredψ1Xt .* filteredψ2Xt ./ (filteredσXt.^2))
end
function calculategirsanovmatrixelement(samplevalueindices1, samplevalueindices2, ψ1Xt, ψ2Xt, σXt, Δt::Number)
    samplevalueindices = samplevalueindices1 .& samplevalueindices2
    filteredψ1Xt = ψ1Xt[samplevalueindices]
    filteredψ2Xt = ψ2Xt[samplevalueindices]
    filteredσXt = σXt[samplevalueindices]
    return Δt * sum(filteredψ1Xt .* filteredψ2Xt ./
        (filteredσXt.^2))
end
function calculategirsanovmatrixelement(samplevalueindices1, samplevalueindices2, ψ1Xt, ψ2Xt, σXt::Number, Δt)
    samplevalueindices = samplevalueindices1 .& samplevalueindices2
    filteredΔt = Δt[samplevalueindices]
    filteredψ1Xt = ψ1Xt[samplevalueindices]
    filteredψ2Xt = ψ2Xt[samplevalueindices]
    return sum(filteredΔt .* filteredψ1Xt .* filteredψ2Xt) / (σXt^2)
end
function calculategirsanovmatrixelement(samplevalueindices1, samplevalueindices2, ψ1Xt, ψ2Xt, σXt::Number, Δt::Number)
    samplevalueindices = samplevalueindices1 .& samplevalueindices2
    filteredψ1Xt = ψ1Xt[samplevalueindices]
    filteredψ2Xt = ψ2Xt[samplevalueindices]
    return Δt * sum(filteredψ1Xt .* filteredψ2Xt) / (σXt^2)
end

"""
    calculateΔt(timeinterval::AbstractRange) = step(timeinterval)
    calculateΔt(timeinterval) = timeinterval[2:end] - timeinterval[1:end-1]

Internal function, not exported!

Calculates Δt. It returns a number when timeinterval is a range object
and an array of increments of the timeinterval vector otherwise.
"""
calculateΔt(timeinterval::AbstractRange) = step(timeinterval)
calculateΔt(timeinterval) = timeinterval[2:end] - timeinterval[1:end-1]
# Test written.

"""
    calculategirsanovmatrix(sizesquarematrix, timeinterval, ψXt, σXt)
    calculategirsanovmatrix(Π, samplevalueindices, timeinterval, ψXt, σXt)

Internal function, not exported!

Calculates the Girsanov matrix. In the case of a FaberSchauderExpansionWithGaussianCoefficients-object
it will make use of the sparsity structure of the Faber-Schauder basis.
"""
function calculategirsanovmatrix(sizesquarematrix, timeinterval, ψXt, σXt)
    girsanovmatrix = Matrix{Float64}(undef, sizesquarematrix, sizesquarematrix)
    Δt = calculateΔt(timeinterval)
    for k in 1:sizesquarematrix
        for l in k:sizesquarematrix
            girsanovmatrix[k,l] = calculategirsanovmatrixelement(
                ψXt[k], ψXt[l], σXt, Δt)
        end
    end
    return Symmetric(girsanovmatrix)
end
function calculategirsanovmatrix(Π, samplevalueindices, timeinterval, ψXt, σXt)
    lengthvectors, rowindices, columnindices = calculatedependentfaberschauderfunctions(Π.higestlevel)
    d = length(Π)
    V = Vector{Float64}(undef, lengthvectors)
    Δt = calculateΔt(timeinterval)
    for i in 1:lengthvectors
        V[i] = calculategirsanovmatrixelement(samplevalueindices[rowindices[i]],
            samplevalueindices[columnindices[i]], ψXt[rowindices[i]], ψXt[columnindices[i]], σXt, Δt)
    end
    return Symmetric(sparse(rowindices, columnindices, V, d, d))
end

"""
    calculategirsanovvectorelement(ΔXt, ψXt, σXt)
    calculategirsanovvectorelement(ΔXt, ψXt, σXt::Number)

    calculategirsanovvectorelement(samplevalueindices, ΔXt, ψXt, σ::Number)

Internal function, used by calculategirsanovvector, not exported!

Calculates the kth Girsanov vector element int_0^T ψ_k(X_t)/(σ^2(X_t)) dX_t
with the Euler-Maruyama method where ψ_k is the kth basis function, and T is the end time.
"""
calculategirsanovvectorelement(ΔXt, ψXt, σXt) = sum(ΔXt .* ψXt ./ (σXt.^2))
#Test written.
calculategirsanovvectorelement(ΔXt, ψXt, σXt::Number) = sum(ΔXt .* ψXt) / (σXt^2)
#Test written.
function calculategirsanovvectorelement(samplevalueindices, ΔXt, ψXt, σXt)
    filteredΔXt = ΔXt[samplevalueindices]
    filterψXt = ψXt[samplevalueindices]
    filteredσXt = σXt[samplevalueindices]
    return sum(filteredΔXt .* filterψXt ./ (filteredσXt.^2))
end
#Test written.
function calculategirsanovvectorelement(samplevalueindices, ΔXt, ψXt, σXt::Number)
    filteredΔXt = ΔXt[samplevalueindices]
    filteredψXt = ψXt[samplevalueindices]
    return sum(filteredΔXt .* filteredψXt) / (σXt^2)
end
#Test written.
"""
    calculategirsanovvector(lengthvector, samplevalues, ψXt, σXt)
    calculategirsanovvector(lengthvector, samplevalueindices, samplevalues, ψXt, σXt)

Internal function, not exported!

Calculates the Girsanov vector.
"""
function calculategirsanovvector(lengthvector, samplevalues, ψXt, σXt)
    ΔXt = samplevalues[2:end] - samplevalues[1:end-1]
    return [calculategirsanovvectorelement(ΔXt, ψXt[k], σXt) for k in
        1:lengthvector]
end
# Test written.
function calculategirsanovvector(lengthvector, samplevalueindices, samplevalues, ψXt, σXt)
    ΔXt = samplevalues[2:end] - samplevalues[1:end-1]
    return [calculategirsanovvectorelement(samplevalueindices[k], ΔXt, ψXt[k],
        σXt) for k in 1:lengthvector]
end
# Test written.

"""
    calculateσXt(σ::Number, v) = σ
    calculateσXt(σ, v) = σ.(v)

Internal function, not exported!.

Returns σ when it is a number, and if σ is a function, it evaluates v in σ.
"""
calculateσXt(σ::Number, v) = σ
calculateσXt(σ, v) = σ.(v)
#Tests written.

"""
    calculateposterior(Π, X, model::SDEModel)
    calculateposterior(Π, X, σ)
    calculateposterior(Π::FaberSchauderExpansionWithGaussianCoefficients, X, σ)

Calculates the posterior distribution Π(⋅∣X) and returns a FaberSchauderExpansionWithGaussianCoefficients
object when Π is a FaberSchauderExpansionWithGaussianCoefficients-object. Otherwise, it returns a GaussianProcess-object.
Uses model to determine σ. When Π is a FaberSchauderExpansionWithGaussianCoefficients-object, then calculateposterior
makes use of the sparsity structure of the Faber-Schauder basis.

# Examples

##Example with Faber-Schauder expansion.
```julia
model = SDEModel(1.0, 0.0, 10000.0, 0.1)
sde = SDE(x->sinpi(2*x), model)
X = rand(sde)
Π = FaberSchauderExpansionWithGaussianCoefficients([2.0^j for j in 0:5])
postΠ = calculateposterior(Π, X, model)
```
##Example with Fourier expansion.
```julia
using Distributions
α = 0.5
model = SDEModel(1.0, 0.0, 10000.0, 0.1)
sde = SDE(x->sinpi(2*x), model)
X = rand(sde)
distribution = MvNormalCanon([k^(α+0.5) for k in 1:50])
Π = GaussianProcess([fourier(k) for k in 1:50], distribution)
postΠ = calculateposterior(Π, X, model)
```
"""
calculateposterior(Π, X, model::SDEModel) = calculateposterior(Π, X, model.σ)
calculateposterior(Π::FaberSchauderExpansionWithGaussianCoefficients, X, model::SDEModel) = calculateposterior(Π, X, model.σ)
function calculateposterior(Π, X, σ)
    σXt = calculateσXt(σ, X.samplevalues[1:end-1])
    ψXt = [f.(X.samplevalues[1:end-1]) for f in Π.basis]
    lengthΠ = length(Π)
    precisionprior = invcov(Π.distribution)
    girsanovvector = calculategirsanovvector(lengthΠ, X.samplevalues, ψXt, σXt)
    girsanovmatrix = calculategirsanovmatrix(lengthΠ, X.timeinterval, ψXt, σXt)
    precisionmatrixposterior = girsanovmatrix + precisionprior
    potentialposterior = girsanovvector + precisionprior * mean(Π.distribution)
    posteriordistributiononcoefficients = MvNormalCanon(potentialposterior, precisionmatrixposterior)
    return GaussianProcess(Π.basis, posteriordistributiononcoefficients)
end
function calculateposterior(Π::FaberSchauderExpansionWithGaussianCoefficients, X, σ)
    samplevaluesmod1 = mod.(X.samplevalues[1:end-1], 1.0)
    σXt = calculateσXt(σ, X.samplevalues[1:end-1])
    ψXt = [ψ.(X.samplevalues[1:end-1]) for ψ in Π.basis]
    # You only need to use the values which are in the support of the function.
    samplevalueindices = [Π.leftboundssupport[i] .≤
        samplevaluesmod1 .≤ Π.rightboundssupport[i] for i in 1:length(Π)]
    girsanovvector = calculategirsanovvector(length(Π), samplevalueindices,
        X.samplevalues, ψXt, σXt)
    girsanovmatrix = calculategirsanovmatrix(Π, samplevalueindices,
        X.timeinterval, ψXt, σXt)
    precisionprior = invcov(Π.distribution)
    precisionmatrixposterior = girsanovmatrix + precisionprior
    potentialposterior = girsanovvector + precisionprior * mean(Π.distribution)
    posteriordistributiononcoefficients = MvNormalCanon(potentialposterior, precisionmatrixposterior)
    return FaberSchauderExpansionWithGaussianCoefficients(Π.higestlevel, posteriordistributiononcoefficients)
end

function calculateposteriorcoefficients(d::GaussianVector,
        X::SamplePath, model::SDEModel)::GaussianProcess
    σXt = calculateσXt(model.σ, X.samplevalues[1:end-1])
    ψXt = [f.(X.samplevalues[1:end-1]) for f in Π.basis]
    lengthΠ = length(d)
    covariancematrixprior = Matrix(cov(Π.distribution))
    precisionprior = inv(covariancematrixprior)
    girsanovvector = calculategirsanovvector(lengthΠ, X.samplevalues, ψXt, σXt)
    girsanovmatrix = calculategirsanovmatrix(lengthΠ, X.timeinterval, ψXt, σXt)
    precisionmatrixposterior = girsanovmatrix + precisionprior
    potentialposterior = girsanovvector + precisionprior * mean(Π.distribution)
    meanposterior =  covariancematrixprior * potentialposterior
    posteriordistribution = GaussianVector(meanposterior,
        precisionmatrixposterior^(-0.5))
    return GaussianProcess(Π.basis, posteriordistribution)
end

function GibbsIter(siteration, A, B, meanGPPost, precGPPost, numberofbasisfunctions, X, ψXt, σXt)
    meanGPPost /= sqrt(siteration)
    precGPPost *= sqrt(siteration)
    girsanovvector = calculategirsanovvector(numberofbasisfunctions, X.samplevalues, ψXt, σXt)
    girsanovmatrix = calculategirsanovmatrix(numberofbasisfunctions, X.timeinterval, ψXt, σXt)
    precGPPost += Matrix(girsanovmatrix)
    # Lijn hieronder moet de rekenfout hersteld worden.
    potentialGPPost = girsanovmatrix + meanGPPost
    meanGPPost = precGPPost \ potentialGPPost
    distGPPost = GaussianVector(meanGPPost, sqrt(inv(precGPPost)))
    ϴ = rand(distGPPost)
    siteration = rand(Gamma(A+numberofbasisfunctions/2, B + 0.5 * (ϴ - meanGPPost)' *
        potentialGPPost * (ϴ - meanGPPost)))
    return (siteration, meanGPPost, precGPPost, ϴ)
end

function GibbsSampler(Π::S, A::Float64, B::Float64, X::SamplePath, model::SDEModel,
        numSamples::Int, S0 = 1.0, burnin::Int = 100, intercept::Int = 1) where {S<:AbstractGaussianProcess}
    meanGPPost = mean(Π.distribution)
    precGPPost = inv(Matrix(cov(Π.distribution)))
    numberofbasisfunctions = convert(Float64, length(Π))

    σXt = calculateσXt(model.σ, X.samplevalues[1:end-1])
    ψXt = [f.(X.samplevalues[1:end-1]) for f in Π.basis]

    siteration = S0

    for i in 1:burnin
        siteration, meanGPPost, precGPPost = GibbsIter(siteration, A, B, meanGPPost,
            precGPPost, numberofbasisfunctions, X, ψXt, σXt)
    end

    samplePostS = Array{Float64}(undef, numSamples)
    samplePostϴ = Array{Float64}(undef, numSamples, numberofbasisfunctions)
    for i in 1:numSamples
        for j in 1:(intercept+1)
            siteration, meanGPPost, precGPPost, ϴ = GibbsIter(siteration, A, B,
                meanGPPost, precGPPost, numberofbasisfunctions, X, ψXt, σXt)
        end
        samplePostS[i], samplePostϴ[i,:] = siteration, ϴ
    end
    return (samplePostS, samplePostϴ)
end


text:
mutable struct expected_holding_cost{T} <: AbstractCost
    h::T
end

function (f::expected_holding_cost)(y, demand_dist::Union{CVNormal,Normal})
    μ = mean(demand_dist)
    σ = std(demand_dist)
    σ == 0 && return linear_holding_cost(f.h)(y-μ)
    h = f.h
    v = σ^2
    e = MathConstants.e
    π = MathConstants.pi
    fh = zero(y)
    if y >= 0
        fh = (((e^(-((y - μ)^2)/2v) -e^(-(μ^2)/2v))*sqrt(v))/(sqrt(2π))) + 1//2*(y-μ)*(erf((y-μ)/(sqrt(2)*sqrt(v))) + erf(μ/(sqrt(2)*sqrt(v))))
        fh *= h
    end
    return fh
end

mutable struct expected_stockout_cost{T} <: AbstractCost
    b::T
end

function (f::expected_stockout_cost)(y, demand_dist::Union{CVNormal,Normal})
    μ = mean(demand_dist)
    σ = std(demand_dist)
    σ == 0 && return linear_stockout_cost(f.h)(y-μ)
    b = f.b
    v = σ^2
    e = MathConstants.e
    π = MathConstants.pi
    fb = zero(y)
    if y >= 0
        fb = ((e^(-((y - μ)^2)/2v))*v + sqrt(π/2) * (μ - y) * sqrt(v) * erfc((y-μ)/(sqrt(2)*sqrt(v))))/(sqrt(2π)*sqrt(v))
        fb *= b
    else
        fb = (e^(-(μ^2)/2v)*v + sqrt(π/2) * (μ - y) * sqrt(v) * (1 + erf(μ/(sqrt(2)*sqrt(v)))))/(sqrt(2π)*sqrt(v))
        fb *= b
    end
    return fb
end


text:
@testset "Parent Probabilities" begin
    bg = 0.5
    kappa = 0.5
    kern(x) = pdf.(Distributions.Exponential(1/0.5), x)
    @testset "Individual" begin
        inferTests && @inferred HawkesProcesses.event_parent_prob(3, bg, kappa, kern, [1, 2])
        testProbs = HawkesProcesses.event_parent_prob(3, bg, kappa, kern, [1, 2])
        @test length(testProbs) == 3
        @test testProbs[1] == bg
    end

    @testset "Individual - Functional Background" begin
        inferTests && @inferred HawkesProcesses.event_parent_prob(3, x->0.5, kappa, kern, [1, 2])
        testProbs = HawkesProcesses.event_parent_prob(3, x->0.5, kappa, kern, [1, 2])
        @test length(testProbs) == 3
        @test testProbs[1] == 0.5
    end

    @testset "Full" begin

        inferTests && @inferred HawkesProcesses.parent_probs(sort(rand(100)), bg, kappa, kern)
        pp = HawkesProcesses.parent_probs(sort(rand(100)), bg, kappa, kern)

        @test length(pp) == 100

    end

end


text:
using ArrayInterface, Test
using Base: setindex
import ArrayInterface: has_sparsestruct, findstructralnz, fast_scalar_indexing
@test ArrayInterface.ismutable(rand(3))

using StaticArrays
@test ArrayInterface.ismutable(@SVector [1,2,3]) == false
@test ArrayInterface.ismutable(@MVector [1,2,3]) == true

using LinearAlgebra, SparseArrays

D=Diagonal([1,2,3,4])
@test has_sparsestruct(D)
rowind,colind=findstructralnz(D)
@test [D[rowind[i],colind[i]] for i in 1:length(rowind)]==[1,2,3,4]
@test length(rowind)==4
@test length(rowind)==length(colind)

Bu = Bidiagonal([1,2,3,4], [7,8,9], :U)
@test has_sparsestruct(Bu)
rowind,colind=findstructralnz(Bu)
@test [Bu[rowind[i],colind[i]] for i in 1:length(rowind)]==[1,7,2,8,3,9,4]
Bl = Bidiagonal([1,2,3,4], [7,8,9], :L)
@test has_sparsestruct(Bl)
rowind,colind=findstructralnz(Bl)
@test [Bl[rowind[i],colind[i]] for i in 1:length(rowind)]==[1,7,2,8,3,9,4]

Tri=Tridiagonal([1,2,3],[1,2,3,4],[4,5,6])
@test has_sparsestruct(Tri)
rowind,colind=findstructralnz(Tri)
@test [Tri[rowind[i],colind[i]] for i in 1:length(rowind)]==[1,2,3,4,4,5,6,1,2,3]
STri=SymTridiagonal([1,2,3,4],[5,6,7])
@test has_sparsestruct(STri)
rowind,colind=findstructralnz(STri)
@test [STri[rowind[i],colind[i]] for i in 1:length(rowind)]==[1,2,3,4,5,6,7,5,6,7]

Sp=sparse([1,2,3],[1,2,3],[1,2,3])
@test has_sparsestruct(Sp)
rowind,colind=findstructralnz(Sp)
@test [Tri[rowind[i],colind[i]] for i in 1:length(rowind)]==[1,2,3]

@test ArrayInterface.ismutable(spzeros(1, 1))
@test ArrayInterface.ismutable(spzeros(1))

@test !fast_scalar_indexing(qr(rand(10, 10)).Q)
@test !fast_scalar_indexing(qr(rand(10, 10), Val(true)).Q)
@test !fast_scalar_indexing(lq(rand(10, 10)).Q)

using BandedMatrices

B=BandedMatrix(Ones(5,5), (-1,2))
B[band(1)].=[1,2,3,4]
B[band(2)].=[5,6,7]
@test has_sparsestruct(B)
rowind,colind=findstructralnz(B)
@test [B[rowind[i],colind[i]] for i in 1:length(rowind)]==[5,6,7,1,2,3,4]
B=BandedMatrix(Ones(4,6), (-1,2))
B[band(1)].=[1,2,3,4]
B[band(2)].=[5,6,7,8]
rowind,colind=findstructralnz(B)
@test [B[rowind[i],colind[i]] for i in 1:length(rowind)]==[5,6,7,8,1,2,3,4]

using BlockBandedMatrices
BB=BlockBandedMatrix(Ones(10,10),([1,2,3,4],[4,3,2,1]),(1,0))
BB[Block(1,1)].=[1 2 3 4]
BB[Block(2,1)].=[5 6 7 8;9 10 11 12]
rowind,colind=findstructralnz(BB)
@test [BB[rowind[i],colind[i]] for i in 1:length(rowind)]==
    [1,5,9,2,6,10,3,7,11,4,8,12,
    1,1,1,1,1,1,1,1,1,1,1,1,1,1,
    1,1,1,1,1,1,1,1,1,1,1,1,1,1,
    1,1,1,1,1]

dense=collect(Ones(8,8))
for i in 1:8
    dense[:,i].=[1,2,3,4,5,6,7,8]
end
BBB=BandedBlockBandedMatrix(dense, ([4, 4] ,[4, 4]), (1, 1), (1, 1))
rowind,colind=findstructralnz(BBB)
@test [BBB[rowind[i],colind[i]] for i in 1:length(rowind)]==
    [1,2,3,1,2,3,4,2,3,4,5,6,7,5,6,7,8,6,7,8,
     1,2,3,1,2,3,4,2,3,4,5,6,7,5,6,7,8,6,7,8]

@testset "setindex" begin
    @testset "$(typeof(x))" for x in [
        zeros(3),
        falses(3),
        spzeros(3),
    ]
        y = setindex(x, true, 1)
        @test iszero(x)  # x is not mutated
        @test y[1] == true
        @test iszero(x[CartesianIndices(size(x)) .== [CartesianIndex(1)]])

        y2 = setindex(x, one.(x), :)
        @test iszero(x)
        @test all(isone, y2)
    end

    @testset "$(typeof(x))" for x in [
        zeros(3, 3),
        falses(3, 3),
        spzeros(3, 3),
    ]
        y = setindex(x, true, 1, 1)
        @test iszero(x)  # x is not mutated
        @test y[1, 1] == true
        @test iszero(x[CartesianIndices(size(x)) .== [CartesianIndex(1, 1)]])

        y2 = setindex(x, one.(x), :, :)
        @test iszero(x)
        @test all(isone, y2)
    end

    @testset "$(typeof(x))" for x in [
        zeros(3, 3, 3),
        falses(3, 3, 3),
    ]
        y = setindex(x, true, 1, 1, 1)
        @test iszero(x)  # x is not mutated
        @test y[1, 1, 1] == true
        @test iszero(x[CartesianIndices(size(x)) .== [CartesianIndex(1, 1, 1)]])

        y2 = setindex(x, one.(x), :, :, :)
        @test iszero(x)
        @test all(isone, y2)
    end
end


text:
module MLPMod

using Flux, Flux.Optimise, Flux.Losses
using StatsBase, StatsModels
using ..DataMod, ..Helpers

export MLPModel
export MLP, trainMLP

struct MLPModel <: RegressionModel
    mlp
    data::Data
    data_test::Data
    best_epoch::Int
    total_epochs::Int
    params::Vector
    loss::Float64
    stopped::Symbol
end

Base.display(m::MLPModel) = Base.display(m.mlp)
StatsBase.predict(m::MLPModel, A::Matrix{Float64}) = Matrix(size(params(m.mlp)[1], 2) == size(A, 1) ? m.mlp(A)' : m.mlp(A')')
StatsBase.predict(m::MLPModel) = Matrix(m.mlp(m.data.xmat')')

function MLP(nodes...)
    length(nodes) < 2 && throw("MLP needs at least an input size and an output size.")
    if nodes[end] > 1
        Chain([[Dense(nodes[i], nodes[i+1], σ) for i in eachindex(nodes)[1:(end-2)]]; Dense(nodes[end-1], nodes[end]); softmax]...)
    else
        Chain([[Dense(nodes[i], nodes[i+1], σ) for i in eachindex(nodes)[1:(end-2)]]; Dense(nodes[end-1], nodes[end])]...)
    end
end
MLP(nodes::Vector{Int}) = MLP(nodes...)

function progress(striploss::Vector{Float64})
    k = length(striploss)
    num = sum(striploss)
    denom = k * minimum(striploss)
    return 1000 * (num / denom - 1)
end

function trainMLP(data::Data, data_test::Data; s::Int = 4, k::Int = 5,
                  stop_progress::Float64 = 0.1, max_epoch::Int = 3000, verbose::Bool = true)::MLPModel
    train = zip(eachrow(data.xmat), eachrow(data.ymat))
    nin = length(data.xs)
    nout = length(data.ys)
    h1 = min(2 * nin, 100)
    h2 = min(floor(Int, 1.3 * nin), 65)
    mlp = MLP(nin, h1, h2, nout)
    if verbose
        println("architecture of MLP:")
        display(mlp)
    end
    loss(x, y) = mse(mlp(x), y)
    ps = params(mlp)
    xtrn, ytrn = data.xmat', data.ymat'
    xtst, ytst = data_test.xmat', data_test.ymat'
    ltrn, ltst = loss(xtrn, ytrn), loss(xtst, ytst)
    pltrn, pltst = xdecim(ltrn, 5), xdecim(ltst, 5)
    verbose && println("[Before training] ~$pltrn training loss; ~$pltst testing loss")
    lastloss = Inf; csc_gain = 0; striploss = Float64[]
    bestloss = Inf; bestprms = []; bestepoch = 0
    total_epochs = 0
    stopped = :maxiter
    for epoch in 1:max_epoch
        total_epochs = epoch
        verbose && print("[Epoch $epoch] ")
        train!(loss, ps, train, ADAM())
        ltrn, ltst = loss(xtrn, ytrn), loss(xtst, ytst)
        pltrn, pltst = xdecim(ltrn, 5), xdecim(ltst, 5)
        verbose && println("~$pltrn training loss; ~$pltst testing loss")
        if ltst < bestloss
            bestprms = collect(params(cpu(mlp)))
            bestloss = ltst
            bestepoch = epoch
        end
        push!(striploss, ltrn)
        if epoch % k == 0
            if ltst > lastloss
                csc_gain += 1
                if csc_gain >= s
                    verbose && println("-- testing loss increased in $s consecutive strips")
                    stopped = :consecincrease
                    break
                elseif verbose
                    println("-- testing loss increased this strip")
                end
            else
                csc_gain = 0
            end
            lastloss = ltst
            if progress(striploss) < stop_progress
                verbose && println("-- progress this strip was < $stop_progress")
                stopped = :noprogress
                break
            end
            empty!(striploss)
        end
    end
    verbose && println("returning model with lowest testing loss ($bestloss) found on epoch $bestepoch")
    Flux.loadparams!(mlp, bestprms)
    return MLPModel(mlp, data, data_test, bestepoch, total_epochs, bestprms, bestloss, stopped)
end

end

using .MLPMod

max_stars_repo_name: KasperNooteboom/thesis-rvfl-fs

max_issues_repo_name: KasperNooteboom/thesis-rvfl-fs

max_forks_repo_name: KasperNooteboom/thesis-rvfl-fs


text:
# Reference article: [1]W.-T. Wu, Y.-H. Liao, and F.-Y. Lin, “Noise suppressions in synchronized chaos lidars,” Opt. Express, vol. 18, no. 25, p. 26155, Dec. 2010, doi: 10.1364/OE.18.026155.

using DifferentialEquations
using Plots
using FFTW

const c = 299792458.0
const ω_T = 2 * pi * c / 1.55e-6
const ω_R = ω_T

# laser parameters
const b = 4  # line width enhancement factor
const γ_n = 0.667e9  # differential carrier relaxation rate
const γ_p = 1.2e9  # nonlinear carrier relaxation rate
const γ_s = 1.458e9  # spontaneous carrier relaxation rate
const γ_c = 2.4e11  # cavity decay rate
const J = 0.333  # normalized dimensionless injection current parameter
const Δω = 0  # angular frequency detuning between the Tx and Rx

function s_clidar_model(du, u, h, p, t)
  η_OF_T, η_OF_R, η_OF_C, η_OEF_T, η_OEF_R, η_OEF_C, τ_T, τ_R, τ_C, σ_1, m = p

  a_T, ϕ_T, n_T, a_R, ϕ_R, n_R = u

  tmp1_T = γ_c * γ_n / γ_s / J * n_T - γ_p * (2a_T + a_T^2)
  # tmp1_R = γ_c * γ_n / γ_s / J * n_R - γ_p * (2a_R + a_R^2)

  # a_C_delayed = h(p, t - τ_C)[1] + σ_1 * randn()
  # ϕ_C_delayed = h(p, t - τ_C)[2] + m * (2π * rand() - π)

  # channel_a = η_OF_C * (1 + a_C_delayed) * cos(ϕ_C_delayed - ϕ_R + ω_T * τ_C)

  # channel_ϕ = η_OF_C * (1 + a_C_delayed) / (1 + a_R) * sin(ϕ_C_delayed - ϕ_R + ω_T * τ_C)

  # channel_n = η_OEF_C * γ_s * (J + 1) * (1 + 2a_C_delayed) + a_C_delayed^2

  du[1] = 0.5tmp1_T*(1 + a_T) + η_OF_T*(1 + h(p, t - τ_T)[1]) * cos(h(p, t - τ_T)[2] - ϕ_T + ω_T*τ_T)

  du[2] = -0.5b*tmp1_T + η_OF_T*(1 + h(p, t - τ_T)[1])/(1 + a_T)*sin(h(p, t - τ_T)[2] - ϕ_T + ω_T*τ_T)

  du[3] = -γ_s*n_T - γ_n*(1 + a_T)^2*n_T - γ_s*J*(2a_T + a_T^2) + γ_s*γ_p/γ_c*J*(2a_T + a_T^2)*(1 + a_T)^2 + η_OEF_T * γ_s * (J + 1) * (1 + 2h(p, t - τ_T)[1]) + h(p, t - τ_T)[1]^2

  # du[4] = 0.5tmp1_R * (1 + a_R) + η_OF_R * (1 + h(p, t - τ_R)[4]) * cos(h(p, t - τ_R)[5] - ϕ_R + ω_R * τ_R) + channel_a

  # du[5] = -0.5b * tmp1_R + η_OF_R * (1 + h(p, t - τ_R)[4]) / (1 + a_R) * sin(h(p, t - τ_R)[5] - ϕ_R + ω_R * τ_R) + channel_ϕ

  # du[6] = -γ_s * n_R - γ_n * (1 + a_R)^2 * n_R - γ_s * J * (2a_R + a_R^2) + γ_s * γ_p / γ_c * J * (2a_R + a_R^2) * (1 + a_R)^2 + η_OEF_R * γ_s * (J + 1) * (1 + 2h(p, t - τ_R)[4]) + h(p, t - τ_R)[4]^2 + channel_n

end

τ_T = 9.5e-9
τ_R = 9.5e-9
τ_C = 15.5e-9
lags = [τ_T, τ_R, τ_C]

η_OF_T, η_OF_R, η_OF_C = 0.2, 0.2, 1.3
η_OEF_T, η_OEF_R, η_OEF_C = 0, 0, 0
σ_1, m = 0, 0

h(p, t) = zeros(6)
u0 = [0, 0, 0, 0, 0, 0]
tspan = (0, 200e-9)
p = (η_OF_T, η_OF_R, η_OF_C, η_OEF_T, η_OEF_R, η_OEF_C, τ_T, τ_R, τ_C, σ_1, m)

prob = DDEProblem(s_clidar_model, u0, h, tspan, p; constant_lags = lags)
alg = MethodOfSteps(Tsit5())
@time sol = solve(prob, alg; dt=0.1e-12, adaptive=true)

fs = (sol.t[end] - sol.t[1]) / (length(sol.t) - 1)
freqs = fftfreq(length(sol.t), fs)
freqs = fftshift(freqs)

a_T = Float64.([u[1] for u in sol.u])
a_R = Float64.([u[4] for u in sol.u])
A_T = fftshift(fft(a_T))
A_R = fftshift(fft(a_R))

p1 = plot(sol.t, a_T, label="\$a_T\$")
p2 = plot(sol.t, a_R, label="\$a_R\$")
p3 = plot(freqs, abs.(A_T), label="\$\\hat{a}_T\$")
p4 = plot(sol)
# p4 = plot(freqs, abs.(A_R), label="\$\\hat{a}_R\$")

plt = plot(p1, p2, p3, p4, layout=(2, 2))

savefig("synchronized_chaos_lidar.png")

display(plt)

max_stars_repo_name: PascalZh/Optics-Simulations

max_issues_repo_name: PascalZh/Optics-Simulations

max_forks_repo_name: PascalZh/Optics-Simulations


text:
struct Student
    ime::String
    prezime::String
    brojIndeksa
end;

s1 = Student("Mitar","Miric",10)

A = [1,2,3,4,5,6,7,8,9,10]
function Sort(A,keys)
    rez = []
    for k = 1:length(A)
        if A[k] == keys
            push!(rez,k)
        end
    end
    return rez
end


text:
using ArgParse

include("BBTNet/models/retinaface.jl")
include("BBTNet/datasets/WIDERFACE.jl")
include("configs.jl")

function parse_cmd()
    s = ArgParseSettings(commands_are_required = false)
    @add_arg_table s begin
        "--batch_size", "-b"
            help = "Batch size to take for training."
            arg_type = Int
            default = 10
            required = false
        "--load_dir"
            help = "If there is a pretrained model, then the path of it."
            arg_type = String
            default = nothing
            required = false
        "--save_dir", "-s"
            help = "Directory path for saving a model after each epoch."
            arg_type = String
            default = "./weights/"
            required = false
        "--log_file"
            help = "Log file directory to write losses after each batch evaluation."
            arg_type = String
            default = nothing
            required = false
        "--mode", "-m"
            help = "Training mode: 0 for only baseline until FPN, 1 for full model, 2 for no cascaded structure."
            arg_type = Int
            default = 1
            required = false
        "--laterals", "-l"
            help = "How many lateral connections will be processed, either 3 or 5. 5 is needed for full model."
            arg_type = Int
            default = 5
            required = false
        "--start_epoch", "-e"
            help = "From which epoch the training will continue."
            arg_type = Int
            default = 1
            required = false
    end
    return parse_args(s)
end

function main()
    parsed_args = parse_cmd()
    mode = parsed_args["mode"]
    scale_cnt = parsed_args["laterals"]
    num_anchors = scale_cnt == 3 ? 2 : 3
    anchor_info = scale_cnt == 3 ? lat_3_anchors : lat_5_anchors

    bs = parsed_args["batch_size"]
    start_epoch = parsed_args["start_epoch"]
    load_path = parsed_args["load_dir"]; save_path = parsed_args["save_dir"]; log_path = parsed_args["log_file"];

    train_dir = wf_path * "train/"
    labels_dir = wf_labels_path * "train/"
    data = WIDER_Data(train_dir, labels_dir, train=true, shuffle=true, batch_size=bs, dtype=atype)
    print("[INFO] Data is loaded!\n")

    model = RetinaFace(
        mode=mode, num_anchors=num_anchors, anchor_info=anchor_info, load_path=load_path, dtype=atype
    )
    print("[INFO] Model is loaded!\n")

    train_model(model, data, save_dir=save_path, start_epoch=start_epoch, log_file=log_path)
end

main()

max_stars_repo_name: barisbatuhan/retinaface

max_issues_repo_name: barisbatuhan/retinaface

max_forks_repo_name: barisbatuhan/retinaface


text:
using Metropolis
using Test

@testset "Metropolis.jl" begin
    include("test_interface.jl")
    include("typetests.jl")
    include("simulationtests.jl")
end


text:
using CytoscapeJS
using Test

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


text:
"""
    SubspaceQs(single_space_strategy::QueryStrategy, combination_fct::Function, scale_fct::Function)

    Query Strategy for informativenss across multiple subspaces.

# Arguments
    - `single_space_strategy`: QueryStrategy that returns informativeness in a subspace
    - `scale_fct`: monotonous scaling function to scale scores within subspace
    - `combination_fct`: reduces subspaces scores to combine scores across subspaces

# Examples
    ```julia-repl
    julia> SubspaceQs(RandomPQs(), [[1,2], [3,4]], +, identity)
    ```
"""
struct SubspaceQs{Q <: QueryStrategy} <: SubspaceQueryStrategy
    single_space_strategy::Q
    subspaces::Vector{Vector{Int}}
    combination_fct::Function
    scale_fct::Function
end

"""
    Consumes params for SubspaceQs and passes the remaining
    params through to initialize the single_space_strategy.
"""
function SubspaceQs{Q}(model,
                    data::Array{T, 2} where T <: Real;
                    subspaces,
                    combination_fct = +,
                    scale_fct = identity,
                    params...) where Q <: QueryStrategy
    params = Dict(params)
    if length(params) == 0
        params = Dict{Symbol, Any}()
    end
    qs = initialize_qs(Q, model, data, params)
    SubspaceQs(qs, subspaces, combination_fct, scale_fct)
end

function qs_score(qs::SubspaceQs,
                  x::Array{T,2},
                  labels::Dict{Symbol, Array{Int, 1}}) where {T <: Real}
    size(x,1) >= maximum(maximum.(qs.subspaces)) || throw(DimensionMismatch("Maximum subspace dimension is
        larger than numver of dimensions in query data x."))
    subspace_scores = qs_score(qs.single_space_strategy, x, labels, qs.subspaces)
    try
        return mapreduce(qs.scale_fct, qs.combination_fct, subspace_scores)
    catch err
        error(getlogger(@__MODULE__), "Failed to reduce subspaces_scores with scale_fct = $(qs.scale_fct) and combination_fct = $(qs.combination_fct).
            Make sure that you use a function that is applicable to reduce (e.g., '+' and not 'sum'). \n $(sprint(io -> showerror(io, err)))")
    end
end

function qs_score(qs::QueryStrategy,
          x::Array{T, 2},
          labels::Dict{Symbol, Array{Int, 1}},
          subspaces::Vector{Vector{Int}}) where {T <: Real}
    error("Query strategy $(typeof(qs)) is not implemented for subspace classifiers.")
end


text:
type DBNSimParams
    sampling_scheme::AbstractSampleMethod
    smoothing::Symbol # :none, :SMA, :WMA
    smoothing_counts::Int    # number of previous counts to use

    function DBNSimParams(
        sampling_scheme::AbstractSampleMethod=SAMPLE_UNIFORM,
        smoothing::Symbol=:none,
        smoothing_counts::Int=1
        )

        @assert(smoothing_counts > 0)
        new(sampling_scheme, smoothing, smoothing_counts)
    end
end
type DynamicBayesianNetworkBehavior <: AbstractVehicleBehavior

    model         :: DBNModel

    ind_lat       :: Int
    ind_lon       :: Int
    symbol_lat    :: Symbol
    symbol_lon    :: Symbol
    simparams_lat :: DBNSimParams
    simparams_lon :: DBNSimParams

    # TODO(tim): remove indicators once we switch completely to RunLogs
    indicators    :: Vector{AbstractFeature}
    extractor     :: FeatureSubsetExtractor
    ordering      :: Vector{Int}
    action_clamper:: DataClamper

    # preallocated memory
    observations  :: Dict{Symbol,Float64}
    assignment    :: Dict{Symbol,Int}
    logPs         :: Dict{Symbol,Float64}
    ind_lat_in_discretizers :: Int
    ind_lon_in_discretizers :: Int
    binprobs_lat  :: Vector{Float64}
    binprobs_lon  :: Vector{Float64}
    temp_binprobs_lat :: Vector{Float64}
    temp_binprobs_lon :: Vector{Float64}

    function DynamicBayesianNetworkBehavior(
        model::DBNModel,
        simparams_lat::DBNSimParams, #OTOD(tim): get rid of these?
        simparams_lon::DBNSimParams,
        action_clamper_orig::DataClamper,
        )

        retval = new()
        retval.model = model

        targets = get_targets(model)

        f_lat = get_target_lat(model, targets)
        retval.ind_lat = indexof(f_lat, model)
        retval.symbol_lat = symbol(f_lat)

        f_lon = get_target_lon(model, targets)
        retval.ind_lon = indexof(f_lon, model)
        retval.symbol_lon = symbol(f_lon)

        retval.simparams_lat = simparams_lat
        retval.simparams_lon = simparams_lon

        retval.indicators = get_indicators(model)
        retval.ordering = topological_sort_by_dfs(model.BN.dag)
        retval.observations = Dict{Symbol,Float64}()
        retval.assignment = Dict{Symbol,Int}()
        retval.logPs = Dict{Symbol,Float64}()

        for f in retval.indicators
            sym = symbol(f)
            retval.observations[sym] = NaN
            retval.assignment[sym] = 0
        end
        retval.logPs[retval.symbol_lat] = NaN
        retval.logPs[retval.symbol_lon] = NaN

        retval.ind_lat_in_discretizers = indexof(retval.symbol_lat, model)
        retval.ind_lon_in_discretizers = indexof(retval.symbol_lon, model)

        old_way = findfirst([node.name for node in model.BN.nodes], retval.symbol_lon)
        @assert(retval.ind_lon_in_discretizers == old_way)

        retval.binprobs_lat = Array(Float64, nlabels(model.discretizers[retval.ind_lat_in_discretizers]))
        retval.binprobs_lon = Array(Float64, nlabels(model.discretizers[retval.ind_lon_in_discretizers]))
        retval.temp_binprobs_lat = deepcopy(retval.binprobs_lat)
        retval.temp_binprobs_lon = deepcopy(retval.binprobs_lon)

        x = Array(Float64, length(retval.indicators))
        retval.extractor = FeatureSubsetExtractor(x, retval.indicators)

        retval.action_clamper = DataClamper(Array(Float64, 2),
                               deepcopy(action_clamper_orig.f_lo),
                               deepcopy(action_clamper_orig.f_hi))

        retval
    end
end

function Base.print(io::IO, BN::DynamicBayesianNetworkBehavior)
    println(io, "DynamicBayesianNetworkBehavior")
    println(io, "\tDBNModel: ")
    print(io, BN.model)
    println(io, "\tind_lat:          ", BN.ind_lat)
    println(io, "\tind_lon:          ", BN.ind_lon)
    println(io, "\tsymbol_lat:       ", BN.symbol_lat)
    println(io, "\tsymbol_lon:       ", BN.symbol_lon)
    println(io, "\tordering:         ", BN.ordering)
    println(io, "\tsubset_extractor: ", BN.extractor)
    println(io, "\taction_clamper:   ", BN.action_clamper)
end
AutomotiveDrivingModels.get_targets(behavior::DynamicBayesianNetworkBehavior) = ModelTargets(symbol2feature(behavior.symbol_lat), Features.symbol2feature(behavior.symbol_lon))
AutomotiveDrivingModels.get_indicators(behavior::DynamicBayesianNetworkBehavior) = behavior.extractor.indicators

# function infer_action_lon_from_input_acceleration(sym::Symbol, accel::Float64, simlog::Matrix{Float64}, frameind::Int, logindexbase::Int)

#     if sym == :f_accel_250ms || sym == :f_accel_500ms
#         return accel
#     elseif sym == :f_des_speed_250ms || sym == :f_des_speed_500ms
#         return accel/Features.KP_DESIRED_SPEED
#     else
#         error("unknown longitudinal target $sym")
#     end
# end
# function infer_action_lat_from_input_turnrate(sym::Symbol, turnrate::Float64, simlog::Matrix{Float64}, frameind::Int, logindexbase::Int)

#     if sym == :f_turnrate_250ms || sym == :f_turnrate_500ms
#         return turnrate
#     elseif sym == :f_des_angle_250ms || sym == :f_des_angle_500ms
#         ϕ = simlog[frameind, logindexbase + LOG_COL_ϕ]
#         return (turnrate / Features.KP_DESIRED_ANGLE) + ϕ
#     else
#         error("unknown lateral target $sym")
#     end
# end

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

sample!(behavior::DynamicBayesianNetworkBehavior, assignment::Dict{Symbol, Int}) = sample!(behavior.model, assignment, behavior.ordering)
sample_and_lopP!(behavior::DynamicBayesianNetworkBehavior, assignment::Dict{Symbol, Int}, logPs::Dict{Symbol, Float64}=Dict{Symbol, Float64}()) = sample!(behavior.model, assignment, logPs, behavior.ordering)

function _copy_extracted_into_obs!(behavior::DynamicBayesianNetworkBehavior)

    extractor = behavior.extractor

    # copy them over into the observation dict
    for (i,f) in enumerate(extractor.indicators)
        sym = symbol(f)
        behavior.observations[sym] = extractor.x[i]
    end

    behavior
end
function _set_and_process_action!(behavior::DynamicBayesianNetworkBehavior, action_lat::Float64, action_lon::Float64)
    behavior.action_clamper.x[1] = action_lat
    behavior.action_clamper.x[2] = action_lon
    process!(behavior.action_clamper)
    behavior.action_clamper
end

function select_action(
    behavior::DynamicBayesianNetworkBehavior,
    runlog::RunLog,
    sn::StreetNetwork,
    colset::UInt,
    frame::Int
    )

    model = behavior.model
    symbol_lat = behavior.symbol_lat
    symbol_lon = behavior.symbol_lon
    extractor = behavior.extractor

    simparams_lat = behavior.simparams_lat
    simparams_lon = behavior.simparams_lon
    samplemethod_lat = simparams_lat.sampling_scheme
    samplemethod_lon = simparams_lon.sampling_scheme

    bmap_lat = model.discretizers[behavior.ind_lat_in_discretizers]
    bmap_lon = model.discretizers[behavior.ind_lon_in_discretizers]

    observations = behavior.observations
    assignment = behavior.assignment

    observe!(extractor, runlog, sn, colset, frame)
    _copy_extracted_into_obs!(behavior)

    encode!(assignment, model, observations)

    assignment[symbol_lat] = -1
    assignment[symbol_lon] = -1
    sample_unset!(model, assignment, behavior.ordering)

    bin_lat = assignment[symbol_lat]
    bin_lon = assignment[symbol_lon]

    behavior.action_clamper.x[1] = decode(bmap_lat, bin_lat, samplemethod_lat)
    behavior.action_clamper.x[2] = decode(bmap_lon, bin_lon, samplemethod_lon)
    process!(behavior.action_clamper)
    action_lat = behavior.action_clamper.x[1]
    action_lon = behavior.action_clamper.x[2]

    (action_lat, action_lon)
end

function _calc_action_loglikelihood(behavior::DynamicBayesianNetworkBehavior)

    #=
    Actions should be in behavior.action_clamper
    Observations should be in behavior.observations
    =#

    model = behavior.model
    symbol_lat = behavior.symbol_lat
    symbol_lon = behavior.symbol_lon
    bmap_lat = model.discretizers[behavior.ind_lat_in_discretizers]
    bmap_lon = model.discretizers[behavior.ind_lon_in_discretizers]

    bin_lat = encode(bmap_lat, behavior.action_clamper.x[1])
    bin_lon = encode(bmap_lon, behavior.action_clamper.x[2])

    observations = behavior.observations # assumed to already be populated
    assignment   = behavior.assignment   # this will be overwritten
    logPs        = behavior.logPs        # this will be overwritten
    binprobs_lat = behavior.binprobs_lat # this will be overwritten
    binprobs_lon = behavior.binprobs_lon # this will be overwritten

    encode!(assignment, model, observations)

    if is_parent(model, symbol_lon, symbol_lat) # lon -> lat
        calc_probability_distribution_over_assignments!(binprobs_lon, model, assignment, symbol_lon)
        fill!(binprobs_lat, 0.0)
        temp = behavior.temp_binprobs_lat
        for (i,p) in enumerate(binprobs_lon)
            assignment[symbol_lon] = i
            calc_probability_distribution_over_assignments!(temp, model, assignment, symbol_lat)
            for (j,v) in enumerate(temp)
                binprobs_lat[j] +=  v * p
            end
        end
    elseif is_parent(model, symbol_lat, symbol_lon) # lat -> lon
        calc_probability_distribution_over_assignments!(binprobs_lat, model, assignment, symbol_lat)
        fill!(binprobs_lon, 0.0)
        temp = behavior.temp_binprobs_lon
        for (i,p) in enumerate(binprobs_lat)
            assignment[symbol_lat] = i
            calc_probability_distribution_over_assignments!(temp, model, assignment, symbol_lon)
            for (j,v) in enumerate(temp)
                binprobs_lon[j] +=  v * p
            end
        end
    else # lat and lon are conditionally independent
        calc_probability_distribution_over_assignments!(binprobs_lat, model, assignment, symbol_lat)
        calc_probability_distribution_over_assignments!(binprobs_lon, model, assignment, symbol_lon)
    end

    P_bin_lat = binprobs_lat[bin_lat]
    P_bin_lon = binprobs_lon[bin_lon]

    p_within_bin_lat = calc_probability_for_uniform_sample_from_bin(P_bin_lat, bmap_lat, bin_lat)
    p_within_bin_lon = calc_probability_for_uniform_sample_from_bin(P_bin_lon, bmap_lon, bin_lon)

    # println("actions: ", action_lat, "  ", action_lon)

    log(p_within_bin_lat) + log(p_within_bin_lon)
end
function calc_action_loglikelihood(
    behavior::DynamicBayesianNetworkBehavior,
    runlog::RunLog,
    sn::StreetNetwork,
    colset::UInt,
    frame::Int,
    action_lat::Float64,
    action_lon::Float64,
    )

    observe!(behavior.extractor, runlog, sn, colset, frame)
    _copy_extracted_into_obs!(behavior)
    _set_and_process_action!(behavior, action_lat, action_lon)

    _calc_action_loglikelihood(behavior)
end
function calc_action_loglikelihood(
    behavior::DynamicBayesianNetworkBehavior,
    features::DataFrame,
    frame::Integer,
    )

    action_lat = features[frame, behavior.symbol_lat]::Float64
    action_lon = features[frame, behavior.symbol_lon]::Float64

    observe!(behavior.extractor, features, frame)
    _copy_extracted_into_obs!(behavior)
    _set_and_process_action!(behavior, action_lat, action_lon)

    _calc_action_loglikelihood(behavior)
end


text:
using StochasticPrograms
using LShapedSolvers
using ProgressiveHedgingSolvers
using LaTeXStrings
using Gurobi
using Plots
pyplot()

KTH_colors = [RGB(25/255,84/255,166/255),
              RGB(157/255,16/255,45/255),
              RGB(98/255,146/255,46/255),
              RGB(36/255,160/255,216/255),
              RGB(228/255,54/255,62/255),
              RGB(176/255,201/255,43/255),
              RGB(216/255,84/255,151/255),
              RGB(250/255,185/255,25/255),
              RGB(101/255,101/255,108/255),
              RGB(189/255,188/255,188/255),
              RGB(227/255,229/255,227/255)]

# Model
@info "Constructing farmer problem"
Crops = [:wheat,:corn,:beets]
Purchased = [:wheat,:corn]
Sold = [:wheat, :corn, :beets_quota, :beets_extra]
Cost = Dict(:wheat=>150,:corn=>230,:beets=>260)
Required = Dict(:wheat=>200,:corn=>240,:beets=>0)
PurchasePrice = Dict(:wheat=>238,:corn=>210)
SellPrice = Dict(:wheat=>170, :corn=>150, :beets_quota=>36, :beets_extra=>10)
Budget = 500
@scenario Yield = begin
    wheat::Float64
    corn::Float64
    beets::Float64
end
ξ₁ = YieldScenario(3.0, 3.6, 24.0, probability = 1/3)
ξ₂ = YieldScenario(2.5, 3.0, 20.0, probability = 1/3)
ξ₃ = YieldScenario(2.0, 2.4, 16.0, probability = 1/3)
farmer_model = StochasticModel((Crops,Cost,Budget), (Required,PurchasePrice,SellPrice), (sp)->begin
    @first_stage sp = begin
        (Crops,Cost,Budget) = stage
        @variable(model, x[c = Crops] >= 0)
        @objective(model, Min, sum(Cost[c]*x[c] for c in Crops))
        @constraint(model, sum(x[c] for c in Crops) <= Budget)
    end
    @second_stage sp = begin
        @decision x
        (Required, PurchasePrice, SellPrice) = stage
        ξ = scenario
        @variable(model, y[p = Purchased] >= 0)
        @variable(model, w[s = Sold] >= 0)
        @objective(model, Min, sum( PurchasePrice[p] * y[p] for p = Purchased) - sum( SellPrice[s] * w[s] for s in Sold))

        @constraint(model, const_minreq[p=Purchased],
            ξ[p] * x[p] + y[p] - w[p] >= Required[p])
        @constraint(model, const_minreq_beets,
            ξ[:beets] * x[:beets] - w[:beets_quota] - w[:beets_extra] >= Required[:beets])
        @constraint(model, const_aux, w[:beets_quota] <= 6000)
    end
end)
farmer_problem = instantiate(farmer_model, [ξ₁,ξ₂,ξ₃])
gurobi = GurobiSolver(OutputFlag=0)
optimize!(farmer_problem, solver = gurobi)
Q = optimal_value(farmer_problem)
@info "Solving farmer problem using progressive-hedging"
sm = StochasticPrograms.StructuredModel(farmer_problem, ProgressiveHedgingSolver(gurobi))
sm()
@info "Plotting results..."
p = plot(sm.Q_history./Q, color_palette = KTH_colors, label="Progressive-hedging", linewidth = 4, xlabel = "Iterations", ylabel = L"\frac{Q_k}{Q^*}")
savefig(p, "farmer_ph.pdf")
@info "Solving farmer problem using L-shaped"
sm = StochasticPrograms.StructuredModel(farmer_problem, LShapedSolver(gurobi))
sm()
@info "Plotting results..."
p = plot(sm.Q_history./Q, label="L-shaped", color_palette = KTH_colors, linewidth = 4, xlabel = "Iterations", ylabel = L"\frac{Q_k}{Q^*}")
savefig(p, "farmer_ls.pdf")

max_stars_repo_name: martinbiel/mpc-reproducibility

max_issues_repo_name: martinbiel/mpc-reproducibility

max_forks_repo_name: martinbiel/mpc-reproducibility


text:
# Base for implementing all algorithms
using Printf
"""
    BanditAlgorithmBase

`BanditAlgorithmBase` defines a set of functions as well as some utility functions for
implementing MAB algorithms. Any algorithm must be a subtype of `BanditAlgorithmBase`.

"""
abstract type BanditAlgorithmBase end

"""
    get_arm_index( agent::BanditAlgorithmBase )

Gets the index of next arm to pull.

## Example
```julia-repl

```
"""

function add_arms!( agent::BanditAlgorithmBase, arms::Array{Int64} )
    error( "No implementation of add_arms() for ", typeof(agent) );
end

function get_arm_index( agent::BanditAlgorithmBase )
    error( "No implementation of getArmIndex() for ", typeof(agent) );
end

"""
    update_reward!( agent::BanditAlgorithmBase, r::Real )

Updates the reward to bandit algorithm `agent`. `r` must a real number within valid range.
"""
function update_reward!( agent::BanditAlgorithmBase, r::Real )
    error( "No implementation of updateReward() for ", typeof(agent) );
end

"""
    reset!( agent::BanditAlgorithmBase )

Resets the internal statistics of the bandit algorithm, except for the number of arms and other
algorithm specific parameters. `reset!()` will make `agent` as fresh as it was first created.
"""
function reset!( agent::BanditAlgorithmBase )
    error( "No implementation of reset() for ", typeof(agent) );
end

"""
    info_str( agent::BanditAlgorithmBase, latex::Bool = false )

Return a information string about the agent and it's parameters.  If flag `latex` is set to
`true`, then the returned string will be a compactable latex string.
"""
function info_str( agent::BanditAlgorithmBase, latex::Bool = false )
    error( "No Implementation of info_str() for ", typeof(agent) );
end

# import Base.show
function Base.show( io::IO, ::MIME"text/plain", agent::BanditAlgorithmBase )
    print( io, @sprintf("Algorithm: %s",info_str(agent)) )
    for param in fieldnames(typeof(agent))
        print( @sprintf("\n    %-16s: ",param), getfield(agent,param) )
    end
end

# Functions to make instances of types
function call( ::Type{T}, params... ) where T <: BanditAlgorithmBase
    return T( params... )
end

"""
    make_agents_with_k( K::Int64, agent_list::Vector{} )

Returns a vector of agents specified in `agent_list` with `K` arms.
```julia-repl

```
"""
function make_agents_with_k( K::Int64, agent_list::Vector{} )
    agents = Vector{BanditAlgorithmBase}()
    for (agent,params) ∈ agent_list
        if !isnothing(params)
            push!( agents, agent(K,params...) )
        else
            push!( agents, agent(K) )
        end
    end
    return agents
end

import Base: ==
"""
    ==

Can be used to compare two agents.
```julia-repl

```
"""
function ==( agent1::BanditAlgorithmBase, agent2::BanditAlgorithmBase )
    if typeof(agent1) != typeof(agent2)
        return false
    else
        for param in fieldnames(agent1)
            if getfield(agent1,param) != getfield(agent2,param)
                # info( "Failed at ", param )
                return false
            end
        end
        return true
    end
    false
end


text:
#=
Contains info specific to finite volume such as
- cell enters
- face centers
- interpolation weights
-

Also has functions useful for FV such as
- cell to node / node to cell interpolation
-
=#

struct FVInfo
    fluxOrder::Int                  # Order for flux reconstruction
    cellCenters::Array{Float64}     # Coordinates of cell centers
    faceCenters::Array{Float64}     # Coordinates of face centers

    cell2node::Array{Array{Int64,1}}            # Neighboring cell indices for each node
    cell2nodeWeight::Array{Array{Float64,1}}    # Cell to node interpolation weights
end

etypetonv = [2, 3, 4, 4, 8, 6, 5, 2, 3, 4, 4, 8, 6, 5, 1, 4, 8, 6, 5]; # number of vertices for each type
etypetodim= [1, 2, 2, 3, 3, 3, 3, 1, 2, 2, 3, 3, 3, 3, 1, 2, 2, 3, 3]; # dimension of each type
etypetonf = [2, 3, 4, 4, 6, 5, 5, 2, 3, 4, 4, 6, 5, 5, 1, 4, 6, 5, 5]; # number of faces for element types
etypetoftype=[0,1, 1, 2, 3, 3, 3, 0, 1, 1, 2, 3, 3, 3, 0, 0, 0, 0, 0]; # type of faces for this element type

# Build the FVInfo from a given grid
function build_FV_info(grid, order=1)
    nel = size(grid.loc2glb, 2);
    nface = size(grid.face2glb, 3);
    nnode = size(grid.allnodes, 2);
    dim = size(grid.allnodes, 1);

    cell_centers = zeros(dim, nel);
    face_centers = zeros(dim, nface);
    cell2node = Array{Array{Int64,1},1}(undef,nnode);
    cell2nodeWeight = Array{Array{Float64,1},1}(undef,nnode);

    # Find cell centers
    for ei=1:nel
        center = zeros(dim);
        nvtx = size(grid.loc2glb, 1);
        for ni=1:nvtx
            center .+= grid.allnodes[:, grid.loc2glb[ni, ei]];
        end
        cell_centers[:, ei] = center ./ nvtx;
    end

    # Find face centers
    for fi=1:nface
        center = zeros(dim);
        nvtx = size(grid.face2glb, 1);
        for ni=1:nvtx
            center .+= grid.allnodes[:, grid.face2glb[ni, 1, fi]];
        end
        face_centers[:, fi] = center ./ nvtx;
    end

    # Build cell2node map
    # Check all elements to find ones touching this node.
    nvtx = size(grid.loc2glb, 1);
    for ni=1:nnode
        cell2node[ni] = zeros(Int,0);
        for ei=1:nel
            for nj=1:nvtx
                if grid.loc2glb[nj, ei] == ni
                    # This cell touches this node
                    push!(cell2node[ni], ei);
                    break;
                end
            end
        end
    end

    # Compute cell2node interpolation weights
    for ni=1:nnode
        cell2nodeWeight[ni] = zeros(length(cell2node[ni]));
        sumweight = 0;
        for ei=1:length(cell2node[ni])
            dist = norm(grid.allnodes[:,ni] - cell_centers[:, cell2node[ni][ei]]);
            sumweight += 1 / dist;
            cell2nodeWeight[ni][ei] = 1 / dist;
        end
        for ei=1:length(cell2node[ni])
            cell2nodeWeight[ni][ei] /= sumweight;
        end
    end

    return FVInfo(order, cell_centers, face_centers, cell2node, cell2nodeWeight);
end

# Interpolate nodal values for one cell.
# For now this just assigns all nodes the cell average.
# This shoould be improved eventually, but needs care to be efficient.
function FV_cell_to_node(cell_value, grid)
    return fill(cell_value, size(grid.loc2glb,1));
end

# Approximates a cell average by simply averaging nodal values.
function FV_node_to_cell(node_values)
    return sum(node_values)/length(node_values);
end

# Interpolate nodal values from neighboring cells.
# Returns an array of nodal values.
function FV_cell_to_node_all(cell_values, node_values; dofs = 1)
    Nnodes = size(fv_grid.allnodes,2);
    if node_values === nothing
        node_values = zeros(Nnodes * dofs);
    end

    for ni=1:Nnodes
        if dofs > 1
            st = (ni-1)*dofs + 1;
            en = ni*dofs;
            node_values[st:en] .= 0;
            ncell = length(fv_info.cell2nodeWeight[ni]);
            for ci=1:ncell
                cind = fv_info.cell2node[ni][ci];
                cst = (cind-1)*dofs + 1;
                cen = cind*dofs;
                node_values[st:en] .+= fv_info.cell2nodeWeight[ni][ci] .* cell_values[cst:cen];
            end

        else
            node_values[ni] = 0;
            ncell = length(fv_info.cell2nodeWeight[ni]);
            for ci=1:ncell
                node_values[ni] += fv_info.cell2nodeWeight[ni][ci] * cell_values[fv_info.cell2node[ni][ci]];
            end
        end
    end

    return node_values;
end

# Find cell averages from nodal values
# Returns an array of cell values.
function FV_node_to_cell_all(node_values, cell_values = nothing)
    Ncells = size(fv_grid.loc2glb,2);
    if cell_values === nothing
        cell_values = zeros(Ncells);
    end

    for ci=1:Ncells
        cell_values[ci] = 0;
        nnode = length(fv_grid.loc2glb[:,ci]);
        for ni=1:nnode
            cell_values[ci] += node_values[fv_grid.loc2glb[ni, ci]];
        end
        cell_values[ci] /= nnode;
    end

    return cell_values;
end

# reconstructs u at x based on a cell set
function FV_reconstruct_value(cellx, cellu, x)
    # If only one cell given, return that value
    if length(cellu) == 1
        return cellu[1];
    end

    # Interpolate value at x
    # Note: if extrapolating from only one side,
    # consider using FV_reconstruct_value_left_right to include slope limiting
    return polyharmonic_interp(x, cellx, cellu)[1];
end

# reconstructs u at x based on left and right cell sets
function FV_reconstruct_value_left_right(leftx, rightx, leftu, rightu, x; limiter="none")
    left = 0.0;
    right = 0.0;
    leftslope = 0.0;
    rightslope = 0.0;
    centerslope = 0.0;

    # if only one cell on either side, just return those values
    # If more, extrapolate
    if length(leftu) == 1
        left = leftu[1];

    elseif length(leftu) > 1
        leftslope = (leftu[1]-leftu[2]) / norm(leftx[:,1] - leftx[:,2]);
        left = polyharmonic_interp(x, leftx, leftu)[1];
    end

    if length(rightu) == 1
        right = rightu[1];

    elseif length(rightu) > 1
        rightslope = (rightu[2]-rightu[1]) / norm(rightx[:,1] - rightx[:,2]);
        right = polyharmonic_interp(x, rightx, rightu)[1];
    end

    # if one side is empty(boundary?), just make them equal
    if length(leftu) == 0 && length(rightu) > 0
        left = right;
    elseif length(rightu) == 0 && length(leftu) > 0
        right = left;

    else
        center_dist = norm(rightx[1] - leftx[1]);
        if center_dist < 1e-16
            # This could also be a boundary.
            centerslope = 0.0;
        else
            centerslope = (rightu[1] - leftu[1]) / center_dist;
        end

    end

    # If a limiter is specified, limit the slope
    if !(limiter == "none")
        left_r = abs(leftslope) > 1e-10 ? centerslope / leftslope : 1;
        right_r = abs(rightslope) > 1e-10 ? centerslope / rightslope : 1;

        left_phi = 2;
        right_phi = 2;
        if limiter == "vanleer"
            left_phi = (left_r + abs(left_r)) / (1+left_r);
            right_phi = (right_r + abs(right_r)) / (1+right_r);
        end

        if length(leftu) > 1
            left = (1-left_phi*0.5)*leftu[1] + (left_phi*0.5)*left;
        end
        if length(rightu) > 1
            right = (1-right_phi*0.5)*rightu[1] + (right_phi*0.5)*right;
        end
    end

    return (left, right);
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
struct IMDP
    states
    actions
    Pmin
    Pmax
    labels
    sink_labels
    initial_state
end

"""
Create an IMDP from a list of result directories.
"""
function create_imdp_from_result_dirs(res_dirs, exp_dir; label_fn=nothing)
    sys_dir, minPr, maxPr = create_switched_system(res_dirs, exp_dir)
    N, M = size(minPr)
    Q = collect(1:M)
    A = collect(1:Int(N/M))
    @debug "Actions: " A
    labels = Dict()
    imdp = IMDP(Q, A, minPr, maxPr, labels, nothing, 1)
    if !isnothing(label_fn)
        create_imdp_labels(label_fn, imdp, "$exp_dir/switched-system/regions.bson")
    end
    return imdp
end

"
Create an IMDP from a set of transition bound matrices.
"
function create_simple_imdp(minPr, maxPr)
    N, _ = size(minPr)
    Q = collect(1:N)
    A = [1]
    labels = Dict()
    for i=1:N-1
        labels[i] = "safe"
    end
    labels[N] = "!safe"
    imdp = IMDP(Q, A, minPr, maxPr, labels, nothing, 1)
    return imdp
end

function create_imdp_labels(labels_fn, imdp, extent_file::String)
    extents, _ = deserialize_region_data(extent_file)
    create_imdp_labels(labels_fn, imdp, extents)
end

function create_imdp_labels(labels_fn, imdp, imdp_state_extents::Dict)
    for state_key in keys(imdp_state_extents)
        state = imdp_state_extents[state_key]
        if state_key == imdp_state_extents.count
            imdp.labels[state_key] = labels_fn(state, unsafe=true)
        else
            imdp.labels[state_key] = labels_fn(state)
        end
    end
end

function write_imdp_to_file_bounded(imdp, Qyes, Qno, filename)
    open(filename, "w") do f
        state_num = length(imdp.states)
        action_num =length(imdp.actions)
        @printf(f, "%d \n", state_num)
        @debug "Length actions: " length(imdp.actions)
        @printf(f, "%d \n", length(imdp.actions))
        # Get number of accepting states from the labels vector
        acc_states = Qyes
        @printf(f, "%d \n", length(acc_states))
        [@printf(f, "%d ", acc_state-1) for acc_state in acc_states]
        @printf(f, "\n")
        sink_states = Qno

        for i=1:state_num
            if isnothing(sink_states) || !(i∈sink_states)
                for action in imdp.actions
                    row_idx = (i-1)*action_num + action
                    ij = findall(>(0.), imdp.Pmax[(i-1)*action_num + action, :])
                    # Something about if the upper bound is less than one? Perhaps for numerical issues?
                    @debug action, i
                    psum = sum(imdp.Pmax[row_idx, :])
                    psum >= 1 ? nothing : throw(AssertionError("Bad max sum: $psum"))
                    for j=ij
                        @printf(f, "%d %d %d %f %f", i-1, action-1, j-1, imdp.Pmin[row_idx, j], imdp.Pmax[row_idx, j])
                        if (i < state_num || j < ij[end] || action < action_num)
                            @printf(f, "\n")
                        end
                    end
                end
            else
                @printf(f, "%d %d %d %f %f", i-1, 0, i-1, 1.0, 1.0)
                if i<state_num
                    @printf(f, "\n")
                end
            end
        end
    end
end

function create_dot_graph(imdp::IMDP, filename::String)
    open(filename, "w") do f
        println(f, "digraph G {")
        println(f, "  rankdir=LR\n  node [shape=\"circle\"]\n  fontname=\"Lato\"\n  node [fontname=\"Lato\"]\n  edge [fontname=\"Lato\"]")
        println(f, "  size=\"8.2,8.2\" node[style=filled,fillcolor=\"#FDEDD3\"] edge[arrowhead=vee, arrowsize=.7]")

        # Initial State
        @printf(f, "  I [label=\"\", style=invis, width=0]\n  I -> %d\n", imdp.initial_state)

        for state in imdp.states
            @printf(f, "  %d [label=<q%d>, xlabel=<%s>]\n", state, state, imdp.labels[state])

            for action in imdp.actions
                row_idx = (state-1)*length(imdp.actions) + action

                for idx in findall(>(0.), imdp.Pmax[row_idx, :])
                    state_p = idx
                    @printf(f, "  %d -> %d [label=<a%d: %.1f-%.1f >]\n", state, state_p, action, imdp.Pmin[row_idx,state_p], imdp.Pmax[row_idx,state_p])
                end
            end
        end
        println(f, "}")
    end
end

"""
Get IMDP Ring states for a gridworld
"""
function get_imdp_ring_states(imdp_state_idx, num_stats_per_col, num_imdp_states; col_ordering=true)
    # TODO: enfore column ordering
    ring_imdp_states = [imdp_state_idx + 1, imdp_state_idx - 1,
                        imdp_state_idx - num_stats_per_col, imdp_state_idx - num_stats_per_col + 1, imdp_state_idx - num_stats_per_col - 1,
                        imdp_state_idx + num_stats_per_col, imdp_state_idx + num_stats_per_col + 1, imdp_state_idx + num_stats_per_col - 1]

   # To do: Better filter

   # Remove the negative states
   ring_imdp_states = ring_imdp_states[ring_imdp_states .> 1]
   ring_imdp_states = ring_imdp_states[ring_imdp_states .< num_imdp_states]
   return ring_imdp_states
end


text:
module Weights
using RecipesBase
using ControlSystems, MonteCarloMeasurements
export makeweight, neglected_delay, gain_and_delay_uncertainty, neglected_lag
export fit_complex_perturbations

"""
    makeweight(low, [mid,] high)

Create a weighting function that goes from gain `low` at zero frequency, through gain `mid` to gain `high` at ∞

# Arguments:
- `low`: A number specifying the DC gain
- `mid`: A number specifying the frequency at which the gain is 1, or a tuple `(freq, gain)`. If `mid` is not specified, the geometric mean of `high` and `low` is used.
- `high`: A number specifying the gain at ∞
"""
function makeweight(low, mid::Number, high)
    makeweight(low, (mid, high < 1 ? sqrt(high*low) : 1), high)
end

makeweight(low, high) = makeweight(low, √(high*low), high)

function makeweight(low, mid, high)
    freq, mag = mid
    p = freq * √(((high/mag)^2-1)/(1-(low/mag)^2))
    W = tf([high, low*p], [1, p])
    W = ss(W)
end

"""
    neglected_delay(Lmax)

Return a multiplicative weight to represent the uncertainty coming from neglecting the dynamics `exp(-s*L)`
where `L ≤ Lmax`.
"Multivariable Feedback Control: Analysis and Design" Ch 7.4.5
"""
function neglected_delay(Lmax)
    gain_and_delay_uncertainty(1, 1, Lmax)
end

"""
    gain_and_delay_uncertainty(kmin, kmax, Lmax)

Return a multiplicative weight to represent the uncertainty coming from neglecting the dynamics `k*exp(-s*L)`
where `k ∈ [kmin, kmax]` and `L ≤ Lmax`.
This weight is slightly optimistic, an expression for a more exact weight appears in eq (7.27), "Multivariable Feedback Control: Analysis and Design"
"""
function gain_and_delay_uncertainty(kmin, kmax, Lmax)
    kb = (kmin+kmax)/2
    rk = (kmax-kmin)/2/kb
    tf([(1+rk/2)*Lmax, rk], [Lmax/2, 1])
end

"""
    neglected_lag(τmax)

Return a multiplicative weight to represent the uncertainty coming from neglecting the dynamics `1/(s*τ + 1)`
where `τ ≤ τmax`.
"Multivariable Feedback Control: Analysis and Design" Ch 7.4.5
"""
function neglected_lag(τmax)
    tf([τmax, 0], [τmax, 1])
end

"""
    centers, radii = fit_complex_perturbations(P, w; relative=true, nominal=:mean)

For each frequency in `w`, fit a circle in the complex plane that contains all models in the model set `P`, represented as an `LTISystem` with `Particles` coefficients. Note, the resulting radii can be quite unstable if the number of particles is small, in particular if the particles are normally distributed instead of uniformly.

If `realtive = true`, circles encompassing `|(P - Pn)/Pn|` will be returned (multiplicative/relative uncertainty).
If `realtive = false`, circles encompassing `|P - Pn|` will be returned (additive uncertainty).

If `nominal = :mean`, the mean of `P` will be used as nominal model. If `nominal = :first`, the first particle will be used. See `MonteCarloMeasurements.with_nominal` to set the nominal value in the first particle.

See also [`nyquistcircles`](@ref) to plot circles (only if relative=false).
"""
function fit_complex_perturbations(P, w; relative=true, nominal=:mean)
    ControlSystems.issiso(P) || throw(ArgumentError("This function only works for SISO systems."))
    r = freqresp(P, w)
    centers = Vector{Complex{eltype(w)}}(undef, length(w))
    radii = Vector{eltype(w)}(undef, length(w))
    for i in axes(r,1)
        ri = r[i,1,1] # NOTE: assumes SISO
        c = nominal === :mean ? pmean(ri) : MonteCarloMeasurements.nominal(ri)
        rad = abs(relative ? pmaximum((ri - c)/c) : pmaximum(ri - c))
        centers[i] = c
        radii[i] = rad
    end
    centers, radii
end

"""
    nyquistcircles(w, centers, radii)

Plot the nyquist curve with circles. It only makes sense to call this function if the circles represent additive uncertainty, i.e., were calculated with `relative=false`.

See also [`fit_complex_perturbations`](@ref)
"""
nyquistcircles
@userplot nyquistcircles
@recipe function nyquistcircles(h::nyquistcircles)
    w, centers, radii = h.args
    t = LinRange(0, 2pi, 100)
    re, im = cos.(t), sin.(t)
    @series begin
        lalbe --> "Center"
        centers
    end
    for i in eachindex(radii)
        @series begin
            alpha --> 0.3
            color --> 1
            fillcolor --> 1
            label --> ""
            re .* radii[i] .+ centers[i].re, im .* radii[i] .+ centers[i].im
        end
    end
end

end


text:
"""
Refines the specificed IMDP states using a uniform discretization of each state.
# Arguments
- `states_to_refine` - Vector of state indeces to be refined
- `params::ExperimentParameters` - Experiment parameters structure
"""
function refine_imdp(states_to_refine, params::ExperimentParameters;
    working_dir=nothing, results_dir=nothing,
    reuse_regions_flag=false, reuse_transition_mats_flag=false,
    label_fcn=nothing)

    experiment_dir = params.experiment_directory
    if isnothing(working_dir)
        working_dir = experiment_dir
    end

    if isnothing(results_dir)
        results_dir = "$experiment_dir/refined"
    end
    !isdir(results_dir) && mkpath(results_dir)

    tmat_filename = "$results_dir/transition_mats.bson"
    # ! Transition matrices being saved incorrectly. Not able to reload.
    if reuse_transition_mats_flag && isfile(tmat_filename) && false
        transition_mats = load_transition_matrices(tmat_filename)
        @goto reuse_transition_mats
    end

    # Load the GPs
    gp_filename = generate_gp_filename(params)
    _, gp_info = load_gp_info(gp_filename)

    region_filename = "$results_dir/regions.bson"
    if reuse_regions_flag && isfile(region_filename)
        region_data = load_region_data(region_filename)
        @goto reuse_regions
    end

    region_data = load_region_data("$working_dir/regions.bson")
    extents = region_data[:extents]
    post_data = region_data[:posts]

    # Generate new discretization
    new_region_ids = []
    for state_to_refine in states_to_refine
        extents, new_ids = uniform_refinement(state_to_refine, extents)
        new_region_ids = new_region_ids ∪ new_ids
    end

    # Generate new posterior bounds
    Threads.@threads for new_id in new_region_ids
        # TODO: Add a check here for "small enough" covariance bounds, and just update mean if is
        σ_ubs = calculate_σ_global_bounds(params.domain, gp_info)
        post_data[new_id] = bound_extent(extents[new_id], gp_info, params.system_params.dependency_dims, known_component=params.system_params.known_dynamics_fcn, σ_ubs=σ_ubs)
    end

    # Fix the extent ids to correspond to order
    new_extent_dict = Dict()
    new_post_dict = Dict()
    sorted_keys = sort(collect(keys(extents)))
    for (i,j) in zip(1:extents.count-1, sorted_keys[1:end-1])
        new_extent_dict[i] = extents[j]
        new_post_dict[i] = post_data[j]
    end
    new_extent_dict[new_extent_dict.count+1] = extents[sorted_keys[end]]

    # Overwrite region data dictionary
    region_data[:extents] = new_extent_dict
    region_data[:posts] = new_post_dict
    save_region_data(new_extent_dict, new_post_dict, region_filename)

    @label reuse_regions
    transition_mats = generate_transition_bounds(params, gp_info, region_data)
    save_transition_matrices(transition_mats, tmat_filename)

    @label reuse_transition_mats
    # Rebuild the IMDP through rebuilding the matrices
    new_imdp = create_simple_imdp(transition_mats["minPr"], transition_mats["maxPr"])

    validate_pimdp(new_imdp)

    if !isnothing(label_fcn)
        create_imdp_labels(label_fcn, new_imdp, new_extent_dict)
    end

    return new_imdp
end


text:
#=

  Calculs d'enfer puzzle in Julia ConstraintSolver.jl

  Problem from Jianyang Zhou "The Manual of NCL version 1.2", page 33
  http://citeseer.ist.psu.edu/161721.html

  The solution is the manual is:
  """
  a = -16, b = -14, c = -13, d = -12, e = -10,
  f = 4, g = 13, h = -1, i = -3, j = -11, k = -9,
  l = 16, m = -8, n = 11, o = 0, p = -6, q = -4,
  r = 15, s = 2, t = 9, u = -15, v = 14, w = -7,
  x = 7, y = -2, z = -5.

  max_{#1\in [1,26]}{|x_{#1}|} minimized to 16
  """

  Also, see the discussion of the Z model:
  http://www.comp.rgu.ac.uk/staff/ha/ZCSP/additional_problems/calculs_enfer/calculs_enfer.ps
  (which shows the same solution).

  This model
  a_max:16
  aa:[-13, -2, -14, -16, -10, 4, 13, -1, -3, -12, -11, 16, -9, 11, 0, -8, -6, 15, 9, 2, -15, 14, -7, 7, -4, -5]
  11.144555 seconds (29.14 M allocations: 5.421 GiB, 15.63% gc time)

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

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

function calculs_d_enfer(print_solutions=true,all_solutions=false)

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

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

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

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

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

                                                            "time_limit"=>13,

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

    nn = 26
    @variable(model, -100 <= aa[1:nn] <= 100, Int)
    a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,t,s,u,v,w,x,y,z = aa

    @variable(model, 0 <= a_abs[1:nn] <= 100, Int)
    for ii in 1:nn
        my_abs(model,aa[ii],a_abs[ii])
    end

    @variable(model, 0 <= a_max <= nn,Int)
    @constraint(model, a_abs .<= a_max)

    @constraint(model, aa in CS.AllDifferent())

    @constraint(model, z+e+r+o     == 0)
    @constraint(model, o+n+e       == 1)
    @constraint(model, t+w+o       == 2)
    @constraint(model, t+h+r+e+e   == 3)
    @constraint(model, f+o+u+r     == 4)
    @constraint(model, f+i+v+e     == 5)
    @constraint(model, s+i+x       == 6)
    @constraint(model, s+e+v+e+n   == 7)
    @constraint(model, e+i+g+h+t   == 8)
    @constraint(model, n+i+n+e     == 9)
    @constraint(model, t+e+n       == 10)
    @constraint(model, e+l+e+v+e+n == 11)
    @constraint(model, t+w+e+l+f   == 12)

    @objective(model,Min,a_max)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                aax = convert.(Integer,JuMP.value.(aa; result=sol))
                a_maxx = convert.(Integer,JuMP.value.(a_max; result=sol))
                println("a_max:$a_maxx")
                println("aa:$aax")

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

    return status
end

@time calculs_d_enfer()

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using AllocationOpt
using Documenter

DocMeta.setdocmeta!(AllocationOpt, :DocTestSetup, :(using AllocationOpt); recursive=true)

makedocs(;
    modules=[AllocationOpt],
    authors="The Graph Foundation",
    repo="https://github.com/graphprotocol/AllocationOpt.jl/blob/{commit}{path}#{line}",
    sitename="AllocationOpt.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://graphprotocol.github.io/AllocationOpt.jl",
        assets=String[],
    ),
    pages=["Home" => "index.md"],
)

deploydocs(; repo="github.com/graphprotocol/AllocationOpt.jl")


text:
using JWebImageDemo
using FileIO
using Test

const TEST_IMAGE = joinpath(@__DIR__, "resources/green_python.jpg")
const NUM = 5

const TMP_PATH = "tmp"
isdir(TMP_PATH) || mkdir(TMP_PATH)

@testset "image permute test" begin
    img = FileIO.load(TEST_IMAGE)
    for i = 2:NUM
        @time img_out = JWebImageDemo.permute_image(img, i)
        # save(Stream{format"PNG"}(io), data...)
        @test size(img) != size(img_out) ||
            any(((i, o),) -> i != o, zip(iterate(img), iterate(img_out)))
        save(joinpath(TMP_PATH, "out_$i.jpg"), img_out)
    end
end