text:
include("obstacles.jl")

mutable struct PlanningMap
    x_range::Vector{Float64}
    y_range::Vector{Float64}
    obstacle_grid::Matrix{Vector{Any}}
    all_obstacles::Vector{Any}
    params::parameters

    function PlanningMap(params::parameters)
        empty_obstacle_grid = Matrix{Vector{Any}}(undef, params.obstacle_grid_n, params.obstacle_grid_n)
        for i = 1:params.obstacle_grid_n, j = 1:params.obstacle_grid_n
            empty_obstacle_grid[i, j] = []
        end
        new(params.x_range, params.y_range, empty_obstacle_grid, [], params)
    end
end

function add_obstacle(this::PlanningMap, obst::obstacle)
    if !(obst in this.all_obstacles)
        push!(this.all_obstacles, obst)
    end
    for ids in obst.occupied_cells
        if !(obst in this.obstacle_grid[ids[1], ids[2]])
            push!(this.obstacle_grid[ids[1], ids[2]], obst)
        end
    end
end

function load_default_map(this::PlanningMap, params::parameters)
    obstacles = Vector{obstacle}()
    line_obst1 = line_obstacle(1, [5.1, 0.0], [5.15, 5.0], params)
    push!(obstacles, line_obst1)
    line_obst2 = line_obstacle(2, [3.0, 5.0], [9.0, 5.0], params)
    push!(obstacles, line_obst2)
    poly1 = convex_polygon(3, [[3.0, 8.0], [5.0, 9.0], [8.0, 7.5], [4.0, 6.5]], params)
    push!(obstacles, poly1)
    poly2 = convex_polygon(4, [[2.0, 1.5], [2.1, 2.6], [4.0, 3.0], [2.5, 2.0], [4.0, 1.0]], params)
    push!(obstacles, poly2)
    for obst in obstacles
        add_obstacle(this, obst)
    end
end

function visualize_map(this::PlanningMap)
    for obst in this.all_obstacles
        plot_obstacle(obst)
    end
end

function edge_collision(this::PlanningMap, p1, p2)
    obstacles_checked = []
    # Intersected cells of obstacle grid
    cells = intersected_cells([p1, p2], this.params)
    for cell_id in cells
        for obst in this.obstacle_grid[cell_id[1], cell_id[2]]
            if !(obst.id in obstacles_checked)
                push!(obstacles_checked, obst.id)
                if collision_check(p1, p2, obst)
                    return true
                end
            end
        end
    end
    return false
end


text:
using Test

@testset "test folder/two" begin
    @test true
    @test true
    @test true
end

ext: jl


text:
#
#      Graphmanifold – an array of points of one manifold,
#          where their vicinity or neighborhoods are given by a graph
#
# Manopt.jl, R. Bergmann, 2018-06-26
import Base: exp, log, show
import SparseArrays: issparse, findnz, sparse
export Graph, GraphVertexPoint, GraphVertexTVector, GraphEdgePoint, GraphEdgeTVector
export distance, dot, exp, log, manifoldDimension, norm, parallelTransport
export zeroTVector
export show, getValue
@doc doc"""
    Graph{M<:Manifold} <: Manifold

The graph manifold models manifold-valued data on a graph
$\mathcal G = (\mathcal V, \mathcal E)$, both on vertices
and edges as well as their interplay. The adjacency is stored in a matrix,
and may contain also the weights.

Since there are two possibilities in dimensions, $\lvert\mathcal V\rvert$ and
$\lvert\mathcal E\rvert$, the manifold itself will refer to the first one, while
depending on the type of `MPoint` one of them is returned.

# Fields
the default values are given in brackets
* `adjacency` – the (sparse) adjacency matrix, might also carry weights, i.e. all
  $a_{ij}>0$ refer to adjacent nodes $i$ and $j$
* `name` – (`A Graph manifold of \$Submanifold.`) name of the manifold
* `manifold` – the internal manifold present at vertices (edges) for
[`GraphVertexPoint`](@ref) ([`GraphEdgePoint`](@ref))
* `dimension` – stores the dimension of the manifold of a `GraphVertexPoint`
* `isDirected` – (`false`) indicates whether the graph is directed or not.
"""
struct Graph{mT<:Manifold} <: Manifold
  name::String
  manifold::mT
  dimension::Int
  adjacency::Mat where {Mat <: AbstractMatrix}
  isDirected::Bool
  abbreviation::String
  Graph{mT}(M::mT, adjacency::Mat, isDir::Bool=false) where {mT <: Manifold, Mat <: AbstractMatrix} = new(
    string("A Graph Manifold of ",M.name,"."),
    M,size(adjacency,1)*manifoldDimension(M),
    adjacency, isDir,
    string("GraphVertex(",M.abbreviation,",",repr(size(adjacency,1)),")") )
end
Graph(M::mT,A::Mat,dir::Bool=false) where {mT <: Manifold, Mat <: AbstractMatrix} = Graph{mT}(M,A,dir)

@doc doc"""
    GraphVertexPoint <: MPoint

A point graph vertex power manifold
$\mathcal M = \mathcal N^{\lvert\mathcal V\rvert}$
represented by a vector of corresponding [`MPoint`](@ref)s.
"""
struct GraphVertexPoint <: MPoint
  value::Vector{P} where {P<:MPoint}
  GraphVertexPoint(v::Vector{P} where {P<:MPoint}) = new(v)
end
getValue(x::GraphVertexPoint) = x.value;

@doc doc"""
    GraphVertexTVector <: TVector

A tangent vector $\xi\in T_x\mathcal M$ to the graph vertex power manifold
$\mathcal M = \mathcal N^{\lvert\mathcal V\rvert}$
represented by a vector of corresponding [`TVector`](@ref)s.
"""
struct GraphVertexTVector <: TVector
  value::Vector{T} where {T <: TVector}
  GraphVertexTVector(value::Vector{T} where {T <: TVector}) = new(value)
end
getValue(ξ::GraphVertexTVector) = ξ.value

@doc doc"""
    GraphEdgePoint <: MPoint

A point graph edge power manifold
$\mathcal M = \mathcal N^{\lvert\mathcal E\rvert}$
represented by a (sparse/not completely filled matrix of corresponding
[`MPoint`](@ref)s.
"""
struct GraphEdgePoint <: MPoint
  value::V where {V <: AbstractMatrix{P} where {P<:MPoint} }
  GraphEdgePoint(v::V where { V <: AbstractMatrix{P} where {P<:MPoint} }) = new(v)
end
getValue(x::GraphEdgePoint) = x.value;

@doc doc"""
    GraphEdgeTVector <: TVector

A tangent vector $\xi\in T_x\mathcal M$ to the graph edge power manifold
$\mathcal M = \mathcal N^{\lvert\mathcal E\rvert}$
represented by a (sparse/not completely filled) matrix of corresponding
[`TVector`](@ref)s.
"""
struct GraphEdgeTVector <: TVector
  value::V where {V <: AbstractMatrix{T} where {T<:TVector}}
  GraphEdgeTVector(v::V where { V <: AbstractMatrix{T} where {T<:TVector} }) = new(v)
end
getValue(ξ::GraphEdgeTVector) = ξ.value

# Functions
# ---
"""
    distance(M,x,y)

compute a vectorized version of distance on the [`Graph`](@ref) manifold `M`
for two [`GraphVertexPoint`](@ref) `x` and `y`.
"""
distance(M::Graph, x::GraphVertexPoint, y::GraphVertexPoint) = sqrt(sum( distance.(Ref(M.manifold), getValue(x), getValue(y) ).^2 ))
"""
    distance(M,x,y)

compute a vectorized version of distance on the [`Graph`](@ref) manifold `M`
for two [`GraphEdgePoint`](@ref) `x` and `y`.
"""
distance(M::Graph, x::GraphEdgePoint, y::GraphEdgePoint) = sqrt(sum( distance.(Ref(M.manifold), getValue(x), getValue(y) ).^2 ))
"""
    dot(M,x,ξ,ν)

computes the inner product as sum of the component inner products on the
[`Graph`](@ref) vertices.
"""
dot(M::Graph, x::GraphVertexPoint, ξ::GraphVertexTVector, ν::GraphVertexTVector) = sum(dot.(Ref(M.manifold),getValue(x), getValue(ξ), getValue(ν) ))
"""
    dot(M,x,ξ,ν)

computes the inner product as sum of the component inner products on the
[`Graph`](@ref) edges.
"""
dot(M::Graph, x::GraphEdgePoint, ξ::GraphEdgeTVector, ν::GraphEdgeTVector) = sum(dot.(Ref(M.manifold),getValue(x), getValue(ξ), getValue(ν) ))
"""
    exp(M,x,ξ[, t=1.0])

computes the product exponential map on the [`Graph`](@ref) vertices and
returns the corresponding [`GraphVertexPoint`](@ref).
"""
exp(M::Graph, x::GraphVertexPoint, ξ::GraphVertexTVector, t::Float64=1.0) = GraphVertexPoint( exp.(Ref(M.manifold), getValue(x) , getValue(ξ) ))
"""
    log(M,x,y)

computes the product logarithmic map on the [`Graph`](@ref) for two
[`GraphVertexPoint`](@ref) `x` and `y` and returns the
corresponding [`GraphVertexTVector`](@ref).
"""
log(M::Graph, x::GraphVertexPoint, y::GraphVertexPoint)::GraphVertexTVector = GraphVertexTVector(log.(Ref(M.manifold), getValue(x), getValue(y) ))
"""
    log(M,x,y)

computes the product logarithmic map on the [`Graph`](@ref) for two
[`GraphEdgePoint`](@ref) `x` and `y` and returns the
corresponding [`GraphEdgeTVector`](@ref).
"""
log(M::Graph, x::GraphEdgePoint, y::GraphEdgePoint)::GraphEdgeTVector = GraphEdgeTVector(log.(Ref(M.manifold), getValue(x), getValue(y) ))
"""
    manifoldDimension(x)

returns the (product of) dimension(s) of the [`Graph`](@ref) manifold `M` the
[`GraphVertexPoint`](@ref) `x` belongs to.
"""
manifoldDimension(x::GraphVertexPoint) = sum(manifoldDimension.( getValue(x) ) )
"""
    manifoldDimension(x)

returns the (product of) dimension(s) of the [`Graph`](@ref) manifold `M` the
[`GraphEdgePoint`](@ref)`x` belongs to.
"""
manifoldDimension(x::GraphEdgePoint) = sum(manifoldDimension.( getValue(x) ) )
"""
    manifoldDimension(M)

returns the (product of) dimension(s) of the [`Graph`](@ref) manifold `M` seen
as a vertex power manifold.
"""
manifoldDimension(M::Graph) = size(M.adjacency,1) * manifoldDimension(M.manifold)
"""
    norm(M,x,ξ)

norm of the [`GraphVertexTVector`](@ref) `ξ` induced by the metric on the
manifold components of the [`Graph`](@ref) manifold `M`.
"""
norm(M::Graph, x::GraphVertexPoint, ξ::GraphVertexTVector) = sqrt(sum(dot.(Ref(M.manifold),getValue(x),getValue(ξ),getValue(ξ))))
"""
    norm(M,x,ξ)

norm of the [`GraphEdgeTVector`](@ref) `ξ` induced by the metric on the manifold
components of the [`Graph`](@ref) manifold `M`.
"""
norm(M::Graph, x::GraphEdgePoint, ξ::GraphEdgeTVector) = sqrt(sum(dot.(Ref(M.manifold),getValue(x), getValue(ξ), getValue(ξ))))

@doc doc"""
    parallelTransport(M,x,ξ)

compute the product parallelTransport map on the [`Graph`](@ref) vertex power
manifold $\mathcal M^{\lvert\mathcal V\rvert}$ and returns
the corresponding [`GraphVertexTVector`](@ref).
"""
parallelTransport(M::Graph, x::GraphVertexPoint, y::GraphVertexPoint, ξ::GraphVertexTVector) = GraphVertexTVector( parallelTransport.(Ref(M.manifold), getValue(x), getValue(y), getValue(ξ)) )
"""
    parallelTransport(M,x,ξ)

compute the product parallelTransport map on the [`Graph`](@ref) edge power
manifold and returns the corresponding [`GraphVertexTVector`](@ref).
"""
parallelTransport(M::Graph, x::GraphEdgePoint, y::GraphEdgePoint, ξ::GraphEdgeTVector) = GraphVertexTVector( parallelTransport.(Ref(M.manifold), getValue(x), getValue(y), getValue(ξ)) )
"""
    randomMPoint(M,[,:Vertex])

compute a random point on the [`Graph`](@ref) manifold, where by default a point
on the vertices is produces, use `:Edge` to generate a `GraphEdgePoint`.
Further optional parameters are passed on to the element wise random point function.
"""
randomMPoint(M::Graph,::Val{:Vertex},options...) = GraphVertexPoint(
    [ randomMPoint( M.manifold, options...) for i=1:size(M.adjacency,1) ]
)
randomMPoint(M::Graph,::Val{:Edge},options...) =  GraphEdgePoint(
    [ randomMPoint( M.manifold, options...) for i=1:sum( M.adjacency .> 0) ]
)

"""
    randomTVector(M,x)

compute a random [`GraphEdgeTVector`](@ref) to the [`GraphVertexPoint`](@ref) `x`
on the [`Graph`](@ref) manifold `M` by calling the inner random vector generation
for every edge
"""
randomTVector(M::Graph, x::GraphVertexPoint,options...) = GraphVertexTVector(
    [
        (M.adjacency[i,j] > 0) ? M.randomTVector(M.manifold, getValue(x)[i],options...) : zeroTVector(M,getValue(x)[i])
        for i = 1:size(M.adjacency,1) for j=1:size(M.adjacency,2)
    ]
)

typicalDistance(M::Graph) = sqrt( size(M.adjacency,1) ) * typicalDistance(M.manifold);
@doc doc"""
    ξ = zeroTVector(M,x)

returns a zero vector in the tangent space $T_x\mathcal M$ of the
[`GraphVertexPoint`](@ref) $x\in\mathcal M$ on the [`Graph`](@ref) vertex manifold `M`.
"""
zeroTVector(M::Graph, x::GraphVertexPoint) = GraphVertexTVector( zeroTVector.(Ref(M.manifold), getValue(x) )  )
@doc doc"""
    ξ = zeroTVector(M,x)

returns a zero vector in the tangent space $T_x\mathcal M$ of the
[`GraphEdgePoint`](@ref) $x\in\mathcal M$ on the [`Graph`](@ref) edge manifold `M`.
"""
zeroTVector(M::Graph, x::GraphEdgePoint) = GraphEdgeTVector( zeroTVector.(Ref(M.manifold), getValue(x) )  )
"""
    startEdgePoint(M,x)

For a [`Graph`](@ref) manifold and a [`GraphVertexPoint`](@ref), this
function constructs the corresponding [`GraphEdgePoint`](@ref), such that each
edge has its start point vertex value assigned.
"""
function startEdgePoint(M::Graph, x::GraphVertexPoint)::GraphEdgePoint
  if length(getValue(x)) == 0
    throw( ErrorException("vertexToStartEdgePoint::No node given"))
  end
  if issparse(M.adjacency)
    (s,e,v) = findnz(M.adjacency)
    return GraphEdgePoint( sparse(s,e, getValue(x)[ s ]) )
  else
    sA = size(M.adjacency)
    A = Matrix(typeof(x[1]),sA)
    for i=1:sA[1]
      for j=1:sA[2]
        if M.adjacency[i,j] > 0
          A[i,j] = getValue(x)[i]
        else
          A[i,j] = missing;
        end
      end
    end
  end
end
"""
    sumEdgeTVector(M,ξ)

return the [`GraphVertexTVector`](@ref) where edge tangents are summed in their
start point.

For an [`GraphEdgeTVector`](@ref) `ξ` on a [`Graph`](@ref) manifold `M`
this function assumes that all edge tangents are attached in a tangent space
corresponding to the same point on the base manifold, i.e. all these vectors can
be summed. This sum per vectex is then returned as a
[`GraphVertexTVector`](@ref).
"""
function sumEdgeTVector(M::Graph, ξ::GraphEdgeTVector, weighted::Bool=false)
  (s,e,v) = findnz(M.adjacency)
  init = falses(size(M.adjacency,1)) # vector initialized?
  lξ = getValue(ξ)
  ν = Vector(  eltype( lξ ),size(M.adjacency,1)  )
  for i=1:length(s) # all edges
    if !init[ s[i] ]
      ν[ s[i] ] = lξ[s[i],e[i]] # initialize (i.e. also carry base if it exists)
      init[ s[i] ] = true
    else
      ν[ s[i] ] = ν[ s[i] ] + lξ[s[i], e[i]]
    end
  end
  return GraphVertexTVector(ν)
end
#
#
# Display functions for the structs
show(io::IO, M::Graph) = print(io,string("The manifold on vertices and edges of a graph of ",repr(M.manifold), " of (vertex manifold) dimension ",repr(M.dimension),".") );
show(io::IO, p::GraphVertexPoint) = print(io,string("GraphVertexV[",join(repr.(p.value),", "),"]"));
show(io::IO, ξ::GraphVertexTVector) = print(io,string("GraphVertexT[",join(repr.(ξ.value),", "),"]"));
show(io::IO, p::GraphEdgePoint) = print(io,string("GraphEdgeV[",join(repr.(p.value),", "),"]"));
show(io::IO, ξ::GraphEdgeTVector) = print(io,string("GraphEdgeT[", join(repr.(ξ.value),", "),"]"));


text:
module words_extractor_jl

using Distributed
using YAML
using Glob

const outdir = "words"

function worker(yaml_path, i, count)
    path = get_filepath(yaml_path)
    words = get_words(yaml_path)
    write(path, join(words, "\n"))
    println("[$(lpad(i, 3, ' '))/$count] $path")
end

function get_words(yaml_path)
    text_path = replace(yaml_path, ".yml" => ".txt")
    text = read(text_path, String) |> lowercase
    split(text, r"[\W\d]+") |> Set |> collect
end

function get_filepath(path)
    meta = YAML.load_file(path)
    """./$outdir/$(meta["lang"])-$(meta["code"]).txt"""
end

function rdir(dir::AbstractString, pat::Glob.FilenameMatch)
    result = String[]
    for (root, dirs, files) in walkdir(dir)
        filepaths = joinpath.(root, files)
        append!(result, filter!(f -> occursin(pat, f), filepaths))
    end
    return result
end

rdir(dir::AbstractString, pat::AbstractString) = rdir(dir, Glob.FilenameMatch(pat))

function main()
    if ispath(outdir)
        rm(outdir, recursive = true)
    end
    mkdir(outdir)
    paths = rdir("../data", fn"../data/??/*.yml")
    count = length(paths)
    i = 1
    Threads.@threads for path in paths
        # println("Spawn $path")
        worker(path, i, count)
        i += 1
    end
end

addprocs()
println(string("Workers ", nworkers()))
println(string("Processing... using ", Threads.nthreads(), " threads"))
@time main()
end # module


text:
export run_opf, run_ac_opf, run_dc_opf

""
function run_ac_opf(file, solver; kwargs...)
    return run_opf(file, ACPPowerModel, solver; kwargs...)
end

""
function run_dc_opf(file, solver; kwargs...)
    return run_opf(file, DCPPowerModel, solver; kwargs...)
end

""
function run_opf(file, model_constructor, solver; kwargs...)
    return run_generic_model(file, model_constructor, solver, post_opf; kwargs...)
end

""
function post_opf(pm::GenericPowerModel)
    variable_voltage(pm)
    variable_generation(pm)
    variable_branch_flow(pm)
    variable_dcline_flow(pm)

    objective_min_fuel_cost(pm)

    constraint_voltage(pm)

    for i in ids(pm, :ref_buses)
        constraint_theta_ref(pm, i)
    end

    for i in ids(pm, :bus)
        constraint_kcl_shunt(pm, i)
    end

    for i in ids(pm, :branch)
        constraint_ohms_yt_from(pm, i)
        constraint_ohms_yt_to(pm, i)

        constraint_voltage_angle_difference(pm, i)

        constraint_thermal_limit_from(pm, i)
        constraint_thermal_limit_to(pm, i)
    end

    for i in ids(pm, :dcline)
        constraint_dcline(pm, i)
    end
end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using Test
using UnitCommitment

push!(Base.LOAD_PATH, @__DIR__)
UnitCommitment._setup_logger()

const ENABLE_LARGE_TESTS = ("UCJL_LARGE_TESTS" in keys(ENV))

@testset "UnitCommitment" begin
    include("usage.jl")
    @testset "import" begin
        include("import/egret_test.jl")
    end
    @testset "instance" begin
        include("instance/read_test.jl")
    end
    @testset "model" begin
        include("model/formulations_test.jl")
    end
    @testset "XavQiuWanThi19" begin
        include("solution/methods/XavQiuWanThi19/filter_test.jl")
        include("solution/methods/XavQiuWanThi19/find_test.jl")
        include("solution/methods/XavQiuWanThi19/sensitivity_test.jl")
    end
    @testset "transform" begin
        include("transform/initcond_test.jl")
        include("transform/slice_test.jl")
        @testset "randomize" begin
            include("transform/randomize/XavQiuAhm2021_test.jl")
        end
    end
    @testset "validation" begin
        include("validation/repair_test.jl")
    end
end


text:
""" Diagrams of a given shape.
"""
module ShapeDiagrams
export Span, Cospan, left, right, apex, base

using ...Theories: dom, codom

""" Span of morphisms in a category.
"""
struct Span{Left,Right}
  left::Left
  right::Right

  function Span(left::Left, right::Right; strict::Bool=true) where {Left,Right}
    if strict && dom(left) != dom(right)
      error("Domains of legs in span do not match: $left vs $right")
    end
    new{Left,Right}(left, right)
  end
end

apex(span::Span) = dom(span.left) # == dom(span.right)
left(span::Span) = span.left
right(span::Span) = span.right

""" Cospan of morphisms in a category.
"""
struct Cospan{Left,Right}
  left::Left
  right::Right

  function Cospan(left::Left, right::Right; strict::Bool=true) where {Left,Right}
    if strict && codom(left) != codom(right)
      error("Codomains of legs in cospan do not match: $left vs $right")
    end
    new{Left,Right}(left, right)
  end
end

base(cospan::Cospan) = codom(cospan.left) # == codom(cospan.right)
left(cospan::Cospan) = cospan.left
right(cospan::Cospan) = cospan.right

end


text:
function getnormal(params, constants, targetdata)
  simdata = rand(Normal(params...), 100)
  ksdist(simdata, targetdata), 1
end

function getuniformdist(params, constants, targetdata)
  params = sort(params)
  simdata = rand(Uniform(params...), 100)
  ksdist(simdata, targetdata), 1
end

srand(1234)
targetdata = rand(Normal(3.5, 0.44), 100)

ABCsetup = ABCSMCModel([getnormal, getuniformdist, getnormal], [2, 2, 2], 0.1, [Prior([Uniform(0.0, 20.0), Uniform(0.0, 2.0)]), Prior([Uniform(0.0, 20.0), Uniform(0.0, 2.0)]), Prior([Uniform(0.0, 20.0), Uniform(0.0, 2.0)])]; nparticles = 100, maxiterations = 10^5)

#test model perturbation kernel
Niterations = 10^6
m = zeros(Int64, Niterations )
mstar = 1
modelprob = [1/3, 1/3, 1/3]
for i in 1:Niterations
  mdoublestar = ApproxBayes.perturbmodel(ABCsetup, mstar, modelprob)
  m[i] = mdoublestar
end

@test isapprox(ABCsetup.modelkern, sum(m.==1)/length(m), rtol = 0.01)

@test isapprox((1 - ABCsetup.modelkern)./(ABCsetup.nmodels - 1), sum(m.==2)/length(m), rtol = 0.01)

#test model perturbation kernel when one model has died out
Niterations = 10^6
m = zeros(Int64, Niterations )
mstar = 1
modelprob = [0.5, 0.0, 0.5]
for i in 1:Niterations
  mdoublestar = ApproxBayes.perturbmodel(ABCsetup, mstar, modelprob)
  m[i] = mdoublestar
end

@test isapprox(ABCsetup.modelkern, sum(m.==1)/length(m), rtol = 0.01)
@test isapprox(0.0, sum(m.==2)/length(m))
@test isapprox(1.0 - ABCsetup.modelkern, sum(m.==3)/length(m), rtol = 0.01)

#test get particle weights
ABCrejresults = runabc(ABCRejectionModel(
          map(x -> x.simfunc, ABCsetup.Models),
          map(x -> x.nparams, ABCsetup.Models),
          ABCsetup.Models[1].ϵ1,
          map(x -> x.prior, ABCsetup.Models);
          nparticles = ABCsetup.Models[1].nparticles,
          maxiterations = ABCsetup.Models[1].maxiterations),
          targetdata);

oldparticles, weightsA = ApproxBayes.setupSMCparticles(ABCrejresults, ABCsetup)
weightsA, modelprob = ApproxBayes.getparticleweights(oldparticles, ABCsetup)

#test if modelprob is the same as from ABCrejresults
@test modelprob[:] ≈ ABCrejresults.modelfreq

for i in 1:ABCsetup.nmodels
  @test (map(x -> x.model, oldparticles).==i) == (weightsA[i, :].> 0.0)
end

#test get proposal function
ABCsetup = ABCSMC(getnormal, 2, 0.1, Prior((Uniform(0.0, 20.0), Exponential(1.0))); nparticles = 100, maxiterations = 10^5)
Niterations = 10^6
p1 = zeros(Float64, Niterations)
p2 = zeros(Float64, Niterations)
@time for i in 1:Niterations
  p1[i], p2[i] = ApproxBayes.getproposal(ABCsetup.prior, 2)
end
@test isapprox(mean(p1), 10.0, rtol = 0.001)
@test isapprox(mean(p2), 1.0, rtol = 0.001)

#test get proposal function with models with different Priors
ABCsetup = ABCSMCModel([getnormal, getuniformdist], [2, 2], 1.0,
[Prior([Exponential(1), Uniform(0.0, 1.0)]), Prior([Uniform(0.0, 3.0), Normal(2.0, 0.1)])]; nparticles = 100, maxiterations = 10^5)
Niterations = 10^6
p1 = Float64[]
p2 = Float64[]
p3 = Float64[]
p4 = Float64[]
@time for i in 1:Niterations
  m = rand(1:2)
  if m == 1
    x = ApproxBayes.getproposal(ABCsetup.Models[m].prior, 2)
    push!(p1, x[1])
    push!(p2, x[2])
  else
    x = ApproxBayes.getproposal(ABCsetup.Models[m].prior, 2)
    push!(p3, x[1])
    push!(p4, x[2])
  end
end
@test isapprox(mean(p1), 1.0, rtol = 0.01)
@test isapprox(mean(p2), 0.5, rtol = 0.01)
@test isapprox(mean(p3), 1.5, rtol = 0.01)
@test isapprox(mean(p4), 2.0, rtol = 0.01)


text:
using DataFrames
using CSV
using Plots
using Statistics
using GLM

dfs = CSV.read("zombiecar2_sweep_res5_steptrue.csv");

for c in [:infected_avg, :infected_std, :infected_lo, :infected_hi,
            :infected_median, :infected_q05, :infected_q95 ]
    dfs[!, Symbol(string("pc_",c))] = dfs[!, c] ./ dfs.n_agents
end

# log(steps) = 2g_nodes*log(n_agents)/n_agents

dfs.e_sim_time = 2dfs.g_nodes.*log.(dfs.n_agents)./dfs.n_agents
dfs.log_n_x_n = log.(dfs.n_agents)./dfs.n_agents

function get_slice(dfs::DataFrame; n_agents=500, discretize_m=50.0, jump=0)
    dfs[ (dfs.n_agents.==n_agents) .& (dfs.discretize_m .== discretize_m) .& (dfs.jump .== jump), :]
end

function plot_res!(p, d, color, lab)
    Plots.plot!(p, d.step, d.pc_infected_avg; linestyle=:solid, lab=lab, color=color)
    Plots.plot!(p, d.step, d.pc_infected_q05; linestyle=:dash, lab="", color=color)
    Plots.plot!(p, d.step, d.pc_infected_q95; linestyle=:dash, lab="", color=color)
end

function plot_sim_res(dfs::DataFrame,
        kw_list::Vector{NamedTuple{(:n_agents, :discretize_m, :jump),Tuple{Int64,Float64,Int64}}})

    ds = [get_slice(dfs; kw...) for kw in kw_list]
    for i in 1:length(kw_list)
        @assert nrow(ds[i]) > 0 "ERRROR no rows found for $(string(kw_list[i]))"
    end
    labs = ["k=$(kw.n_agents), Jump-over=$(kw.jump==1 ? "Yes" : "No"), discretization=$(round(Int, kw.discretize_m))m" for kw in kw_list]
    ordering = sortperm(ds; lt = (x,y) -> x.step[end] < y.step[end] )
    ds = ds[ordering]
    labs = labs[ordering]

    ncolor = max(length(kw_list),2)
    cols = reshape( range(colorant"blue", stop=colorant"red",length=ncolor), 1, ncolor );
    p = Plots.plot(;xlabel = "Simulation steps", ylabel="Share of agents who received the message", ylim=(0.0,1.1),
        legend=:bottomright, size=(500,300))

    for i in 1:length(labs)
        plot!(p,[0],[0],color=cols[i], lab=labs[i])
    end
    for i in 1:length(kw_list)
        plot_res!(p, ds[i], cols[i], "")
    end
    p
end

Plots.pgfplotsx()

plt = plot_sim_res(dfs, [
    (n_agents=100, discretize_m=50.0, jump=1 ),
    (n_agents=250, discretize_m=50.0, jump=1 ),
    (n_agents=400, discretize_m=50.0, jump=1 ),
    (n_agents=700, discretize_m=50.0, jump=1 ),
    (n_agents=1000, discretize_m=50.0, jump=1 ),
])
savefig(plt, "../5d63f3b5c3791641ba87e19f/message_prop_med.pdf")

plt = plot_sim_res(dfs, [
    (n_agents=1000, discretize_m=50.0, jump=1 ),
    (n_agents=2000, discretize_m=50.0, jump=1 ),
    (n_agents=5000, discretize_m=50.0, jump=1 ),
    (n_agents=10000, discretize_m=50.0, jump=1 ),
    (n_agents=20000, discretize_m=50.0, jump=1 ),
])
savefig(plt, "../5d63f3b5c3791641ba87e19f/message_prop_high.pdf")


text:
using Base.Broadcast: AbstractArrayStyle, Broadcasted

struct LayeredArrayStyle{N} <: AbstractArrayStyle{N} end
(::Type{<: LayeredArrayStyle})(::Val{N}) where {N} = LayeredArrayStyle{N}()

Broadcast.BroadcastStyle(::Type{<: AbstractLayeredArray{<: Any, <: Any, N}}) where {N} = LayeredArrayStyle{N}()

Broadcast.broadcastable(x::AbstractLayeredArray) = x
Broadcast.broadcastable(x::Broadcasted{<: LayeredArrayStyle}) = copy(x)
Broadcast.instantiate(x::Broadcasted{<: LayeredArrayStyle}) = x

function Base.copy(bc::Broadcasted{<: LayeredArrayStyle})
    f, args = bc.f, bc.args
    layer, args′ = lazyables(f, args...)
    T = return_eltype(f, args′...)
    LazyLayeredArray{layer, T}(Broadcast.broadcasted(f, args′...))
end

function Base.copyto!(dest::AbstractLayeredArray, src::Broadcasted{<: LayeredArrayStyle})
    set!(dest, copy(src))
end


text:
module SimpleARules

include("arules.jl")

export apriori, frequent

end # module


text:
InputFolder = './Images/Vessels/';
OutputFolder = './Results/Images/Vessels/';
Fill = -1;

@iA = '*.tif';

@fxg_mRadStdThr [iA] > [M];
params.GRad = 1;
params.RSet = 5;
params.Fracstd = 0.75;
params.SegThr = 0.075;
/endf

@fxm_kSklMorph [M] > [C];
params.Mode = 'thin';
params.MinBrcLgth = 3;
/endf

/show iA > C;
/keep C > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
## Conversion

"""
Convert a domain object to a domain object with the given element type `T`.
"""
convert_domain(::Type{T}, d::Domain{T}) where {T} = d
convert_domain(::Type{T}, d::Domain{S}) where {S,T} = convert(Domain{T}, d)

convert_domain(::Type{T}, d) where {T} = _convert_domain(T, d, eltype(d))
_convert_domain(::Type{T}, d, ::Type{T}) where {T} = d
_convert_domain(::Type{T}, d, ::Type{S}) where {S,T} =
    error("Don't know how to convert the element type of $(typeof(d)) to $S. Please implement `convert_domain` for this combination of types.")

convert_domain(::Type{T}, d::Number) where {T} = convert(T, d)

_convert_domain(::Type{T}, d::AbstractArray{S}, ::Type{S}) where {S,T} = (b = similar(d, T); b .= d)
_convert_domain(::Type{T}, d::Set{S}, ::Type{S}) where {S,T} = convert(Set{T}, d)

## Promotion

"""
Promote the given domains to have a common element type.
"""
promote_domain() = ()
promote_domain(domains...) = promote_domains(domains)

"""
Promote an iterable list of domains to have a common element type.
"""
promote_domains(domains) = _promote_domains(mapreduce(eltype, promote_type, domains), domains)

# We simply rely on Julia's promotion system to promote the element types, and
# we invoke `convert_domain` on the domains with the resulting type.
_promote_domains(::Type{T}, domains) where {T} = convert_domain.(T, domains)


text:
module Client
using HTTP, JSON3, StructTypes, Dates, JLD2, Distributed, REPL.TerminalMenus, Reexport
using ..DFControl
using ..Utils
using ..FileIO
@reexport using ..Servers
@reexport using ..Structures
@reexport using ..Calculations
@reexport using ..Jobs

@inline function JSON3.read(::StructTypes.Mutable, buf, pos, len, b, ::Type{Calculation};
                            kw...)
    x = Calculation{Calculations.NoPackage}(""; exec = Exec())
    pos, x = JSON3.read!(StructTypes.Mutable(), buf, pos, len, b, Calculation, x; kw...)
    if x.exec.exec ∈ Calculations.WAN_EXECS
        p = Wannier90
    elseif x.exec.exec ∈ Calculations.QE_EXECS
        p = QE
    elseif x.exec.exec ∈ Calculations.ELK_EXECS
        p = ELK
    else
        @warn "Package not identified from execs $(x.exec)."
    end

    t = Calculation{p}(x.name, x.flags, x.data, x.exec, x.run,  x.infile,
                       x.outfile)
    return pos, t
end

using ..Calculations: set_name!, set_kpoints!, data; export set_name!, set_kpoints!, data
using ..Structures: set_pseudos!, element; export set_pseudos!, element

include("job.jl")
export submit, save, isrunning, state, versions, last_version, switch_version!, rm_version!, outputdata,
       registered_jobs, running_jobs,
       environment_from_jobscript, get_environment, add_environment, rm_environment!,
       known_execs, get_exec

include("pseudos.jl")
export configure_pseudoset, rm_pseudoset!

end


text:
@testset "basics" begin
header_ref = BGEN.Header(example_10bits)

@testset "header" begin
    header_test = BGEN.Header(example_10bits)
    @test header_test.offset == 0x0000178c # 6028
    @test header_test.header_length == 0x00000014 # 20
    @test header_test.n_variants == 0x000000c7 # 199
    @test header_test.n_samples == 0x000001f4 # 500
    @test header_test.compression == 1
    @test header_test.layout == 2
    @test header_test.has_sample_ids == 1
    end

    @testset "samples_separate" begin
    n_samples = 500
    samples_test2 = BGEN.get_samples(example_sample, n_samples)
    samples_correct = [(@sprintf "sample_%03d" i) for i in 1:n_samples]
    @test all(samples_correct .== samples_test2)
    samples_test3 = BGEN.get_samples(n_samples)
    @test all([string(i) for i in 1:n_samples] .== samples_test3)
end

bgen = BGEN.Bgen(example_10bits)
@testset "bgen" begin
    @test bgen.fsize == 223646
    @test bgen.header == header_ref
    n_samples = bgen.header.n_samples
    samples_correct = [(@sprintf "sample_%03d" i) for i in 1:n_samples]
    @test all(samples_correct .== bgen.samples)
    variants = parse_variants(bgen)
    var = variants[4]
    @test length(variants) == 199
    @test var.offset == 0x0000000000002488
    @test var.geno_offset == 0x00000000000024b1
    @test var.next_var_offset == 0x0000000000002902
    @test var.geno_block_size == 0x00000451
    @test var.n_samples == 0x000001f4
    @test var.varid == "SNPID_5"
    @test var.rsid == "RSID_5"
    @test var.chrom == "01"
    @test var.pos == 0x00001388
    @test var.n_alleles == 2
    @test all(var.alleles.== ["A", "G"])
end

@testset "preamble" begin
    io, v, h = bgen.io, parse_variants(bgen)[1], bgen.header
    decompressed = BGEN.decompress(io, v, h)
    preamble = BGEN.parse_preamble(decompressed, h, v)
    @test preamble.phased == 0
    @test preamble.min_ploidy == 2
    @test preamble.max_ploidy == 2
    @test all(preamble.ploidy .== 2)
    @test preamble.bit_depth == 10
    @test preamble.max_probs == 3
    @test length(preamble.missings) == 1
    @test preamble.missings[1] == 1
end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule GAP_pkg_digraphs_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("GAP_pkg_digraphs")
JLLWrappers.@generate_main_file("GAP_pkg_digraphs", UUID("ce5c30f6-a98c-53bc-b643-089f5717f8db"))
end  # module GAP_pkg_digraphs_jll


text:
using Documenter, BitInformation

makedocs(
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename="BitInformation.jl",
    authors="M Klöwer",
    modules=[BitInformation],
    pages=["Home"=>"index.md",
            "Bitwise information"=>"bitinformation.md",
            "Transformations"=>"transformations.md",
            "Rounding"=>"rounding.md",
            "Function index"=>"functions.md"]
)

deploydocs(
    repo = "github.com/milankl/BitInformation.jl.git",
    devbranch = "main"
)


text:
module Rqa

default_min_length = 2

"""Helper struct keeping sizes of diagonal lines and total area
of recurrence plot. In case of regular recurrence plot only upper triangle
(without the diagonal) is considered; in case of cross recurrence plot it is
the whole rectanle"""
type RQAdata
  area :: Int64
  lines :: Dict{Int64, Int64}
end

"Computes recurrence rate - probability of recurrence"
function RR(data :: RQAdata)
  if isempty(data.lines)
    0
  else
    sum(kv -> kv[1]*kv[2], data.lines) / data.area
  end
end

"Computes ratio of points lying on lines; by line we understand long enough
sequence of recurrences"
function DET(data :: RQAdata, min_length :: Int = default_min_length)
  if isempty(data.lines)
    NaN
  else
    longblack = sum(data.lines) do kv
      if kv[1] >= min_length
        kv[1]*kv[2]
      else
        0
      end
    end

    allblack = sum(kv -> kv[1]*kv[2], data.lines)
    longblack / allblack
  end
end

"Computes mean length of line; by line we understand long enough sequence
of recurrences"
function LMEAN(data :: RQAdata, min_length :: Int = default_min_length)
  llines = filter((k,v) -> k >= min_length, data.lines)
  if isempty(llines)
    NaN
  else
    length_sum = sum(kv -> kv[1] * kv[2], llines)
    count = sum(kv -> kv[2], llines)
    length_sum / count
  end
end

"Computes entropy of lines; by line we understand long enough sequence
of recurrences"
function LENT(data :: RQAdata, min_length :: Int = default_min_length)
  llines = filter((k,v) -> k >= min_length, data.lines)
  if isempty(llines)
    NaN
  else
    cnt = sum(kv -> kv[2], llines)
    sum(kv -> kv[2] / cnt * log2(cnt / kv[2]), llines)
  end
end

"""
Calculates the lines count for RQA. Takes into consideration only one triangle
of the plot, without the diagonal"""
function rqadata(sig :: AbstractVector, eq :: Function = (x, y) -> x == y)
  n = length(sig)
  lines = Dict{Int64, Int64}()

  for lag = 1:(n - 1)
    l = 0
    for i in (lag + 1):n
      if eq(sig[i], sig[i - lag])
        l += 1
      else
        if l > 0
          lines[l] = get(lines, l, 0) + 1
        end
        l = 0
      end
    end
    if l > 0
      lines[l] = get(lines, l, 0) + 1
    end
  end

  RQAdata(div(n*(n - 1), 2), lines)
end

"Computes the recurrence data for Cross Recurrence Quantification Analysis"

function crqadata(
    sig1 :: AbstractVector,
    sig2 :: AbstractVector,
    eq :: Function = (x, y) -> x == y)
  n = length(sig1)
  m = length(sig2)
  lines = Dict{Int64, Int64}()

  # iterates over all diagonals in the rectangle
  # lag 0 is the "main" diagonal - one passing through the top left corner
  for lag = (1 - n):(m - 1)
    l = 0
    for i in max(1, 1-lag):min(n, m - lag)
      if eq(sig1[i], sig2[i + lag])
        l += 1
      else
        if l > 0
          lines[l] = get(lines, l, 0) + 1
        end
        l = 0
      end
    end
    if l > 0
      lines[l] = get(lines, l, 0) + 1
    end
  end

  RQAdata(n*m, lines)
end

end


text:
#------------------------------------------------------------------------------
# Author: Edward Tate <edward.tate@erlang-solutions.com>
#------------------------------------------------------------------------------

export read_erlang_term, read_binary_stream
export write_erlang_term, write_binary_stream

#------------------------------------------------------------------------------
# Byte Reading Functions
#------------------------------------------------------------------------------

function read_erlang_term (s::TCPSocket)
    bytes = Array(Uint8, 0)
    length_header = readbytes(s, 4)
    nbytes = decode_integer(length_header)
    nbytes_read = 0
    while nbytes_read < nbytes
        next_bytes = readbytes(s, nbytes)
        nbytes_read += length(next_bytes)
        bytes = vcat(next_bytes, bytes)
    end
    vcat(length_header, bytes)
end

function read_binary_stream (s::TCPSocket)
    bytes = Array(Uint8, 0)
    length_header = readbytes(s, 4)
    nbytes = decode_integer(length_header)
    nbytes_read_total = 0
    while nbytes > 0
        nbytes_read = 0
        while nbytes_read < nbytes
            next_bytes = readbytes(s, nbytes)
            nbytes_read += length(next_bytes)
            bytes = vcat(next_bytes, bytes)
        end
        nbytes_read_total += nbytes_read
        length_header = readbytes(s, 4)
        nbytes = decode_integer(length_header)
    end
    bytes, nbytes_read_total
end

#------------------------------------------------------------------------------
# Byte Writing Functions
#------------------------------------------------------------------------------

function write_erlang_term (s::TCPSocket, bytes::Array{Uint8})
    nbytes = length(bytes)
    nbytes_written = 0
    nbytes_written_total = 0
    while nbytes_written < nbytes
        nbytes_written = write(s, bytes)
        bytes = bytes[nbytes_written:end]
        nbytes_written_total += nbytes_written
    end
    flush(s)
    nbytes_written_total
end

function write_binary_stream (s::TCPSocket, bytes::Array{Uint8})
    length_header = encode_integer(uint32(length(bytes)))
    bytes = vcat(length_header, bytes)
    bytes = vcat(bytes, uint8([0,0,0,0]))
    nbytes = length(bytes)
    nbytes_written = 0
    while nbytes_written < nbytes
        nbytes_written = write(s, bytes)
        bytes = bytes[nbytes_written:end]
        nbytes_written_total += nbytes_written
    end
    flush(s)
    nbytes_written_total
end

max_issues_repo_name: esl/FsServer

max_forks_repo_name: esl/FsServer


text:
# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors    #src
# This Source Code Form is subject to the terms of the Mozilla Public License   #src
# v.2.0. If a copy of the MPL was not distributed with this file, You can       #src
# obtain one at https://mozilla.org/MPL/2.0/.                                   #src

# # The correlation problem

# Given three random variables A, B, C and given bounds on two of the three
# correlation coefficients:
#
#     -0.2 <= ρ_AB <= -0.1
#     0.4 <= ρ_BC <=  0.5
#
# We can use the following property of the correlations to determine bounds on
# ρ_AC by solving a SDP:
#
#     |  1    ρ_AB  ρ_AC |
#     | ρ_AB   1    ρ_BC |  ≽ 0
#     | ρ_AC  ρ_BC   1   |

using JuMP
import SCS
import Test  #src

function example_corr_sdp()
    model = Model(SCS.Optimizer)
    set_silent(model)
    @variable(model, X[1:3, 1:3], PSD)
    ## Diagonal is 1s
    @constraint(model, X[1, 1] == 1)
    @constraint(model, X[2, 2] == 1)
    @constraint(model, X[3, 3] == 1)
    ## Bounds on the known correlations
    @constraint(model, X[1, 2] >= -0.2)
    @constraint(model, X[1, 2] <= -0.1)
    @constraint(model, X[2, 3] >=  0.4)
    @constraint(model, X[2, 3] <=  0.5)
    ## Find upper bound
    @objective(model, Max, X[1, 3])
    optimize!(model)
    println("An upper bound for X[1, 3] is $(value(X[1, 3]))")
    Test.@test value(X[1, 3]) ≈ 0.87195 atol = 1e-4  #src
    ## Find lower bound
    @objective(model, Min, X[1, 3])
    optimize!(model)
    println("A lower bound for X[1, 3] is $(value(X[1, 3]))")
    Test.@test value(X[1, 3]) ≈ -0.978 atol = 1e-3  #src
    return
end

example_corr_sdp()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using DataStructures: Deque

#=------------------------------------------------------------------------------
This section contains the `Tile` data structure, `TileCircle` data structure,
and functions that operate on these data structures. The `TileCircle` provides
and interface to an underlying double-ended queue and represents the circle
of tiles placed during a game. Tiles are placed in skip-wise fashion, that is,
each tile is placed between the tiles 1 and 2 spaces clockwise from the last
placed tile.
------------------------------------------------------------------------------=#

"""
A `Tile` is a data structure containing the value of the tile and
the ID of the player who played it.
"""
struct Tile
    number::Int
    played_by::Int
end

"""
TileCircle only contains a double-ended queue, a data structure with efficent
push/pop operations on both ends. I'm wrapping the deque in a data structure
so I can specialize functions for it.
"""
struct TileCircle
    inner::Deque{Tile}
end
TileCircle() = TileCircle(Deque{Tile}()) # Default empty constructor

"""
    rotate!(circle::TileCircle)

Rotate the entire circle counter-clockwise by one 'tick', which has the effect
of moving the insertion point for new values forward by one.
"""
function rotate!(circle::TileCircle)
    isempty(circle.inner) && return
    push!(circle.inner, popfirst!(circle.inner))
end

"""
    push!(circle::TileCircle, tile::Tile)

Specialization of the base `push!()` function for a `TileCircle`, always inserts
new values after skipping the next clockwise value.
"""
function Base.push!(circle::TileCircle, tile::Tile)
    rotate!(circle)
    push!(circle.inner, tile)
end

"Dispatches for methods on the inner double-ended queue."
Base.isempty((; inner)::TileCircle) = isempty(inner)
Base.length((; inner)::TileCircle)  = length(inner)

"""
    next_tile(circle::TileCircle)

Returns the next tile in sequence from the current position.
"""
function next_tile(circle::TileCircle)
    isempty(circle) && error("This `TileCircle` is empty!")
    return first(circle.inner)
end

"""
    rotate_to!(circle::TileCircle)

Rotate the circle until the next tile is the tile whose value matches `n`.
"""
function rotate_to!(circle::TileCircle, n::Int=0)
    0 <= n < length(circle) || error("There's no `$n` tile in this `TileCircle`!")
    while next_tile(circle).number != n
        rotate!(circle)
    end
end

"""
    collect(circle::TileCircle)

Collect the contents of the circle into a Vector, with tiles ordered by
position (i.e., tile `0` in the first position, tile `1` in the second
position, etc.).
"""
function Base.collect(circle::TileCircle)
    new_copy = deepcopy(circle)
    isempty(circle) || rotate_to!(new_copy)
    return collect(new_copy.inner)
end

#=------------------------------------------------------------------------------
This section contains the `TileGame` data structure and functions that operate
on that data structure. The `TileGame` models and provides an interface for
an ongoing game, tracking the circle of tiles laid, the number of players,
and the number of tiles played and remaining. Functions are provided for
playing one or all rounds of the game and calculating the scores given the
current state of the game.
------------------------------------------------------------------------------=#

"""
A `TileGame` represents the current state of an ongoing game, including
the circle of tiles, the number of players, the number of tiles played,
and the number of tiles left to be played.
"""
mutable struct TileGame
    circle::TileCircle
    players::Int
    played::Int
    remaining::Int
end

"""
    TileGame(players::Int, tiles::Int)

Constructor that takes a number of players and number of tiles, creating
a `TileGame`.
"""
function TileGame(players::Int, tiles::Int)
    circle = TileCircle()
    return TileGame(circle, players, 0, tiles)
end

"""
    play_round!(game::TileGame)

Given a `TileGame`, play a single round, laying the next tile into the
circle and updating the number of tiles played and the number of tiles
remaining.
"""
function play_round!(game::TileGame)
    game.remaining <= 0 && return
    tile = Tile(game.played, game.played % game.players)
    push!(game.circle, tile)
    game.played += 1
    game.remaining -= 1
end

"""
    play_all!(game::TileGame)

Given a `TileGame`, play all the remaining rounds, updating the game to
its final state with all tiles in the circle and none remaining to be played.
"""
function play_all!(game::TileGame)
    while game.remaining > 0
        play_round!(game)
    end
end

"""
    scores((; circle, players)::TileGame)

Calculate and return the list of scores for each player given the current
state of a `TileGame`. Returns a list of scores, in order by player ID.
"""
function scores((; circle, players)::TileGame)
    score_list = Dict()
    for (idx, tile) in enumerate(collect(circle))
        (; number, played_by) = tile
        position      = idx - 1
        current_score = get!(score_list, played_by, 0)
        score_list[played_by] = current_score + (number * position)
    end
    return score_list
end

#=------------------------------------------------------------------------------
This section contains the function to complete the puzzle, demonstrating
the ergonomics of the previously implemented data structures.
------------------------------------------------------------------------------=#

"""
    game_winning_score(players::Int, tiles::Int)

Given the number of players and tiles, return the highest score at the
end of the game.
"""
function game_winning_score(players::Int, tiles::Int)
    game = TileGame(players, tiles)
    play_all!(game)                                 # 1, Lay all the tiles
    scorecard = scores(game)                        # 2. Determine player scores
    return reduce(max, values(scorecard), init = 0) # 3. Return maximum score
end

#=------------------------------------------------------------------------------
This section contains the tests to verify the operation of the main function,
proving puzzle completion.
------------------------------------------------------------------------------=#

using Test

@testset "Should return zero when there are no tiles." begin
    @test game_winning_score(5, 0) == 0
end

@testset "Should return high score for a one-player game" begin
    @test game_winning_score(1, 10) == 211
end

@testset "Should return high score for a two-player game" begin
    @test game_winning_score(2, 15) == 424
end

@testset "Should return high score for a many-player game" begin
    @test game_winning_score(100, 5000) == 314997735
end

max_stars_repo_name: BsideAlgorithms/CodingDojo

max_issues_repo_name: BsideAlgorithms/CodingDojo

max_forks_repo_name: BsideAlgorithms/CodingDojo


text:
@testset "Kalman filter" begin
    @testset "Time update inplace functions" begin
        x = randn(2)
        F = randn(2,2)
        x_apri = similar(x)
        @test KalmanFilter.calc_apriori_state!(x_apri, x, F) ≈ F * x

        Q = randn(2,2)
        P = randn(2,2)
        FP = similar(P)
        P_apri = similar(P)
        @test KalmanFilter.calc_apriori_covariance!(P_apri, FP, P, F, Q) ≈ F * P * F' .+ Q
    end

    @testset "Measurement update inplace functions" begin

        x = randn(2)
        y = randn(2)
        H = randn(2,2)
        ỹ = similar(y)
        @test KalmanFilter.calc_innovation!(ỹ, H, x, y) ≈ y .- H * x

        H = randn(2,2)
        R = randn(2,2)
        PHᵀ = randn(2,2)
        S = similar(R)
        @test KalmanFilter.calc_innovation_covariance!(S, H, PHᵀ, R) ≈ H * PHᵀ .+ R

        PHᵀ = randn(2,2)
        SL = randn(2,2)
        S = SL*SL'
        S_lu = similar(S)
        K = similar(PHᵀ)
        @test KalmanFilter.calc_kalman_gain!(S_lu, K, PHᵀ, S) ≈ PHᵀ / S

        x = randn(2)
        x_posterior = similar(x)
        K = randn(2,2)
        ỹ = randn(2)
        @test KalmanFilter.calc_posterior_state!(x_posterior, x, K, ỹ) ≈ x .+ K * ỹ

        P = randn(2,2)
        P_posterior = similar(P)
        PHᵀ = randn(2,2)
        K = randn(2,2)
        @test KalmanFilter.calc_posterior_covariance!(P_posterior, P, PHᵀ, K) ≈ P .- PHᵀ * K'

    end

    @testset "Time update" begin

        x = [1., 1.]
        P = [1. 0.; 0. 1.]
        F = [1. 0.; 0. 2.]
        Q = [1. 0.; 0. 1.]

        tu = time_update(x, P, F, Q)
        @test get_state(tu) == [1., 2.]
        @test get_covariance(tu) == [2. 0.; 0. 5.]

        tu_inter = KFTUIntermediate(2)
        tu = time_update!(tu_inter, x, P, F, Q)
        @test get_state(tu) == [1., 2.]
        @test get_covariance(tu) == [2. 0.; 0. 5.]
    end

    @testset "Measurement update" begin

        y = [1., 1.]
        x = [1., 1.]
        P = [1. 0.; 0. 1.]
        H = [1. 0.; 0. 1.]
        R = [1. 0.; 0. 1.]

        mu = measurement_update(x, P, y, H, R)
        @test get_state(mu) == [1., 1.]
        @test get_covariance(mu) == [0.5 0.; 0. 0.5]
        @test get_innovation(mu) == [0.0, 0.0]
        @test get_innovation_covariance(mu) == [2.0 0.0; 0.0 2.0]
        @test get_kalman_gain(mu) == [0.5 0.0; 0.0 0.5]

        mu_inter = KFMUIntermediate(2,2)
        mu = measurement_update!(mu_inter, x, P, y, H, R)
        @test get_state(mu) == [1., 1.]
        @test get_covariance(mu) ≈ [0.5 0.; 0. 0.5]
        @test get_innovation(mu) == [0.0, 0.0]
        @test get_innovation_covariance(mu) == [2.0 0.0; 0.0 2.0]
        @test get_kalman_gain(mu) ≈ [0.5 0.0; 0.0 0.5]

        y = 1.
        x = [1., 1.]
        P = [1. 0.; 0. 1.]
        H = [1., 0.]'
        R = 1.

        mu = measurement_update(x, P, y, H, R)
        @test get_state(mu) == [1., 1.]
        @test get_covariance(mu) == [0.5 0.; 0. 1.]
        @test get_innovation(mu) == 0.
        @test get_innovation_covariance(mu) == 2.
        @test get_kalman_gain(mu) == [0.5, 0.]

        H_temp = Matrix(H)
        y_temp = [y]
        R_temp = [R]
        mu_inter = KFMUIntermediate(2,1)
        mu = measurement_update!(mu_inter, x, P, y_temp, H_temp, R_temp)
        @test get_state(mu) == [1., 1.]
        @test get_covariance(mu) ≈ [0.5 0.; 0. 1.]
        @test get_innovation(mu) == [0.]
        @test get_innovation_covariance(mu) ≈ [2.]
        @test get_kalman_gain(mu) ≈ [0.5, 0.]

        y = [1., 1.]
        x = 1.
        P = 1.
        H = [1., 0.]
        R = [1. 0.; 0. 1.]

        mu = measurement_update(x, P, y, H, R)
        @test get_state(mu) == 1.
        @test get_covariance(mu) == 0.5
        @test get_innovation(mu) == [0.0, 1.0]
        @test get_innovation_covariance(mu) == [2.0 0.0; 0.0 1.0]
        @test get_kalman_gain(mu) == [0.5 0.0]

        x_temp = [x]
        P_temp = Matrix(adjoint([P]))
        H_temp = Matrix{eltype(H)}(undef, 2, 1)
        H_temp[1] = H[1]
        H_temp[2] = H[2]
        mu_inter = KFMUIntermediate(1,2)
        mu = measurement_update!(mu_inter, x_temp, P_temp, y, H_temp, R)
        @test get_state(mu) == [1.]
        @test get_covariance(mu) ≈ [0.5]
        @test get_innovation(mu) == [0.0, 1.0]
        @test get_innovation_covariance(mu) == [2.0 0.0; 0.0 1.0]
        @test get_kalman_gain(mu) ≈ [0.5 0.0]
    end
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 8fb71700-b4a0-11eb-1a77-57de26bec3ff
begin
	using Pkg,PlutoUI,Plots
	plotlyjs()
	Pkg.add("ForwardDiff")
	using ForwardDiff
end

# ╔═╡ 4d1975a9-08d7-43d3-86aa-5a1604a477f7
f(x) = x^2

# ╔═╡ 5474a79c-5e5a-4afe-92db-2614cc32f8a0
md"""
x₀ = $(@bind x0 Slider(-3:0.1:3,show_value=true,default=0))

Mostrar derivada: $(@bind mostrar_derivada CheckBox(default=false))
Mostrar recta tangente: $(@bind mostrar_tangente CheckBox(default=false))
Mostrar recta normal: $(@bind mostrar_normal CheckBox(default=false))
"""

# ╔═╡ a0ac637a-2b32-4568-9ffa-3ad5c69dbbd1
begin
	y0 = f(x0)
	df(x) = ForwardDiff.derivative(f,x)
	l(x) = df(x0)*(x-x0) + f(x0)
	n(x) = (-1/df(x0))*(x-x0) + f(x0)
	p = plot(f,label="f(x)",ylim=[-3,3],color=:blue)
	if mostrar_derivada
	plot!(df,label="f'(x)",color=:green)
	end
	if mostrar_tangente
	plot!(l,label=false,color=:red)
	end
	if mostrar_normal
	plot!(n,label=false,color=:black)
	end
	scatter!([(x0,y0)],label=false)
	p
end

# ╔═╡ e309ed4b-abf7-47e9-b79d-4db213eacbb0
md"""
Un par que preguntas a revisar cuando analicen una funcion

- Que signo tiene la derivada en ese punto?

- Cuales son los ceros y los conjuntos de positividad y negatividad de la derivada?

- Cuales son los maximos y minimos de la funcion?

- Cuales son los intervalis de crecimiento y decrecimiento de la funcion?
"""

# ╔═╡ Cell order:
# ╟─8fb71700-b4a0-11eb-1a77-57de26bec3ff
# ╠═4d1975a9-08d7-43d3-86aa-5a1604a477f7
# ╟─5474a79c-5e5a-4afe-92db-2614cc32f8a0
# ╟─a0ac637a-2b32-4568-9ffa-3ad5c69dbbd1
# ╟─e309ed4b-abf7-47e9-b79d-4db213eacbb0

max_stars_repo_name: akielbowicz/cbc-analisis-ingenieria

max_issues_repo_name: akielbowicz/cbc-analisis-ingenieria

max_forks_repo_name: akielbowicz/cbc-analisis-ingenieria


text:
################################################################################
#                   CORE DISPATCHVARIABLEREF METHOD EXTENSIONS
################################################################################
# Extend dispatch_variable_ref
function dispatch_variable_ref(model::InfiniteModel,
                               index::HoldVariableIndex
                               )::HoldVariableRef
    return HoldVariableRef(model, index)
end

# Extend _add_data_object
function _add_data_object(model::InfiniteModel,
                          object::VariableData{<:HoldVariable}
                          )::HoldVariableIndex
    return MOIUC.add_item(model.hold_vars, object)
end

# Extend _data_dictionary (type based)
function _data_dictionary(model::InfiniteModel, ::Type{HoldVariable}
    )::MOIUC.CleverDict{HoldVariableIndex, VariableData{HoldVariable{GeneralVariableRef}}}
    return model.hold_vars
end

# Extend _data_dictionary (reference based)
function _data_dictionary(vref::HoldVariableRef
    )::MOIUC.CleverDict{HoldVariableIndex, VariableData{HoldVariable{GeneralVariableRef}}}
    return JuMP.owner_model(vref).hold_vars
end

# Extend _data_object
function _data_object(vref::HoldVariableRef
    )::VariableData{HoldVariable{GeneralVariableRef}}
  object = get(_data_dictionary(vref), JuMP.index(vref), nothing)
  object === nothing && error("Invalid hold variable reference, cannot find " *
                        "corresponding variable in the model. This is likely " *
                        "caused by using the reference of a deleted variable.")
  return object
end

# Extend _core_variable_object
function _core_variable_object(vref::HoldVariableRef
    )::HoldVariable{GeneralVariableRef}
    return _data_object(vref).variable
end

################################################################################
#                          DEFINTION HELPER METHODS
################################################################################
# Check that parameter_bounds argument is valid
function _check_bounds(bounds::ParameterBounds{GeneralVariableRef};
                       _error = error)::Nothing
    depend_counter = Dict{DependentParameterRef, Int}()
    for (pref, set) in bounds
        # check that pref is an infinite parameter
        if !(_index_type(pref) <: InfiniteParameterIndex)
            _error("Can only specify infinite parameters for parameter bounds.")
        end
        # check that respects lower bound
        if JuMP.has_lower_bound(pref) && (JuMP.lower_bound(set) < JuMP.lower_bound(pref))
                _error("Specified parameter lower bound exceeds that defined " *
                       "for $pref.")
        end
        # check that respects upper bound
        if JuMP.has_upper_bound(pref) && (JuMP.upper_bound(set) > JuMP.upper_bound(pref))
                _error("Specified parameter upper bound exceeds that defined " *
                       "for $pref.")
        end
        # keep track of dependent parameters with equality conditions to ensure completeness
        if (_index_type(pref) == DependentParameterIndex) &&
           (JuMP.lower_bound(set) == JuMP.upper_bound(set))
            index = DependentParameterIndex(DependentParametersIndex(_raw_index(pref)), 1)
            dumby_pref = dispatch_variable_ref(pref.model, index)
            if haskey(depend_counter, dumby_pref)
                depend_counter[dumby_pref] += 1
            else
                depend_counter[dumby_pref] = 1
            end
        end
    end
    # check dimensions of dependent parameter equalities (if any are provided)
    for (pref, count) in depend_counter
        if _num_parameters(pref) != count
            _error("Cannot specify equality parameter bounds for a subset of " *
                   "dependent infinite parameters.")
        end
    end
    return
end

# Define _make_variable (used by JuMP.build_variable)
function _make_variable(_error::Function, info::JuMP.VariableInfo, ::Type{Hold};
    parameter_bounds::ParameterBounds{GeneralVariableRef} = ParameterBounds(),
    extra_kw_args...)::HoldVariable{GeneralVariableRef}
    # check for unneeded keywords
    for (kwarg, _) in extra_kw_args
        _error("Keyword argument $kwarg is not for use with hold variables.")
    end
    # check that the bounds don't violate parameter domains
    _check_bounds(parameter_bounds)
    # make variable and return
    return HoldVariable(_make_float_info(info), parameter_bounds)
end

# Validate parameter bounds and add support(s) if needed
function _validate_bounds(model::InfiniteModel,
                          bounds::ParameterBounds{GeneralVariableRef};
                          _error = error)::Nothing
    depend_supps = Dict{DependentParametersIndex, Matrix{Float64}}()
    for (pref, set) in bounds
        # check validity
        JuMP.check_belongs_to_model(pref, model)
        # ensure has a support if a point constraint was given
        if (JuMP.lower_bound(set) == JuMP.upper_bound(set))
            if _index_type(pref) == IndependentParameterIndex
                # label will be UserDefined
                add_supports(pref, JuMP.lower_bound(set), check = false)
            else
                index = DependentParametersIndex(_raw_index(pref))
                if !haskey(depend_supps, index)
                    dumby_pref = dispatch_variable_ref(model, DependentParameterIndex(index, 1))
                    depend_supps[index] = Matrix{Float64}(undef, _num_parameters(dumby_pref), 1)
                end
                depend_supps[index][_param_index(pref)] = JuMP.lower_bound(set)
            end
        end
    end
    # add dependent supports if any are given
    for (index, supp) in depend_supps
        prefs = [dispatch_variable_ref(model, DependentParameterIndex(index, i))
                 for i in 1:length(supp)]
        add_supports(prefs, supp, check = false) # label will be UserDefined
    end
    return
end

# Define _check_and_make_variable_ref (used by JuMP.add_variable)
function _check_and_make_variable_ref(model::InfiniteModel,
                                      v::HoldVariable,
                                      name::String)::HoldVariableRef
    _validate_bounds(model, v.parameter_bounds)
    data_object = VariableData(v, name)
    vindex = _add_data_object(model, data_object)
    vref = HoldVariableRef(model, vindex)
    if !isempty(v.parameter_bounds)
        model.has_hold_bounds = true
    end
    return vref
end

################################################################################
#                           VARIABLE INFO METHODS
################################################################################
# Set info for hold variables
function _update_variable_info(vref::HoldVariableRef,
                               info::JuMP.VariableInfo)::Nothing
    bounds = _core_variable_object(vref).parameter_bounds
    _set_core_variable_object(vref, HoldVariable(info, bounds))
    return
end

################################################################################
#                           PARAMETER BOUND METHODS
################################################################################
"""
    parameter_bounds(vref::HoldVariableRef)::ParameterBounds

Return the [`ParameterBounds`](@ref) object associated with the hold variable
`vref`. It contains a dictionary where each key is a parameter reference
which points to an `IntervalSet` that that defines a sub-domain for `vref`
relative to that parameter reference.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel())
julia> @infinite_parameter(model, t in [0, 10])
t

julia> @hold_variable(model, vref, parameter_bounds = (t in [0, 2]))
vref

julia> parameter_bounds(vref)
Subdomain bounds (1): t ∈ [0, 2]
```
"""
function parameter_bounds(vref::HoldVariableRef
    )::ParameterBounds{GeneralVariableRef}
    return _core_variable_object(vref).parameter_bounds
end

"""
    has_parameter_bounds(vref::HoldVariableRef)::Bool

Return a `Bool` indicating if `vref` is limited to a sub-domain as defined
by parameter bound.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel())
julia> @infinite_parameter(model, t in [0, 10])
t

julia> @hold_variable(model, vref, parameter_bounds = (t in [0, 2]))
vref

julia> has_parameter_bounds(vref)
true
```
"""
function has_parameter_bounds(vref::HoldVariableRef)::Bool
    return !isempty(parameter_bounds(vref))
end

# Other variable types
function has_parameter_bounds(vref::DispatchVariableRef)::Bool
    return false
end

# Internal function used to change the parameter bounds of a hold variable
function _update_variable_param_bounds(vref::HoldVariableRef,
                                       bounds::ParameterBounds{GeneralVariableRef}
                                       )::Nothing
    info = _variable_info(vref)
    _set_core_variable_object(vref, HoldVariable(info, bounds))
    JuMP.owner_model(vref).has_hold_bounds = true
    if is_used(vref)
        set_optimizer_model_ready(JuMP.owner_model(vref), false)
    end
    return
end

# Check that the bounds dictionary is compadable with existing dependent measures
function _check_meas_bounds(bounds::ParameterBounds{GeneralVariableRef},
                            data::AbstractMeasureData;
                            _error::Function = error
                            )::Nothing
    if !measure_data_in_hold_bounds(data, bounds)
        _error("New bounds don't span existing dependent measure bounds.")
    end
    return
end

# Update the current bounds to overlap with the new bounds if possible
function _update_bounds(bounds1::ParameterBounds{GeneralVariableRef},
                        bounds2::ParameterBounds{GeneralVariableRef};
                        _error::Function = error
                        )::Nothing
    # check each new bound
    for (pref, set) in bounds2
        # we have a new bound
        if !haskey(bounds1, pref)
            bounds1[pref] = set
        # the previous set and the new one do not overlap
        elseif (JuMP.lower_bound(set) > JuMP.upper_bound(bounds1[pref])) ||
               (JuMP.upper_bound(set) < JuMP.lower_bound(bounds1[pref]))
            _error("Sub-domains of constraint and/or hold variable(s) do not" *
                   " overlap. Consider changing the parameter bounds of the" *
                   " constraint and/or hold variable(s).")
        # we have an existing bound
        else
            # we have a new stricter lower bound to update with
            if JuMP.lower_bound(set) > JuMP.lower_bound(bounds1[pref])
                bounds1[pref] = IntervalSet(JuMP.lower_bound(set),
                                            JuMP.upper_bound(bounds1[pref]))
            end
            # we have a new stricter upper bound to update with
            if JuMP.upper_bound(set) < JuMP.upper_bound(bounds1[pref])
                bounds1[pref] = IntervalSet(JuMP.lower_bound(bounds1[pref]),
                                            JuMP.upper_bound(set))
            end
        end
    end
    return
end

## Update the variable bounds if it has any
# DispatchVariableRef
function _update_var_bounds(vref::DispatchVariableRef,
                            constr_bounds::ParameterBounds{GeneralVariableRef}
                            )::Nothing
    return
end

# HoldVariableRef
function _update_var_bounds(vref::HoldVariableRef,
                            constr_bounds::ParameterBounds{GeneralVariableRef}
                            )::Nothing
    if has_parameter_bounds(vref)
        _update_bounds(constr_bounds, parameter_bounds(vref))
    end
    return
end

# MeasureRef
function _update_var_bounds(mref::MeasureRef,
                            constr_bounds::ParameterBounds{GeneralVariableRef}
                            )::Nothing
    vrefs = _all_function_variables(measure_function(mref))
    for vref in vrefs
        _update_var_bounds(vref, constr_bounds)
    end
    return
end

# GeneralVariableRef
function _update_var_bounds(vref::GeneralVariableRef,
                            constr_bounds::ParameterBounds{GeneralVariableRef}
                            )::Nothing
    return _update_var_bounds(dispatch_variable_ref(vref), constr_bounds)
end

## Rebuild the constraint bounds (don't change in case of error)
# BoundedScalarConstraint
function _rebuild_constr_bounds(c::BoundedScalarConstraint,
                                var_bounds::ParameterBounds{GeneralVariableRef};
                                _error::Function = error
                                )::JuMP.AbstractConstraint
    # prepare new constraint
    func = JuMP.jump_function(c)
    set = JuMP.moi_set(c)
    vrefs = _all_function_variables(func)
    orig_bounds = original_parameter_bounds(c)
    c_new = BoundedScalarConstraint(func, set, copy(orig_bounds), orig_bounds)
    # look for bounded hold variables and update bounds
    for vref in vrefs
        _update_var_bounds(vref, parameter_bounds(c_new))
    end
    # check if the constraint has bounds and update if doesn't
    if isempty(parameter_bounds(c_new))
        c_new = JuMP.ScalarConstraint(func, set)
    end
    return c_new
end

# ScalarConstraint
function _rebuild_constr_bounds(c::JuMP.ScalarConstraint,
                                var_bounds::ParameterBounds{GeneralVariableRef};
                                _error::Function = error
                                )::BoundedScalarConstraint
    func = JuMP.jump_function(c)
    set = JuMP.moi_set(c)
    return BoundedScalarConstraint(func, set, var_bounds, ParameterBounds())
end

"""
    set_parameter_bounds(vref::HoldVariableRef,
                         bounds::ParameterBounds{GeneralVariableRef};
                         [force::Bool = false])::Nothing

Specify a new dictionary of parameter bounds `bounds` for the hold variable `vref`.
These are stored in a [`ParameterBounds`](@ref) object which contains a dictionary.
Note the dictionary keys must be infinite parameter references and the values
must be `IntervalSet`s that indicate a particular sub-domain for which `vref`
is defined. This is meant to be primarily used by
[`@set_parameter_bounds`](@ref) which provides a more intuitive syntax.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel())
julia> @infinite_parameter(model, t in [0, 10])
t

julia> @hold_variable(model, vref)
vref

julia> set_parameter_bounds(vref, ParameterBounds(Dict(t => IntervalSet(0, 2))))

julia> parameter_bounds(vref)
Subdomain bounds (1): t ∈ [0, 2]
```
"""
function set_parameter_bounds(vref::HoldVariableRef,
                              bounds::ParameterBounds{GeneralVariableRef};
                              force::Bool = false, _error::Function = error
                              )::Nothing
    if has_parameter_bounds(vref) && !force
        _error("$vref already has parameter bounds. Consider adding more using " *
               "`add_parameter_bounds` or overwriting them by setting " *
               "the keyword argument `force = true`")
    else
        # check that bounds are valid
        _check_bounds(bounds, _error = _error)
        # check dependent measures
        model = JuMP.owner_model(vref)
        cindices = ConstraintIndex[]
        for mindex in _measure_dependencies(vref)
            mref = dispatch_variable_ref(model, mindex)
            _check_meas_bounds(bounds, _core_variable_object(mref).data,
                               _error = _error)
            append!(cindices, _constraint_dependencies(mref))
        end
        # set the new bounds
        _validate_bounds(model, bounds, _error = _error)
        _update_variable_param_bounds(vref, bounds)
        # check and update dependent constraints
        union!(cindices, _constraint_dependencies(vref))
        for cindex in cindices
            cref = _temp_constraint_ref(model, cindex)
            constr = _core_constraint_object(cref)
            new_constr = _rebuild_constr_bounds(constr, bounds, _error = _error)
            _set_core_constraint_object(cref, new_constr)
        end
    end
    return
end

## Check and update the constraint bounds (don't change in case of error)
# BoundedScalarConstraint
function _update_constr_bounds(bounds::ParameterBounds{GeneralVariableRef},
                               c::BoundedScalarConstraint;
                               _error::Function = error
                               )::BoundedScalarConstraint
    new_bounds = copy(parameter_bounds(c))
    _update_bounds(new_bounds, bounds, _error = _error)
    func = JuMP.jump_function(c)
    set = JuMP.moi_set(c)
    orig_bounds = original_parameter_bounds(c)
    return BoundedScalarConstraint(func, set, new_bounds, orig_bounds)
end

# ScalarConstraint
function _update_constr_bounds(bounds::ParameterBounds{GeneralVariableRef},
                               c::JuMP.ScalarConstraint;
                               _error::Function = error
                               )::BoundedScalarConstraint
    func = JuMP.jump_function(c)
    set = JuMP.moi_set(c)
    return BoundedScalarConstraint(func, set, bounds, ParameterBounds())
end

"""
    add_parameter_bounds(vref::HoldVariableRef,
                         new_bounds::ParameterBounds{GeneralVariableRef}
                         )::Nothing

Add additional parameter bounds to `vref` such that it is defined over the
sub-domain based on the intersection of the existing bounds and `new_bounds`.
This is primarily meant to be used by [`@add_parameter_bounds`](@ref).

```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel())
julia> @infinite_parameter(model, t in [0, 10])
t

julia> @hold_variable(model, vref)
vref

julia> add_parameter_bounds(vref, ParameterBounds(t => IntervalSet(0, 2)))

julia> parameter_bounds(vref)
Subdomain bounds (1): t ∈ [0, 2]
```
"""
function add_parameter_bounds(vref::HoldVariableRef,
                              new_bounds::ParameterBounds{GeneralVariableRef};
                              _error::Function = error
                              )::Nothing
    # if overriding any existing parameter bounds, dispatch to set new ones
    if any(haskey(parameter_bounds(vref), pref) for pref in keys(new_bounds))
        merged_bounds = merge(parameter_bounds(vref), new_bounds)
        return set_parameter_bounds(vref, merged_bounds, force = true,
                                    _error = _error)
    end
    # check the new bounds
    _check_bounds(new_bounds, _error = _error)
    # check dependent measures
    model = JuMP.owner_model(vref)
    cindices = ConstraintIndex[]
    for mindex in _measure_dependencies(vref)
        mref = dispatch_variable_ref(model, mindex)
        _check_meas_bounds(new_bounds, _core_variable_object(mref).data,
                           _error = _error)
        append!(cindices, _constraint_dependencies(mref))
    end
    # check and update dependent constraints
    union!(cindices, _constraint_dependencies(vref))
    for cindex in cindices
        cref = _temp_constraint_ref(model, cindex)
        constr = _core_constraint_object(cref)
        new_constr = _update_constr_bounds(new_bounds, constr, _error = _error)
        _set_core_constraint_object(cref, new_constr)
    end
    _validate_bounds(model, new_bounds, _error = _error)
    # add the bounds
    merge!(parameter_bounds(vref), new_bounds)
    # update status
    JuMP.owner_model(vref).has_hold_bounds = true
    if is_used(vref)
        set_optimizer_model_ready(model, false)
    end
    return
end

"""
    delete_parameter_bounds(vref::HoldVariableRef)::Nothing

Delete all the parameter bounds of the hold variable `vref`. Any constraints
that employ `vref` will be updated accordingly.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel())
julia> @infinite_parameter(model, x[1:2] in [0, 10])
2-element Array{ParameterRef,1}:
 x[1]
 x[2]

julia> @hold_variable(model, z, parameter_bounds = (x in [0, 1]))
z

julia> delete_parameter_bounds(z)

julia> parameter_bounds(z)
Subdomain bounds (0):
```
"""
function delete_parameter_bounds(vref::HoldVariableRef)::Nothing
    # get the current bounds
    bounds = parameter_bounds(vref)
    # check if there are bounds and act accordingly
    if !isempty(bounds)
        _update_variable_param_bounds(vref, ParameterBounds())
        # check for dependent measures that are used by constraints
        model = JuMP.owner_model(vref)
        cindices = ConstraintIndex[]
        for mindex in _measure_dependencies(vref)
            mref = dispatch_variable_ref(model, mindex)
            _check_meas_bounds(bounds, _core_variable_object(mref).data,
                               _error = error)
            append!(cindices, _constraint_dependencies(mref))
        end
        # check and update dependent constraints
        union!(cindices, _constraint_dependencies(vref))
        for cindex in cindices
            cref = _temp_constraint_ref(model, cindex)
            constr = _core_constraint_object(cref)
            new_constr = _rebuild_constr_bounds(constr, bounds, _error = error)
            _set_core_constraint_object(cref, new_constr)
        end
        # update status
        if is_used(vref)
            set_optimizer_model_ready(model, false)
        end
    end
    return
end

################################################################################
#                                 DELETION
################################################################################
# Extend _delete_variable_dependencies (for use with JuMP.delete)
function _delete_variable_dependencies(vref::HoldVariableRef)::Nothing
    # remove variable info constraints associated with vref
    _delete_info_constraints(vref)
    # update dependent constraint bounds by deleting the variable bounds
    delete_parameter_bounds(vref)
    return
end


text:
"""
From: CMES, vol.49, no.2, pp.81-110, 2009

The problem considered in this section is that of a hyperbolic paraboloid shell,
clamped along one side and free on three edges and loaded by self-weight
(Figure 16). This is a pure bending dominated problem and known to be a very
hard test for locking behaviour as suggested in References (Chapelle and Bathe,
1998; Bathe, Iosilevich, and Chapelle, 2000).

Table 2 in Bathe, Iosilevich, and Chapelle, 2000 (computed with MITC16).
t/L   Strain energy  Displacement
1/100 1.6790e-3 9.3355e-5
1/1000 1.1013e-2 6.3941e-3
1/10,000 8.9867e-2 5.2988e-1

See also the table in
An Improved Quadrilateral Flat Element with Drilling
Degrees of Freedom for Shell Structural Analysis
H. Nguyen-Van1 , N. Mai-Duy1 and T. Tran-Cong1
CMES, vol.49, no.2, pp.81-110, 2009.

The shell geometry is described by the
equation: z = x^2 −y^2 ; (x,y) ∈ − L/2 ; L/2
"""
module clamped_hypar_examples

using LinearAlgebra
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FESetShellQ4Module: FESetShellQ4
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json
using FinEtools.MeshExportModule.VTKWrite: vtkwrite

function _execute_full(tL_ratio = 1/100, g = 80*0.1^0, analyt_sol=-9.3355e-5, n = 32, visualize = false)
    # analytical solution for the vertical deflection and the midpoint of the
    # free edge
    # Parameters:
    E = 2.0e11;
    nu = 0.3;
    L = 1.0
    thickness = tL_ratio * L
    # Bathe, Iosilevich, and Chapelle (2000) with a refined mesh of
    # high-order element MITC16
    formul = FEMMShellT3FFModule

    @info "Mesh: $n elements per side"

    tolerance = L/n/1000
    fens, fes = T3block(L,L,n,n);
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        x=fens.xyz[i, 1]-L/2; y=fens.xyz[i, 2]-L/2;
        fens.xyz[i, :] .= (x, y, x^2 - y^2)
    end

    mater = MatDeforElastIso(DeforModelRed3D, E, nu)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    # femm.drilling_stiffness_scale = 0.1
    # femm.mult_el_size = 0.2
    stiffness = formul.stiffness
    associategeometry! = formul.associategeometry!

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Clamped edge
    l1 = selectnode(fens; box = Float64[-L/2 -L/2 -Inf Inf -Inf Inf], inflate = tolerance)
    for i in [1,2,3,4,5,6]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    # Midpoint of the free edge
    nl = selectnode(fens; box = Float64[L/2 L/2 0 0 -Inf Inf], inflate = tolerance)
    lfemm = FEMMBase(IntegDomain(fes, TriRule(1)))
    # computeforce!(forceout, XYZ, tangents, fe_label) = let
    #     n = cross(tangents[:, 1], tangents[:, 2])
    #     n = n / norm(n)
    #     forceout .= 0.0
    #     forceout[3] = -g * n[3]
    #     return forceout
    # end
    # fi = ForceIntensity(FFlt, 6, computeforce!)
    fi = ForceIntensity(FFlt[0, 0, -g, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 3);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    targetu = dchi.values[nl, 3][1]
    @info "Solution: $(round(targetu/analyt_sol, digits = 4)*100)%"

        # Visualization
    if visualize
        scattersysvec!(dchi, (L/4)/abs(targetu).*U)
        update_rotation_field!(Rfield0, dchi)
        plots = cat(plot_space_box([[0 0 -L/2]; [L/2 L/2 L/2]]),
        #plot_nodes(fens),
            plot_midsurface(fens, fes; x = geom0.values, facecolor = "rgb(12, 12, 123)"),
            plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
            dims = 1)
        pl = render(plots)
    end
    return targetu
end

function _execute_half(orientation = :a, tL_ratio = 1/100, g = 80*0.1^0, analyt_sol=-9.3355e-5, n = 32, visualize = false)
    # analytical solution for the vertical deflection and the midpoint of the
    # free edge
    # Parameters:
    E = 2.0e11;
    nu = 0.3;
    L = 1.0
    thickness = tL_ratio * L
    # Bathe, Iosilevich, and Chapelle (2000) with a refined mesh of
    # high-order element MITC16
    formul = FEMMShellT3FFModule

    @info "Mesh: $n elements per side"

    tolerance = L/n/1000
    fens, fes = T3block(L,L/2,n,Int(round(n/2)),orientation);
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        x=fens.xyz[i, 1]-L/2; y=fens.xyz[i, 2]-L/2;
        fens.xyz[i, :] .= (x, y, x^2 - y^2)
    end

    mater = MatDeforElastIso(DeforModelRed3D, E, nu)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    # femm.drilling_stiffness_scale = 0.1
    # femm.mult_el_size = 0.2
    stiffness = formul.stiffness
    associategeometry! = formul.associategeometry!

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Clamped edge
    l1 = selectnode(fens; box = Float64[-L/2 -L/2 -Inf Inf -Inf Inf], inflate = tolerance)
    for i in [1,2,3,4,5,6]
        setebc!(dchi, l1, true, i)
    end
    # Symmetry Plane
    l1 = selectnode(fens; box = Float64[-Inf Inf 0 0 -Inf Inf], inflate = tolerance)
    for i in [2,4,6]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    # Midpoint of the free edge
    nl = selectnode(fens; box = Float64[L/2 L/2 0 0 -Inf Inf], inflate = tolerance)
    lfemm = FEMMBase(IntegDomain(fes, TriRule(1)))
    fi = ForceIntensity(FFlt[0, 0, -g, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 3);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    targetu = dchi.values[nl, 3][1]
    @info "Solution: $(round(targetu/analyt_sol, digits = 4)*100)%"

        # Visualization
    if visualize

        vtkwrite("clamped_hypar-$(orientation)-$(n).vtu", fens, fes; vectors = [("u", dchi.values[:, 1:3])])

        scattersysvec!(dchi, (L/4)/abs(targetu).*U)
        update_rotation_field!(Rfield0, dchi)
        plots = cat(plot_space_box([[0 0 -L/2]; [L/2 L/2 L/2]]),
        #plot_nodes(fens),
            plot_midsurface(fens, fes; x = geom0.values, facecolor = "rgb(12, 12, 123)"),
            plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
            dims = 1)
        pl = render(plots)
    end
    return targetu
end

function test_convergence(orientation = :a)

    tL_ratios = [1/100, 1/1000, 1/10000];
    gs = [80*0.1^0, 80*0.1^1, 80*0.1^2]
    analyt_sols = [-9.3355e-5, -6.3941e-3, -5.2988e-1];

    ns = [4, 8, 16, 32, 64, 128, 256, 512]
    all_results = []
    for (tL_ratio, g, analyt_sol) in zip(tL_ratios, gs, analyt_sols)
        @info "Clamped hypar, t/L=$(tL_ratio)"
        results = Float64[]
        for n in ns
        # for n in [4, 8, 16, 32, 64, ]
            r = _execute_half(orientation, tL_ratio, g, analyt_sol, n, false)
            push!(results, r/analyt_sol)
        end
        push!(all_results, results)
    end

    return ns, all_results
end

function test_0_001(orientation = :a, ns = [8, 16, 32, 48, ])
    tL_ratio = 1/1000
    g = 80*0.1^1
    analyt_sol = -6.3941e-3

    @info "Clamped hypar, t/L=$(tL_ratio)"
    results = Float64[]
    for n in ns
        # for n in [4, 8, 16, 32, 64, ]
        r = _execute_half(orientation, tL_ratio, g, analyt_sol, n, false)
        push!(results, r/analyt_sol)
    end

    return ns, results
end

# using Gnuplot

# ns = 1 ./ [4, 8, 16, 32, 64, 128, 256, ]
# results = [0.8839348674712617, 0.8750949157612452, 0.9199805802913757, 0.9619508790573108, 0.9856803572479892, 0.9955727622499687, 0.9993169031485688]
# @gp ns results "with lp"      :-
# results = [1.0429797613488236, 0.9314984628085947, 0.9365905225801154, 0.9565506281799385, 0.9764476699285441, 0.9902805329751646, 0.9968920296205528]
# @gp :- ns results "with lp"    :-
# results = [1.251888013432877, 1.0155533090845452, 0.9678060658415124, 0.9718188010061173, 0.9812934066246979, 0.9889499817887738, 0.9953521405300628]
# @gp :- ns results "with lp"

end # module

using PGFPlotsX

using .clamped_hypar_examples

let
    for orientation in (:a, :b)
        ns, all_results = clamped_hypar_examples.test_convergence(orientation)

        objects = []

        results = all_results[1]
        @pgf p = PGFPlotsX.Plot(
        {
        color = "black",
        line_width  = 0.7
        },
        Coordinates([v for v in  zip(1 ./ ns, results)])
        )
        push!(objects, p)
        push!(objects, LegendEntry("t/L=1/100"))

        results = all_results[2]
        @pgf p = PGFPlotsX.Plot(
        {
        color = "black",
        style = "dashed",
        line_width  = 0.7
        },
        Coordinates([v for v in  zip(1 ./ ns, results)])
        )
        push!(objects, p)
        push!(objects, LegendEntry("t/L=1/1000"))

        results = all_results[3]
        @pgf p = PGFPlotsX.Plot(
        {
        color = "black",
        style = "dotted",
        line_width  = 0.7
        },
        Coordinates([v for v in  zip(1 ./ ns, results)])
        )
        push!(objects, p)
        push!(objects, LegendEntry("t/L=1/10000"))

        @pgf ax = Axis(
        {
        xlabel = "Relative Element Size [ND]",
        ylabel = "Normalized Displacement [ND]",
        ymin = 0.7,
        ymax = 1.1,
        xmode = "log",
        ymode = "linear",
        yminorgrids = "true",
        grid = "both",
        legend_style = {
        at = Coordinate(0.5, 1.05),
        anchor = "south",
        legend_columns = -1
        },
        },
        objects...
        )

        display(ax)
        pgfsave("clamped_hypar_examples-dependence-on-t-L-$(orientation).pdf", ax)
    end
end

using PGFPlotsX

objects = []

# What is the source of these data points?
# Development of a strain-smoothed three-node triangular flat shell
# element with drilling degrees of freedom
# C.M. Shin, B.C. Lee / Finite Elements in Analysis and Design 86 (2014) 71–80
# For some reason the results there are scaled with -5.958e-3.

# Allman = [0.07, 0.3074, 0.829, 0.9545] .* (-5.958e-3) ./ (-6.3941e-3)
# Providas_Kattis = [1.07, 1.00, 0.9859, 0.9857] .* (-5.958e-3) ./ (-6.3941e-3)
# Cook_flat_stiffened = [0.5632, 0.9399, 0.9849, 0.9913] .* (-5.958e-3) ./ (-6.3941e-3)
# Shin_Lee = [1.0826 1.0144 0.9986 0.9991]  .* (-5.958e-3) ./ (-6.3941e-3)

# all_results = [("Present", ns, results, "*"), ("Allman", Allman, "x"), ("Providas, Kattis", Providas_Kattis, "triangle"), ("Cook, flat, stiffened", Cook_flat_stiffened, "square"), ("Shin, Lee", Shin_Lee, "diamond")]

# Here we take the reference results from Bathe, Iosilevich, Chapelle.

# An efficient three‑node triangular Mindlin–Reissner flat shell element
# Hosein Sangtarash1 · Hamed Ghohani Arab1 · Mohammad R. Sohrabi1 · Mohammad R. Ghasemi1
# TMRFS = ("TMRFS", [4, 8, 16, 24], [0.533 0.857 0.990 0.998], "diamond")
# MITC3p = ("MITC3+", [4, 8, 16, ], [0.407 0.768 0.93], "triangle")
# all_results = [("Present", ns, results, "*"), TMRFS, MITC3p]

# Performance of the MITC3+ and MITC4+ shell elements in widely-used
# benchmark problems
# Yeongbin Ko a, Youngyu Lee b, Phill-Seung Lee a,⇑, Klaus-Jürgen Bathe c
# a Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology, 291 Daehak-ro, Yuseong-gu,
MITC3p = ("MITC3+", [4, 8, 16, 32, 64], [0.9533 0.9589 0.9728 0.9868 0.9951], "diamond")

ns, resultsa = clamped_hypar_examples.test_0_001(:a, [4, 8, 16, 32, 64, 128, 256])
ns, resultsb = clamped_hypar_examples.test_0_001(:b, [4, 8, 16, 32, 64, 128, 256])
all_results = [("Present (a)", ns, resultsa, "*"), ("Present (b)", ns, resultsb, "o"), MITC3p]

let
    for r in  all_results
        @pgf p = PGFPlotsX.Plot(
        {
        color = "black",
        line_width  = 0.7,
        style = "solid",
        mark = "$(r[4])"
        },
        Coordinates([v for v in  zip(1 ./ (r[2]), abs.(1 .- r[3]))])
        )
        push!(objects, p)
        push!(objects, LegendEntry("$(r[1])"))
    end

    @pgf ax = Axis(
    {
    xlabel = "Relative element size [ND]",
    ylabel = "Normalized Displacement Error [ND]",
        # xmin = range[1],
        # xmax = range[2],
    xmode = "log",
    ymode = "log",
    yminorgrids = "true",
    grid = "both",
    legend_style = {
    at = Coordinate(0.5, 1.05),
    anchor = "south",
    legend_columns = -1
    },
    },
    objects...
    )

    display(ax)
    pgfsave("clamped_hypar_examples-0_001-errors.pdf", ax)
end

using .clamped_hypar_examples
ns, all_results = clamped_hypar_examples.test_0_001(:a, 4 .* [24, 48, 96])
q1, q2, q3 = all_results[:]
qtrue = (q2^2 - q1 * q3) / (2*q2 - q1 - q3)
@show (q1, q2, q3) .* -6.3941, qtrue .* -6.3941
ns, all_results = clamped_hypar_examples.test_0_001(:b, 4 .* [24, 48, 96])
q1, q2, q3 = all_results[:]
qtrue = (q2^2 - q1 * q3) / (2*q2 - q1 - q3)
@show (q1, q2, q3) .* -6.3941, qtrue .* -6.3941

q1, q2, q3 = MITC3p[3][3:end]
qtrue = (q2^2 - q1 * q3) / (2*q2 - q1 - q3)
@show (q1, q2, q3) .* -6.3941, qtrue .* -6.3941


text:
using Documenter
using Glob
using Literate

using Distributions
using HMMBase
using Random

Random.seed!(2019)

# https://discourse.julialang.org/t/deactivate-plot-display-to-avoid-need-for-x-server/19359
ENV["GKSwstype"] = "nul"

examples = []

if !("SKIP_EXAMPLES" in keys(ENV))
    examples = map(glob("*.jl", "examples/")) do x
        joinpath("examples/", split(basename(x), ".")[1])
    end
    for example in examples
        Literate.markdown("$(example).jl", "docs/src/examples", documenter = true)
    end
end

makedocs(
    sitename = "HMMBase",
    modules = [HMMBase],
    format = Documenter.HTML(
        assets = ["assets/goatcounter.js"]
    ),
    pages = [
        "index.md",
        "Manual" => ["basics.md", "models.md", "algorithms.md", "utilities.md"],
        "Examples" => map(example -> "$(example).md", examples),
        "internals.md",
        "migration.md",
        "_index.md",
    ],
)

deploydocs(repo = "github.com/maxmouchet/HMMBase.jl.git")