text:
using PatternDispatch

degree_scale(degree, scale) = error()
degree_scale(degree, scale, tuning) = error()

@patterns begin
  function (@inverse degree_scale(degree, scale::Scale))(e)
  	e::Event
  	@guard haskey(e, :degree)
  	@guard haskey(e, :scale)
  	degree = e[:degree]
  	scale = e[:scale]
  end
  function (@inverse degree_scale(degree, scale::Scale, tuning::Tuning))(e)
  	e::Event
  	@guard haskey(e, :degree)
  	@guard haskey(e, :scale)
  	@guard haskey(e, :tuning)
  	degree = e[:degree]
  	scale = e[:scale]
  	tuning = e[:tuning]
  end

  to_freq(x) = error("Cannot convert to frequency: $(x)")
  to_freq(degree_scale(d, s)) = dst_to_freq(d, s, Tuning(:edo_12))
  to_freq(degree_scale(d, s, t)) = dst_to_freq(d, s, t)
end

function dst_to_freq(d::Int, s::Scale, t::Tuning)
  # TODO make this more robust
  t[s[d]]
end


text:
using SubHunt
using Test
using Random
using POMDPs
using POMDPPolicies
using POMDPSimulators
using DiscreteValueIteration
using ParticleFilters
using POMDPModelTools
using QMDP

@testset "VI" begin
    rng = MersenneTwister(6)
    pomdp = SubHuntPOMDP()
    # show(STDOUT, MIME("text/plain"), SubVis(pomdp))

    solver = ValueIterationSolver(verbose=true)
    vipolicy = solve(solver, UnderlyingMDP(pomdp))
    #
    s = rand(rng, initialstate(pomdp))
    @show value(vipolicy, s)
    @show value(vipolicy, SubState(s.own, s.target, s.goal, true))
end

@testset "QMDP and PF" begin
    rng = MersenneTwister(6)
    pomdp = SubHuntPOMDP()
    policy = solve(QMDPSolver(verbose=true), pomdp)

    # policy = RandomPolicy(pomdp, rng=rng)

    filter = BootstrapFilter(pomdp, 10000, rng)

    for (s, a, r, sp) in stepthrough(pomdp, policy, filter, "s,a,r,sp", max_steps=200, rng=rng)
        v = SubVis(pomdp, s=s, a=a, r=r)
        show(stdout, MIME("text/plain"), v)
    end
end

@testset "DPOMDP and PF" begin
    rng = MersenneTwister(6)
    dpomdp = DSubHuntPOMDP(SubHuntPOMDP(), 1.0)
    policy = RandomPolicy(dpomdp, rng=rng)
    filter = BootstrapFilter(dpomdp, 10000, rng)

    for (s, a, r, sp) in stepthrough(dpomdp, policy, filter, "s,a,r,sp", max_steps=200, rng=rng)
        v = SubVis(dpomdp.cp, s=s, a=a, r=r)
        show(stdout, MIME("text/plain"), v)
    end
end

@testset "Visualization" begin
    rng = MersenneTwister(6)
    pomdp = SubHuntPOMDP()
    policy = solve(QMDPSolver(verbose=true), pomdp)

    # policy = RandomPolicy(pomdp, rng=rng)

    filter = BootstrapFilter(pomdp, 10000, rng)

    for step in stepthrough(pomdp, policy, filter, "s,a,r,sp", max_steps=200, rng=rng)
        show(stdout, MIME("text/plain"), render(pomdp, step))
    end
end

ext: jl


text:
immutable LogNormal <: ContinuousUnivariateDistribution
    meanlog::Float64
    sdlog::Float64
    function LogNormal(ml::Real, sdl::Real)
    	sdl > zero(sdl) || error("sdlog must be positive")
    	new(float64(ml), float64(sdl))
    end
end

LogNormal(ml::Real) = LogNormal(ml, 1.0)
LogNormal() = LogNormal(0.0, 1.0)

@_jl_dist_2p LogNormal lnorm

entropy(d::LogNormal) = 0.5 + 0.5 * log(2.0 * pi * d.sdlog^2) + d.meanlog

insupport(::LogNormal, x::Real) = zero(x) < x < Inf
insupport(::Type{LogNormal}, x::Real) = zero(x) < x < Inf

function kurtosis(d::LogNormal)
   exp(4.0 * d.sdlog^2) + 2.0 * exp(3.0 * d.sdlog^2) +
        3.0 * exp(2.0 * d.sdlog^2) - 6.0
end

mean(d::LogNormal) = exp(d.meanlog + d.sdlog^2 / 2)

median(d::LogNormal) = exp(d.meanlog)

# mgf(d::LogNormal)
# cf(d::LogNormal)

mode(d::LogNormal) = exp(d.meanlog - d.sdlog^2)
modes(d::LogNormal) = [mode(d)]

function skewness(d::LogNormal)
    (exp(d.sdlog^2) + 2.0) * sqrt(exp(d.sdlog^2) - 1.0)
end

function var(d::LogNormal)
    sigsq = d.sdlog^2
    (exp(sigsq) - 1) * exp(2d.meanlog + sigsq)
end

function fit_mle{T <: Real}(::Type{LogNormal}, x::Array{T})
    lx = log(x)
    LogNormal(mean(lx), std(lx))
end


text:
struct Handler{P<:AbstractPath}
    path::P
    settings        # Could be Vector or Pairs on 0.6 or 1.0 respectively
end

"""
    Handler(path::Union{String, AbstractPath}; kwargs...)
    Handler(bucket::String, prefix::String; kwargs...)

Handles iteratively saving JLSO file to the specified path location.
FilePath are used to abstract away differences between paths on S3 or locally.
"""
Handler(path::AbstractPath; kwargs...) = Handler(path, kwargs)
Handler(path::String; kwargs...) = Handler(Path(path), kwargs)
Handler(bucket::String, prefix::String; kwargs...) = Handler(S3Path("s3://$bucket/$prefix"), kwargs)

"""
    path(handler, name)

Determines the path to save to based on the handlers path prefix, name, and context.
Tags are used to dynamically prefix the named file with the handler's path.
Names with a '.' separators will be used to form subdirectories
(e.g., "Foo.bar.x" will be saved to "\$prefix/Foo/bar/x.jlso").
"""
function path(handler::Handler{P}, name::String) where P
    prefix = ["$key=$val" for (key,val) in CONTEXT_TAGS[]]

    # Split up the name by '.' and add the jlso extension
    parts = split(name, '.')
    parts[end] = string(parts[end], ".jlso")

    return join(handler.path, prefix..., parts...)
end

"""
    stage!(handler::Handler, jlso::JLSOFIle, data::Dict{Symbol})

Update the JLSOFile with the new data.
"""
function stage!(handler::Handler, jlso::JLSO.JLSOFile, data::Dict{Symbol})
    for (k, v) in data
        jlso[k] = v
    end

    return jlso
end

"""
    commit!(handler, path, jlso)

Write the JLSOFile to the path as bytes.
"""
function commit!(handler::Handler{P}, path::P, jlso::JLSO.JLSOFile) where P <: AbstractPath
    # NOTE: This is only necessary because FilePathsBase.FileBuffer needs to support
    # write(::FileBuffer, ::UInt8)
    # https://github.com/rofinn/FilePathsBase.jl/issues/45
    io = IOBuffer()
    write(io, jlso)
    bytes = take!(io)
    mkdir(parent(path); recursive=true, exist_ok=true)
    write(path, bytes)
end

function checkpoint(handler::Handler, name::String, data::Dict{Symbol}; tags...)
    checkpoint_deprecation(tags...)
    with_checkpoint_tags(tags...) do
        debug(LOGGER, "Checkpoint $name triggered, with context: $(join(CONTEXT_TAGS[], ", ")).")
        jlso = JLSO.JLSOFile(Dict{Symbol, Vector{UInt8}}(); handler.settings...)
        p = path(handler, name)
        stage!(handler, jlso, data)
        commit!(handler, p, jlso)
    end
end

#=
Define our no-op conditions just to be safe
=#
function checkpoint(handler::Nothing, name::String, data::Dict{Symbol}; tags...)
    checkpoint_deprecation(tags...)
    with_checkpoint_tags(tags...) do
        debug(LOGGER, "Checkpoint $name triggered, but no handler has been set.")
        nothing
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 98ed9a78-f50c-11ea-1d20-7bcaee689a99
using Pkg, DrWatson

# ╔═╡ 98edd6be-f50c-11ea-26c3-ff6ac7641060
begin
	@quickactivate "StatisticalRethinkingTuring"
	using StatisticalRethinking
end

# ╔═╡ 6e665478-f50c-11ea-1861-3f509f92ac5c
md"## Fig3.6t.jl"

# ╔═╡ 98ee52ec-f50c-11ea-22bc-c1075a60bfed
md"### snippet 3.23"

# ╔═╡ 98f67062-f50c-11ea-23f3-539ae4f7e67b
begin
	N = 10000
	p =
	for j in 1:9

	  prob = j * 0.1
	  d = rand(Binomial(9, prob), N);
	  h = fit(Histogram, d, -0.5:1:9.5)

	  p[j] = plot(xlim=(0,9), xticks=0:9)
	  for (i, w) in enumerate(h.weights)
		plot!(p[j], [i-1, i-1], [0.0, w], color=:blue,
		  leg=false, title="prob=$(round(prob, digits=1))")
	  end
	end
end

# ╔═╡ 98f93450-f50c-11ea-3767-e526c9109464
plot(p..., layout=(3,3))

# ╔═╡ 99019636-f50c-11ea-3ac3-4107ebb330c5
md"## End of Fig3.6t.jl"

# ╔═╡ Cell order:
# ╟─6e665478-f50c-11ea-1861-3f509f92ac5c
# ╠═98ed9a78-f50c-11ea-1d20-7bcaee689a99
# ╠═98edd6be-f50c-11ea-26c3-ff6ac7641060
# ╟─98ee52ec-f50c-11ea-22bc-c1075a60bfed
# ╠═98f67062-f50c-11ea-23f3-539ae4f7e67b
# ╠═98f93450-f50c-11ea-3767-e526c9109464
# ╟─99019636-f50c-11ea-3ac3-4107ebb330c5


text:
@testset "mask_area" begin
    for i in 1:5
        mask = rand(0:1, 156, 85)
        @test mask_area(MUMBAI_COORDINATE_SYSTEM, mask) >= 0
    end
    for i in 1:5
        mask = zeros(156, 85)
        @test mask_area(MUMBAI_COORDINATE_SYSTEM, mask) == 0
    end
end


text:
export default_type, default_library, similar_library, library

"""
    default_type(d::FullDim, ::Type{T}) where {T}

Returns the default polyhedron type for `d`-dimensional polyhedron of coefficient type `T`.
"""
function default_type end

function default_type(d::StaticArrays.Size{N}, ::Type{T}) where {N, T}
    return DefaultPolyhedron{T, Intersection{T, StaticArrays.SVector{N[1], T}, typeof(d)}, Hull{T, StaticArrays.SVector{N[1], T}, typeof(d)}}
end
function default_type(d::StaticArrays.Size{(1,)}, ::Type{T}) where T
    return Interval{T, StaticArrays.SVector{1, T}, typeof(d)}
end
function default_type(d::Int, ::Type{T}) where T
    if d == 1
        return Interval{T, Vector{T}, typeof(d)}
    else
        return DefaultPolyhedron{T, Intersection{T, Vector{T}, typeof(d)}, Hull{T, Vector{T}, typeof(d)}}
    end
end

"""
    default_library(d::FullDim, ::Type{T}) where {T}

Returns the default polyhedral library for `d`-dimensional polyhedron of coefficient type `T`.

## Examples

To obtain the default library for 2-dimensional polyhedra of eltype `Float64`,
do `default_library(2, Float64)`.

Given an `StaticArrays.SVector` `v`, to obtain a default library for points
of the type of `v` in a type stable way, do
`default_library(Polyhedra.FullDim(v), eltype(v))`.
"""
function default_library end

_default_type(::Type{T}) where T = T
# See https://github.com/JuliaPolyhedra/Polyhedra.jl/issues/35
_default_type(::Type{AbstractFloat}) = Float64
default_library(::StaticArrays.Size, T::Type) = DefaultLibrary{_default_type(T)}()
default_library(::StaticArrays.Size{(1,)}, T::Type) = IntervalLibrary{_default_type(T)}()
function default_library(d::Int, ::Type{T}) where T
    if d == 1
        return IntervalLibrary{_default_type(T)}()
    else
        return DefaultLibrary{_default_type(T)}()
    end
end

"""
    similar_library(lib::Library, d::FullDim, T::Type)

Returns a library that supports polyhedra of full dimension `T` with coefficient type `T`. If `lib` does not support it, this commonly calls `default_library(d, T)`.
"""
function similar_library end

# Shortcuts
similar_library(p::Polyhedron, d::FullDim, ::Type{T}) where T = similar_library(library(p), d, T)
similar_library(p::Polyhedron, ::Type{T}) where T = similar_library(library(p), T)
similar_library(p::Polyhedron{T}, d::FullDim) where T = similar_library(p, d, T)

"""
    library(p::Polyhedron)

Returns the library used by `p`.
"""
function library end

"""
    default_solver(p::Rep)

Returns a default linear programming solver for the polyhedron `p` (e.g. CDD has an internal solver which is used by default).
"""
default_solver(p::Rep) = nothing
function default_solver(p::Rep, ps::Rep...)
    s = default_solver(p)
    if s === nothing
        default_solver(ps...)
    else
        s
    end
end

"""
    solver(p::Rep, solver::MathProgBase.AbstractMathProgSolver=default_solver(p))

If the V-representation of `p` has been computed, returns `VRepSolver()`, otherwise, returns `solver`.
"""
function solver(p::Rep, solver::Union{Nothing, MPB.AbstractMathProgSolver}=default_solver(p))
    if vrepiscomputed(p)
        VRepSolver()
    else
        solver
    end
end

solver(v::VRepresentation, solver::Union{Nothing, MPB.AbstractMathProgSolver}=default_solver(v)) = VRepSolver()
solver(h::HRepresentation, solver::Union{Nothing, MPB.AbstractMathProgSolver}=default_solver(h)) = solver

_promote_reptype(P::Type{<:HRep}, ::Type{<:HRep}) = P
_promote_reptype(P::Type{<:VRep}, ::Type{<:VRep}) = P
# Breaks ambiguity with above two methods
_promote_reptype(P::Type{<:Polyhedron}, ::Type{<:Polyhedron}) = P

function promote_reptype(P1::Type{<:Rep}, P2::Type{<:Rep})
    _promote_reptype(P1, P2)
end
promote_reptype(P1::Type{<:Rep}, P2::Type{<:Rep}, P::Type{<:Rep}...) = promote_reptype(promote_reptype(P1, P2), P...)
promote_reptype(P::Type{<:Rep}) = P

# whether Rep has the constructor Rep(::It...; solver=...)
supportssolver(::Type{<:Rep}) = false

function constructpolyhedron(RepT::Type{<:Rep{T}}, d::FullDim, p::Tuple{Vararg{Rep}}, it::It{T}...) where T
    if supportssolver(RepT)
        solver = default_solver(p...)
        if solver !== nothing
            return RepT(d, it..., solver=solver)
        end
    end
    RepT(d, it...)::RepT # FIXME without this type annotation even convexhull(::PointsHull{2,Int64,Array{Int64,1}}, ::PointsHull{2,Int64,Array{Int64,1}}) is not type stable, why ?
end

function default_similar(p::Tuple{Vararg{Rep}}, d::FullDim, ::Type{T}, it::It{T}...) where T
    # Some types in p may not support `d` or `T` so we call `similar_type` after `promote_reptype`
    RepT = similar_type(promote_reptype(typeof.(p)...), d, T)
    constructpolyhedron(RepT, d, p, it...)
end

"""
    similar(p::Tuple{Vararg{Polyhedra.Rep}}, d::Polyhedra.FullDim, ::Type{T}, it::Polyhedra.It{T}...)

Creates a representation with a type similar to `p` of a polyhedron of full dimension `d`, element type `T` and initialize it with the iterators `it`.
The type of the result will be chosen closer to the type of `p[1]`.
"""
Base.similar(p::Tuple{Vararg{Rep}}, d::FullDim, ::Type{T}, it::It{T}...) where {T} = default_similar(p, d, T, it...)
function promote_coefficient_type(p::Tuple{Vararg{Rep}})
    promote_type(coefficient_type.(p)...)
end
function Base.similar(p::Tuple{Vararg{Rep}}, d::FullDim, it::It...)
    T = promote_coefficient_type(p)
    similar(p, d, T, it...)
end
function Base.similar(p::Tuple{Vararg{Rep}}, it::It...)
    similar(p, FullDim(p[1]), it...)
end
Base.similar(p::Rep, args...) = similar((p,), args...)


text:
# Struct of distribution, corresponding parameters, and a sample.
struct DistSpec{VF<:VariateForm,VS<:ValueSupport,F,T,X,G,B<:Tuple}
    name::Symbol
    f::F
    "Distribution parameters."
    θ::T
    "Sample."
    x::X
    "Transformation of sample `x`."
    xtrans::G
    "Broken backends"
    broken::B
end

function DistSpec(f, θ, x, xtrans=nothing; broken=())
    name = f isa Distribution ? nameof(typeof(f)) : nameof(typeof(f(θ...)))
    return DistSpec(name, f, θ, x, xtrans; broken=broken)
end

function DistSpec(name::Symbol, f, θ, x, xtrans=nothing; broken=())
    F = f isa Distribution ? typeof(f) : typeof(f(θ...))
    VF = Distributions.variate_form(F)
    VS = Distributions.value_support(F)
    return DistSpec{VF,VS,typeof(f),typeof(θ),typeof(x),typeof(xtrans),typeof(broken)}(
        name, f, θ, x, xtrans, broken,
    )
end

Distributions.variate_form(::Type{<:DistSpec{VF}}) where VF = VF
Distributions.value_support(::Type{<:DistSpec{VF,VS}}) where {VF,VS} = VS

# Auxiliary method for vectorizing parameters and samples
vectorize(v::Number) = [v]
vectorize(v::Diagonal) = v.diag
vectorize(v::AbstractVector{<:AbstractMatrix}) = mapreduce(vectorize, vcat, v)
vectorize(v) = vec(v)

"""
    unpack(x, inds, original...)

Return a tuple of unpacked parameters and samples in vector `x`.

Here `original` are the original full set of parameters and samples, and
`inds` contains the indices of the original parameters and samples for which
a possibly different value is given in `x`. If no value is provided in `x`,
the original value of the parameter is returned. The values are returned
in the same order as the original parameters.
"""
function unpack(x, inds, original...)
    offset = 0
    newvals = ntuple(length(original)) do i
        if i in inds
            v, offset = unpack_offset(x, offset, original[i])
        else
            v = original[i]
        end
        return v
    end
    offset == length(x) || throw(ArgumentError())

    return newvals
end

# Auxiliary methods for unpacking numbers and arrays
function unpack_offset(x, offset, original::Number)
    newoffset = offset + 1
    val = x[newoffset]
    return val, newoffset
end
function unpack_offset(x, offset, original::AbstractArray)
    newoffset = offset + length(original)
    val = reshape(x[(offset + 1):newoffset], size(original))
    return val, newoffset
end
function unpack_offset(x, offset, original::AbstractArray{<:AbstractArray})
    newoffset = offset
    val = map(original) do orig
        out, newoffset = unpack_offset(x, newoffset, orig)
        return out
    end
    return val, newoffset
end

# Run AD tests of a
function test_ad(dist::DistSpec; kwargs...)
    @info "Testing: $(dist.name)"

    f = dist.f
    θ = dist.θ
    x = dist.x
    g = dist.xtrans
    broken = dist.broken

    # Create functions with all possible arguments
    f_loglik_allargs = let f=f, g=g
        function (x, θ...)
            dist = f(θ...)
            xtilde = g === nothing ? x : g(x)
            return loglikelihood(dist, xtilde)
        end
    end
    f_logpdf_allargs = let f=f, g=g
        function (x, θ...)
            dist = f(θ...)
            xtilde = g === nothing ? x : g(x)
            if dist isa UnivariateDistribution && xtilde isa AbstractArray
                return sum(logpdf.(dist, xtilde))
            else
                return sum(logpdf(dist, xtilde))
            end
        end
    end

    # For all combinations of distribution parameters `θ`
    for inds in powerset(2:(length(θ) + 1))
        # Test only distribution parameters
        if !isempty(inds)
            xtest = mapreduce(vcat, inds) do i
                vectorize(θ[i - 1])
            end
            f_loglik_test = let xorig=x, θorig=θ, inds=inds
                x -> f_loglik_allargs(unpack(x, inds, xorig, θorig...)...)
            end
            f_logpdf_test = let xorig=x, θorig=θ, inds=inds
                x -> f_logpdf_allargs(unpack(x, inds, xorig, θorig...)...)
            end

            @test f_loglik_test(xtest) ≈ f_logpdf_test(xtest)

            test_ad(f_loglik_test, xtest, broken; kwargs...)
            test_ad(f_logpdf_test, xtest, broken; kwargs...)
        end

        # Test derivative with respect to location `x` as well
        # if the distribution is continuous
        if Distributions.value_support(typeof(dist)) === Continuous
            xtest = isempty(inds) ? vectorize(x) : vcat(vectorize(x), xtest)
            push!(inds, 1)
            f_loglik_test = let xorig=x, θorig=θ, inds=inds
                x -> f_loglik_allargs(unpack(x, inds, xorig, θorig...)...)
            end
            f_logpdf_test = let xorig=x, θorig=θ, inds=inds
                x -> f_logpdf_allargs(unpack(x, inds, xorig, θorig...)...)
            end

            @test f_loglik_test(xtest) ≈ f_logpdf_test(xtest)

            test_ad(f_loglik_test, xtest, broken; kwargs...)
            test_ad(f_logpdf_test, xtest, broken; kwargs...)
        end
    end
end

function test_ad(f, x, broken = (); rtol = 1e-6, atol = 1e-6)
    finitediff = FDM.grad(central_fdm(5, 1), f, x)[1]

    if AD == "All" || AD == "Tracker"
        if :Tracker in broken
            @test_broken Tracker.data(Tracker.gradient(f, x)[1]) ≈ finitediff rtol=rtol atol=atol
        else
            @test Tracker.data(Tracker.gradient(f, x)[1]) ≈ finitediff rtol=rtol atol=atol
        end
    end

    if AD == "All" || AD == "ForwardDiff"
        if :ForwardDiff in broken
            @test_broken ForwardDiff.gradient(f, x) ≈ finitediff rtol=rtol atol=atol
        else
            @test ForwardDiff.gradient(f, x) ≈ finitediff rtol=rtol atol=atol
        end
    end

    if AD == "All" || AD == "Zygote"
        if :Zygote in broken
            @test_broken Zygote.gradient(f, x)[1] ≈ finitediff rtol=rtol atol=atol
        else
            @test Zygote.gradient(f, x)[1] ≈ finitediff rtol=rtol atol=atol
        end
    end

    if AD == "All" || AD == "ReverseDiff"
        if :ReverseDiff in broken
            @test_broken ReverseDiff.gradient(f, x) ≈ finitediff rtol=rtol atol=atol
        else
            @test ReverseDiff.gradient(f, x) ≈ finitediff rtol=rtol atol=atol
        end
    end

    return
end


text:
# indexing tests
x = create_var(4, 6)
y = x[1:3, 5:end]
z = x[:]
t = [1, 2, 3]
f = t * x[2, 3]
row = x[1, :]
col = x[:, 1]
entry = x[4]
test_expr(y, (x,))
test_expr(z, (x,))
test_expr(f, (x,))
test_expr(row, (x,))
test_expr(col, (x,))
test_expr(entry, (x,))

info("All index tests passed")


text:
module Move_Pendulum

using  Modia3D
import Modia3D.ModiaMath

# Properties
vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)
vmat2 = Modia3D.Material(color="Red")
Lx = 1.6
Ly = 0.2*Lx
Lz = 0.2*Lx
m  = 1.0

# Pendulum
world  = Modia3D.Object3D(Modia3D.CoordinateSystem(0.5*Lx))
body   = Modia3D.Object3D(Modia3D.Solid(Modia3D.SolidBeam(Lx,Ly,Lz), Modia3D.MassProperties(m=m), vmat1))
frame1 = Modia3D.Object3D(body; r=[-Lx/2, 0.0, 0.0])
cyl    = Modia3D.Object3D(frame1,Modia3D.Cylinder(Ly/2,1.2*Ly; material=vmat2))
rev    = Modia3D.Revolute(world, frame1)

# Visualize pendulum
# Modia3D.visualizeWorld!(world, sceneOptions=Modia3D.SceneOptions(visualizeFrames=true, defaultFrameLength=0.3*Lx))

scene  = Modia3D.Scene(Modia3D.SceneOptions(visualizeFrames=true, defaultFrameLength=0.3*Lx))
Modia3D.initAnalysis!(world, scene)

tStart = 0.0
tEnd   = 3.0

for time = range(tStart, stop=tEnd, length=101)
  delta_phi = Modia3D.linearMovement(2*pi, tStart, tEnd, time)
  Modia3D.setAngle!(rev, delta_phi)
  Modia3D.updatePosition!(world)
  Modia3D.visualize!(scene,time)
end
Modia3D.closeAnalysis!(scene)

println("... success of Move_Pendulum.jl!")
end


text:
module GenStdLib

import Base.Docs: Binding, DocStr

# Constants.
const DOCSTRING_DIRECTIVE = r"^(.. (function|type|data):: ).*"

# Types.

typealias Signature Tuple{Binding, Type}

type State
    files      :: Dict{String, Vector{String}}
    validdocs  :: Dict{String, Tuple{Module, Signature, DocStr}}
    baddocs    :: Dict{Signature, Tuple{Module, DocStr}}
    documented :: Dict{Signature, Tuple{Module, DocStr, String}}
    errorlevel :: Int
    debug      :: Bool
    State() = new(Dict(), Dict(), Dict(), Dict(), 0, "JULIA_GENSTDLIB_DEBUG" in keys(ENV))
end

# Documentation Translator.

function translate(dirs::Vector; root = dirname(@__FILE__))
    cd(root) do
        state = State()
        println("\n# GENSTDLIB\n")
        # Find all available docstrings within `Base` and its submodules.
        info("loading docstrings from modules.")
        loaddocs!(state)
        # Read in rst pages and update all docstrings.
        info("parsing external documentation.")
        for dir in dirs, file in readdir(dir)
            translate(state, joinpath(dir, file))
        end
        # Write the newly updated rst docs back to their files -- replaces original content.
        if state.errorlevel < 2
            info("writing documentation back to file.")
            for (file, contents) in state.files
                open(file, "w") do io
                    for line in contents
                        println(io, line)
                    end
                end
            end
        else
            warn("errors found while generating documentation. Aborting file writing.")
        end
        # Report statistics about 'Base' docs.
        summarise(state)
    end
end
function translate(state::State, file::AbstractString)
    input, output = split(readstring(file), '\n')[1:end-1], []
    while !isempty(input)
        if ismatch(DOCSTRING_DIRECTIVE, first(input))
            append!(output, getdoc(state, file, input))
        else
            push!(output, shift!(input))
        end
    end
    state.files[file] = output
    return state
end

# Documentation Summary Report.

function summarise(state::State)
    println("\n# SUMMARY\n")
    println("   ", length(state.validdocs),  " docstrings found in modules.")
    println("   ", length(state.documented), " docstrings found in external docs.")

    # Of the missing docstrings, which are public and which are private?
    public, private = [], []
    for (sig, (mod, signature, docstr)) in state.validdocs
        # The docstring is valid, but not found in the external docs.
        if !haskey(state.documented, signature)
            binding, typesig = signature
            exported = binding.var in names(binding.mod)
            push!(exported ? public : private, (binding, typesig))
        end
    end
    println("   ", length(public), " public docstrings missing from external docs.")
    state.debug && (println(); foreach(printmissing, public); println())
    println("   ", length(private), " private docstrings missing from external docs.")
    state.debug && (println(); foreach(printmissing, private); println())

    # Expected format is: first element in docstring should be a code block.
    println("   ", length(state.baddocs),    " docstrings not matching expected format. Skipped.")
    println("   ", length(state.files),      " files parsed.\n")
    state.debug || println("Set 'JULIA_GENSTDLIB_DEBUG' ENV for additional debug info.\n")
end
printmissing(x) = println("     X ", x[1], " :: ", x[2])

# Documentation Loader.

function loaddocs!(state::State)
    for mod in Base.Docs.modules
        for (binding, multidoc) in Base.Docs.meta(mod)
            for (typesig, docstr) in multidoc.docs
                loaddocs!(state, mod, binding, typesig, docstr)
            end
        end
    end
    for (keyword, docstr) in Base.Docs.keywords
        loaddocs!(state, Main, keyword, Union{}, docstr)
    end
    return state
end
function loaddocs!(state::State, mod, binding::Binding, typesig, docstr)
    markdown = Base.Docs.parsedoc(docstr)
    if validdocstr(markdown)
        code = rstrip(markdown.content[1].code)
        if haskey(state.validdocs, code) && state.validdocs[code][3] !== docstr
            code = indent(code)
            warn("duplicate signature found for '$binding' in module '$mod':\n\n$code\n")
            state.errorlevel = 2
        else
            state.validdocs[code] = (mod, (binding, typesig), docstr)
        end
    else
        if haskey(state.baddocs, (binding, typesig))
            warn("duplicate binding '$binding :: $typesig' found in module '$mod'.")
            state.errorlevel = 2
        else
            state.baddocs[(binding, typesig)] = (mod, docstr)
        end
    end
    return state
end
function loaddocs!(state::State, mod, keyword::Symbol, typesig, docstr)
    binding = Base.Docs.Binding(mod, keyword)
    loaddocs!(state, mod, binding, typesig, docstr)
end

# Retrieve and format docstrings.

function getdoc(state::State, file::AbstractString, input::Vector)
    # Capture the lines containing the docstring signature.
    output = []
    while !isempty(input)
        line = rstrip(shift!(input))
        push!(output, line)
        ismatch(r"^$", line) && break
    end
    # Recover the unindented version of the signature.
    n = length(match(DOCSTRING_DIRECTIVE, first(output))[1])
    b = IOBuffer()
    for line in output[1:end-1]
        println(b, line[(n + 1):end])
    end
    # The signature may contain `\` characaters, which must be unescaped.
    signature = unescape_string(rstrip(String(take!(b))))
    # Splice the correct docstring into the output after the signature.
    if haskey(state.validdocs, signature)
        # Push the rst text for the docstring into the output.
        mod, (binding, typesig), docstr = state.validdocs[signature]
        md  = Markdown.MD(Base.Docs.parsedoc(docstr).content[2:end])
        rst = Base.Markdown.rst(dropheaders(md))
        push!(output, "   .. Docstring generated from Julia source", "")
        for line in split(rst, '\n')
            line = isempty(line) ? "" : string(" "^3, line)
            push!(output, rstrip(line))
        end
        # Consume all indented lines from the current docstring.
        while !isempty(input)
            line = first(input)
            ismatch(r"^[^\s]", line) ? break : shift!(input)
        end
        # Track which docstrings have been found in the external docs.
        state.documented[(binding, typesig)] = (mod, docstr, signature)
    else
        signature = indent(signature)
        warn("missing docs for signature:\n\n$signature\n")
        state.errorlevel = 1
    end
    return output
end

# Replace headers in docs with bold since Sphinx does not allow headers inside docstrings.
dropheaders(md) = Markdown.MD(map(bold, md.content))
bold(x::Markdown.Header) = Markdown.Paragraph(Markdown.Bold(x.text))
bold(other) = other

# Utilities.

function indent(str::AbstractString, indent = 4)
    buf = IOBuffer()
    for line in split(str, '\n')
        println(buf, " "^indent, line)
    end
    String(take!(buf))
end

function validdocstr(markdown::Base.Markdown.MD)
    content = markdown.content
    !isempty(content) && isa(first(content), Base.Markdown.Code)
end
validdocstr(other) = false

end

# The docstring for `build_sysimg` is defined is this file and included within the
# `devdocs/sysimg.rst` file, so we include it here to make it visible to the docsystem.
include(joinpath("..", "contrib", "build_sysimg.jl"))

GenStdLib.translate(["manual", "stdlib", "devdocs"])

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
using Test
import Logging
import Pluto: Configuration, Notebook, ServerSession, ClientSession, update_run!, Cell, WorkspaceManager, save_notebook

@testset "Frontmatter" begin
    d = Dict{String,Any}("title" => 123, "wow" => ["a", "b"])
    e = Dict{String,Any}()
    cases = [
        1 => (
            Notebook([
                Cell("""
                frontmatter = Dict{Any,Any}("title" => 123, "wow" => ["a", "b"])
                """),
            ]),
            d,
        ),
        2 => (
            Notebook([
                Cell("""
                frontmatter = (title = 123, wow = ["a", "b"])
                """),
            ]),
            d
        ),
        3 => (
            Notebook([
                Cell("""
                frontmatter = 123
                """),
            ]),
            e
        ),
        4 => (
            Notebook([
                Cell("""
                frontmatter = 123
                """),
                Cell("""
                frontmatter = 123
                """),
            ]),
            e
        ),
        5 => (
            Notebook([
                Cell("""
                """),
            ]),
            e
        ),
    ]

    @testset "Notebook $i" for (i, (nb, expected)) in cases
        save_notebook(nb)

        Logging.with_logger(Logging.NullLogger()) do
            @test Pluto.frontmatter(nb.path) == expected
            @test Pluto.frontmatter(nb) == expected
            @test typeof(Pluto.frontmatter(nb)) === Dict{String,Any}
        end

        if i ∈ (3,4)
            @test_throws Exception Pluto.frontmatter(nb.path, raise=true)
        end

    end
end


text:
## serialization helper type: Packed31 ##

"""
Helper type to read/write 31 bits in 1/2/5 bytes.

Why 31 bits and not 32?

 * encoding is "dense", no redundancy
 * category+length in FlyToken have exactly 31 bits

Properties defined:

 * .nibble: lowest 4 bits as a Nibble

 * .all: whole value as UInt32

 * .len: Bits 4..30  as UInt32

"""
primitive type Packed31 32 end

function Packed31(v::Union{UInt32,Int32})
    reinterpret(Packed31,v)
end

Base.UInt32(v::Packed31) = reinterpret(UInt32,v)

Base.show(io::IO, v::Packed31) = show(io,UInt32(v))

"extracts the lowest 4 bits from  an unsigned value as an UInt8"
function bits0_3(v)
    (v%UInt8) & 0x0f
end
bits0_3(v::Packed31) = bits0_3(UInt32(v))

nibble(v::Packed31) = Nibble(bits0_3(v))

"right shift by a nibble (4 bits)"
bits4_30(v::Packed31) = UInt32(v) >> 4

function Packed31(bits0_3::UInt8, bits4_30::UInt32)
    @boundscheck checkulimit(bits0_3,15)
    @boundscheck checkulimit(bits4_30,(1<<27)-1)
    Packed31((bits4_30<<4) + bits0_3)
end

function Packed31(bits0_3::Nibble, bits4_30::UInt32)
    @boundscheck checkulimit(bits4_30,(1<<27)-1)
    Packed31((bits4_30<<4) + UInt8(bits0_3))
end

function Base.getproperty(p::Packed31, s::Symbol)
    if s === :nibble
        return Nibble(bits0_3(p))
    elseif s === :all
        return UInt32(p)
    elseif s === :len
        return bits4_30(p)%UInt32
    else
        return Core.getproperty(p, s)
    end
end

"""
    read(io::IO,::Packed31)

read 31 bits in a variable length encoding format of 1/2/5 bytes.
"""
function Base.read(io::IO,::Type{Packed31})
    v = read(io,UInt8) %UInt32
    if v <= 127
        return Packed31(v)
    end
    v = ((v-128) <<8 ) | (read(io,UInt8)%UInt32)
    if v <= 127
        # 5 byte encoding
        v = (value << 24) |  ((read(io,UInt8)%UInt32)<<16)  |  ((read(io,UInt8)%UInt32)<<8)  |  ((read(io,UInt8)%UInt32))
    end
    return Packed31(v)
end

"""
    write(io::IO, v::Packed31)

write 31 bits in variable length encoding in an platform independent format.
0..127 are encoded in 1 byte, 128..32767 in 2 bytes, all other need 5 bytes.

This encoding is recommended for serialization of nonnegative Int32
values if encoded size matters and if small values dominate,
e.g. 50% are below 32768 or 25% are below 128.

"""
function Base.write(io::IO, v::Packed31)
    u = UInt32(v)
    if u<=127
        write(io,u%UInt8)
        return nothing
    end
    if  u <= 32767
        write(io,hton(u%UInt16+0x8000))
        return nothing
    end
    write(io,0x80)
    write(io,hton(u))
    nothing
end


text:
using Example2018
using Test

# Polynomials for a ring:
# 1. addition is associative (done)
# 2. additive identity (done)
# 3. additive inverse (done)
# 4. addition commutes
# 5. multiplication.. (not implemented yet)

@testset "Polynomials form a ring" begin

z=Polynomial([0,0,0])

for n in 1:100
    deg=rand(0:100)
    p=Polynomial(rand(-1000000:1000000,deg))
    q=-p

    @test p+z==p # 2.
    @test z+p==p # 2.
    @test p+q==z # 3.

    dega=rand(0:100); degb=rand(0:100); degc=rand(0:100)
    a=Polynomial(rand(-1000000:1000000,dega))
    b=Polynomial(rand(-1000000:1000000,degb))
    c=Polynomial(rand(-1000000:1000000,degc))
    prov1=a+b
    prov2=b+c
    @test prov1+c==a+prov2 # 1.
    @test a+b==b+a #4.

end

    end


text:
#= ######################################################
Exemplary models
###################################################### =#

function model_two_intesecting_circles(d, seed, leftout)
    # two intersecting circles
    m = create_empty_model()
    center1 = Point2D(4,3)
    add!(m, circle_open(d, center1, seed = seed, leftout = leftout, open = "right", dir = "pos"))
    dist_x = m.nodes[end].x - center1.x
    center2 = Point2D(center1.x + 2 * dist_x, center1.y)
    add!(m, circle_open(d, center2, seed = seed, leftout = leftout, open = "left", dir = "pos"))
    offset_model!(m)
    mi = make_model_immutable(m)
    println("model has ", mi.nelem, " elements")
    return mi
end

function model_rectangle(a, b, elemsize)
    # single rectangle
    m = create_empty_model()
    seeda = round(Integer,a/elemsize)
    seedb = round(Integer,b/elemsize)
    center = Point2D(0.0,0.0)
    add!(m, edge(center, center+Point2D(a,0.0), seed = seeda, dir = "pos"))
    add!(m, edge(center+Point2D(a,0.0), center+Point2D(a,b), seed = seedb, dir = "pos"))
    add!(m, edge(center+Point2D(a,b), center+Point2D(0.0,b), seed = seeda, dir = "pos"))
    add!(m, edge(center+Point2D(0.0,b), center, seed = seedb, dir = "pos"))
    offset_model!(m)
    mi = make_model_immutable(m)
    println("model has ", mi.nelem, " elements")
    return mi
end

function model_circle_in_circle(d_in, d_out, elemsize)
    # circle in circle
    m = create_empty_model()
    seed_in = round(Integer,(pi*d_in)/elemsize)
    seed_out = round(Integer,(pi*d_out)/elemsize)
    center = Point2D(0.0,0.0)
    add!(m, circle(d_in, center, seed = seed_in, dir = "neg"))
    add!(m, circle(d_out, center, seed = seed_out, dir = "pos"))
    offset_model!(m)
    mi = make_model_immutable(m)
    println("model has ", mi.nelem, " elements")
    return mi
end

function model_circle_with_opening_line(d, seed, leftout)
    # open circle closed with line
    m = create_empty_model()
    center = Point2D(0.0,0.0)
    add!(m, circle_open(d, center, seed = seed, leftout = leftout, open = "right", dir = "pos"))
    p1 = m.nodes[1]
    p2 = m.nodes[end]
    seed_edge = leftout + 1
    add!(m, edge(p1, p2, seed = seed_edge, dir = "neg"))
    offset_model!(m)
    mi = make_model_immutable(m)
    println("model has ", mi.nelem, " elements")
    return mi
end

function model_isosceles_triangle(a, alpha_deg, elemsize)
    # isocles triangle
    m = create_empty_model()
    alpha_rad = alpha_deg*pi/180
    point = Point2D(sin(alpha_rad) * a, cos(alpha_rad) * a)
    l = sqrt((point.x-0)^2 + (point.y-a)^2)
    seeda = round(Integer,a/elemsize)
    seedl = round(Integer,l/elemsize)
    center = Point2D(0.0,0.0)
    add!(m, edge(center, center+Point2D(0.0,a), seed = seeda, dir = "neg"))
    add!(m, edge(center, point, seed = seeda, dir = "pos"))
    add!(m, edge(center+Point2D(0.0,a), point, seed = seedl, dir = "neg"))
    offset_model!(m)
    mi = make_model_immutable(m)
    println("model has ", mi.nelem, " elements")
    return mi
end

function model_two_rectangles_with_holes(elemsize)
    # two rectangles connected with a pinhole
    m = create_empty_model()
    l_spalt = 0.1
    x_offset = 0.01
    # left rect with hole
    a = 0.6
    b = 0.6
    seeda = round(Integer,a/elemsize)
    seedb = round(Integer,b/elemsize)
    seedb2 = round(Integer,((b-l_spalt)/2)/elemsize)
    center = Point2D(0.0,0.0)
    add!(m, edge(center, center+Point2D(a,0.0), seed = seeda, dir = "pos"))
    add!(m, edge(center+Point2D(a,b), center+Point2D(0.0,b), seed = seeda, dir = "pos"))
    add!(m, edge(center+Point2D(0.0,b), center, seed = seedb, dir = "pos"))
    add!(m, edge(center+Point2D(a,0.0), center+Point2D(a,(b-l_spalt)/2), seed = seedb2, dir = "pos"))
    add!(m, edge(center+Point2D(a,(b+l_spalt)/2), center+Point2D(a,b), seed = seedb2, dir = "pos"))
    # right rect with hole
    c = 0.6
    d = 1.0
    y_offset = (b-d)/2
    seedc = round(Integer,c/elemsize)
    seedd = round(Integer,d/elemsize)
    seedd2 = round(Integer,((d-l_spalt)/2)/elemsize)
    center = Point2D(a+x_offset,y_offset)
    add!(m, edge(center, center+Point2D(c,0.0), seed = seedc, dir = "pos"))
    add!(m, edge(center+Point2D(c,d), center+Point2D(0.0,d), seed = seedc, dir = "pos"))
    add!(m, edge(center+Point2D(c,0.0), center+Point2D(c,d), seed = seedd, dir = "pos"))
    add!(m, edge(center+Point2D(0.0,0.0), center+Point2D(0.0,(d-l_spalt)/2), seed = seedd2, dir = "neg"))
    add!(m, edge(center+Point2D(0.0,(d+l_spalt)/2), center+Point2D(0.0,d), seed = seedd2, dir = "neg"))
    offset_model!(m)
    mi = make_model_immutable(m)
    println("model has ", mi.nelem, " elements")
    return mi
end

function model_lines()
    # two opposing lines
    m = create_empty_model()
    seed = 10
    a = 3.0
    d = 0.5
    add!(m, edge(Point2D(0.0,0.0), Point2D(a,0.0), seed = seed, dir = "pos"))
    add!(m, edge(Point2D(0.0,d), Point2D(a,d), seed = seed, dir = "neg"))
    offset_model!(m)
    mi = make_model_immutable(m)
    println("model has ", mi.nelem, " elements")
    return mi
end

function model_cosinus(a, b, seed)
    # cosinus
    m = create_empty_model()
    add!(m, cosinus(a, b, Point2D(0.0,0.0), seed = seed, dir = "pos"))
    elemsize = m.elem[1].area
    l_edge = get_length(m.nodes[1], m.nodes[end])
    seed2 = round(Integer,(l_edge/elemsize))
    add!(m, edge(m.nodes[1], m.nodes[end], seed = seed2, dir = "neg"))
    offset_model!(m)
    mi = make_model_immutable(m)
    println("model has ", mi.nelem, " elements")
    return mi
end

function model_trapez(a, b, h, elemsize)
    # create model for trapez
    p1 = Point2D((b-a)/2, 0.0)
    p2 = Point2D((b-a)/2 + a, 0.0)
    p3 = Point2D(0.0, h)
    p4 = Point2D(b, h)
    l = get_length(p3, p1)
    m = create_empty_model()
    add!(m, edge(p3, p1, seed = round(Integer,l/elemsize), dir = "pos"))
    add!(m, edge(p1, p2, seed = round(Integer,a/elemsize), dir = "pos"))
    add!(m, edge(p2, p4, seed = round(Integer,l/elemsize), dir = "pos"))
    add!(m, edge(p4, p3, seed = round(Integer,b/elemsize), dir = "pos"))
    offset_model!(m)
    mi = make_model_immutable(m)
    println("model has ", mi.nelem, " elements")
    return mi
end

function model_right_triangle(a, b, elemsize)
    # create model
    center = Point2D(0.0, 0.0)
    m = create_empty_model()
    add!(m, edge(center, center+Point2D(0.0,a), seed = round(Integer,a/elemsize), dir = "neg"))
    add!(m, edge(center, center+Point2D(b,0.0), seed = round(Integer,b/elemsize), dir = "pos"))
    l = get_length(Point2D(0.0,a), Point2D(b,0.0))
    add!(m, edge(center+Point2D(0.0,a), center+Point2D(b,0.0), seed = round(Integer,l/elemsize), dir = "neg"))
    offset_model!(m)
    mi = make_model_immutable(m)
    println("model has ", mi.nelem, " elements")
    return mi
end

function model_furnace1()
    # furnace after WS shr build info broschure
    m = create_empty_model()
    # user input
    d = 0.195 # d of shr
    hg = 0.06 # height of gut
    wg = 1.6 # width of gut
    hb = 0.05 # distance gut to furnace bottom
    # shr
    center_shr = Point2D(0.0,(0.5*d + 0.5*d + hg + hb))
    elemsize_shr = 0.02
    seedshr = round(Integer,(pi*d)/elemsize_shr)
    x = d + 0.5*d
    add!(m, circle(d, center_shr+Point2D(x,0.0), seed = seedshr, dir = "neg"))
    x = d + d + 2*d + 0.5*d
    add!(m, circle(d, center_shr+Point2D(x,0.0), seed = seedshr, dir = "neg"))
    x = d + d + 2*d + d + 2*d + 0.5*d
    add!(m, circle(d, center_shr+Point2D(x,0.0), seed = seedshr, dir = "neg"))
    # gut
    elemsize_g = 0.03
    seedwg = round(Integer,wg/elemsize_g)
    seedhg = round(Integer,hg/elemsize_g)
    center_g = Point2D(((d + d + 2*d + 0.5*d) - (0.5*wg)),hb)
    add!(m, edge(center_g, center_g+Point2D(wg,0.0), seed = seedwg, dir = "neg"))
    add!(m, edge(center_g+Point2D(wg,0.0), center_g+Point2D(wg,hg), seed = seedhg, dir = "neg"))
    add!(m, edge(center_g+Point2D(wg,hg), center_g+Point2D(0.0,hg), seed = seedwg, dir = "neg"))
    add!(m, edge(center_g+Point2D(0.0,hg), center_g, seed = seedhg, dir = "neg"))
    # furnace chamber
    w = d + d + 2*d + d + 2*d + d + d
    h = 0.5*d + d + 0.5*d + hg + hb
    elemsize_f = 0.05
    seedw = round(Integer,w/elemsize_f)
    seedh = round(Integer,h/elemsize_f)
    center_f = Point2D(0.0,0.0)
    add!(m, edge(center_f, center_f+Point2D(w,0.0), seed = seedw, dir = "pos"))
    add!(m, edge(center_f+Point2D(w,0.0), center_f+Point2D(w,h), seed = seedh, dir = "pos"))
    add!(m, edge(center_f+Point2D(w,h), center_f+Point2D(0.0,h), seed = seedw, dir = "pos"))
    add!(m, edge(center_f+Point2D(0.0,h), center_f, seed = seedh, dir = "pos"))
    # finish
    offset_model!(m)
    mi = make_model_immutable(m)
    println("model has ", mi.nelem, " elements")
    return mi
end

function model_furnace2(d, t, n, h)
    # simplified furnace for specht book shr study
    m = create_empty_model()
    # user input
    # d = 0.195 # d of shr
    # t = 0.5 # distance between centers of shr
    # n = 3 # number of shr
    # h = 0.5 # height of furnace
    w = n*t # width of furnace
    # shr
    center_shr = Point2D(0.0,h/2)
    elemsize_shr = 0.02
    seedshr = round(Integer,(pi*d)/elemsize_shr)
    for i = 1:n
        x = 0.5*t + (i-1)*t
        add!(m, circle(d, center_shr+Point2D(x,0.0), seed = seedshr, dir = "neg"))
    end
    # furnace chamber
    elemsize_f = 0.05
    seedw = round(Integer,w/elemsize_f)
    seedh = round(Integer,h/elemsize_f)
    center_f = Point2D(0.0,0.0)
    add!(m, edge(center_f+Point2D(w,0.0), center_f+Point2D(w,h), seed = seedh, dir = "pos"))
    add!(m, edge(center_f+Point2D(w,h), center_f+Point2D(0.0,h), seed = seedw, dir = "pos"))
    add!(m, edge(center_f+Point2D(0.0,h), center_f, seed = seedh, dir = "pos"))
    # gut
    add!(m, edge(center_f, center_f+Point2D(w,0.0), seed = seedw, dir = "pos"))
    # finish
    offset_model!(m)
    mi = make_model_immutable(m)
    println("model has ", mi.nelem, " elements")
    return mi
end

function model_labyrinth(path, d, elemsize)
    # labyrinth constructed with a path
    m = create_empty_model()
    # elemsize = 0.02
    # initial path
    # path = [0.4, -0.4, 0.4, 0.2, 0.5, -0.3, -0.2]
    # d = 0.1
    # offset path
    path_o = zeros(size(path,1))
    for i = 1:(size(path,1))
        if i == 1 || i == size(path,1)
            if i == 1
                if path[i+1] > 0
                    path_o[i] = path[i] + d
                else
                    path_o[i] = path[i] - d
                end
            else
                if path[i-1] > 0
                    if path[i] > 0
                        path_o[i] = path[i] + d
                    else
                        path_o[i] = path[i] + d
                    end
                else
                    if path[i] > 0
                        path_o[i] = path[i] + d
                    else
                        path_o[i] = path[i] + d
                    end
                end
            end
        else
            if path[i-1] > 0
                if path[i+1] > 0
                    path_o[i] = path[i]
                else
                    if path[i] > 0
                        path_o[i] = path[i] - 2*d
                    else
                        path_o[i] = path[i] + 2*d
                    end
                end
            else
                if path[i+1] < 0
                    path_o[i] = path[i]
                else
                    if path[i] > 0
                        path_o[i] = path[i] + 2*d
                    else
                        path_o[i] = path[i] - 2*d
                    end
                end
            end
        end
    end
    # start
    point1s = Point2D(0.0,0.0)
    point2s = Point2D(0.0,d)
    add!(m, edge(point1s, point2s, seed = round(Integer,d/elemsize), dir = "neg"))
    # path
    point1 = point1s
    point2 = point2s
    for i = 1:size(path,1)
        point1 = point2
        if isodd(i)
            # move in x
            point2 = Point2D(point1[1]+path[i],point1[2])
        else
            # move in y
            point2 = Point2D(point1[1],point1[2]+path[i])
        end
        # println(point1, " -> ",point2)
        l = point2 - point1
        seed = round(Integer,norm(l)/elemsize)
        add!(m, edge(point1, point2, seed = seed, dir = "neg"))
    end
    point1f = point2
    # path_o
    point1 = point2s
    point2 = point1s
    for i = 1:size(path_o,1)
        point1 = point2
        if isodd(i)
            # move in x
            point2 = Point2D(point1[1]+path_o[i],point1[2])
        else
            # move in y
            point2 = Point2D(point1[1],point1[2]+path_o[i])
        end
        # println(point1, " -> ",point2)
        l = point2 - point1
        seed = round(Integer,norm(l)/elemsize)
        add!(m, edge(point1, point2, seed = seed, dir = "pos"))
    end
    point2f = point2
    # finish
    add!(m, edge(point1f, point2f, seed = round(Integer,d/elemsize), dir = "neg"))
    offset_model!(m)
    mi = make_model_immutable(m)
    println("model has ", mi.nelem, " elements")
    return mi
end

function model_square_in_square(a, b, elemsize)
    # square in square
    # all edges are individual parts
    m = create_empty_model()
    seeda = round(Integer,a/elemsize)
    seedb = round(Integer,b/elemsize)
    # inner square
    center = Point2D(-0.5*a,-0.5*a)
    add!(m, edge(center, center+Point2D(a,0.0), seed = seeda, dir = "neg"))
    add!(m, edge(center+Point2D(a,0.0), center+Point2D(a,a), seed = seeda, dir = "neg"))
    add!(m, edge(center+Point2D(a,a), center+Point2D(0.0,a), seed = seeda, dir = "neg"))
    add!(m, edge(center+Point2D(0.0,a), center, seed = seeda, dir = "neg"))
    # outer square
    center = Point2D(-0.5*b,-0.5*b)
    add!(m, edge(center, center+Point2D(b,0.0), seed = seedb, dir = "pos"))
    add!(m, edge(center+Point2D(b,0.0), center+Point2D(b,b), seed = seedb, dir = "pos"))
    add!(m, edge(center+Point2D(b,b), center+Point2D(0.0,b), seed = seedb, dir = "pos"))
    add!(m, edge(center+Point2D(0.0,b), center, seed = seedb, dir = "pos"))
    offset_model!(m)
    mi = make_model_immutable(m)
    println("model has ", mi.nelem, " elements")
    return mi
end

function model_rect_in_rect(r1x, r1y, r2x, r2y, elemsize)
    # rectangle in rectangle
    m = create_empty_model()
    center = Point2D(0.0,0.0)
    seedx = round(Integer,r1x/elemsize)
    seedy = round(Integer,r1y/elemsize)
    add!(m, rectangle(r1x, r1y, center, seedx = seedx, seedy = seedy, dir = "pos"))
    seedx = round(Integer,r2x/elemsize)
    seedy = round(Integer,r2y/elemsize)
    add!(m, rectangle(r2x, r2y, center, seedx = seedx, seedy = seedy, dir = "neg"))
    offset_model!(m)
    mi = make_model_immutable(m)
    println("model has ", mi.nelem, " elements")
    return mi
end

function model_circles_in_circle_cross(d_in, d_out, n_diag, elemsize)
    # circles in circle - distribution: cross
    if iseven(n_diag)
        error("n_diag is an even number which is not allowed")
    end
    m = create_empty_model()
    seed_in = round(Integer,(pi*d_in)/elemsize)
    seed_out = round(Integer,(pi*d_out)/elemsize)
    center = Point2D(0.0,0.0)
    add!(m, circle(d_in, center, seed = seed_in, dir = "neg"))
    dist = (d_out/2) / ((n_diag-1)/2 + 1)
    distadd = dist
    for i = 1:((n_diag-1)/2)
        pos = Point2D(center.x+distadd,center.y)
        add!(m, circle(d_in, pos, seed = seed_in, dir = "neg"))
        pos = Point2D(center.x,center.y+distadd)
        add!(m, circle(d_in, pos, seed = seed_in, dir = "neg"))
        pos = Point2D(center.x-distadd,center.y)
        add!(m, circle(d_in, pos, seed = seed_in, dir = "neg"))
        pos = Point2D(center.x,center.y-distadd)
        add!(m, circle(d_in, pos, seed = seed_in, dir = "neg"))
        distadd += dist
    end
    add!(m, circle(d_out, center, seed = seed_out, dir = "pos"))
    offset_model!(m)
    mi = make_model_immutable(m)
    println("model has ", mi.nelem, " elements")
    return mi
end

function get_rand_pos_cirlce(d_in, d_out, center, n; sector = 1.0)
    # get random positions in circle
    # circlepart: 1.0 -> full circle 360 deg
    pos = Vector{Point2D{Float64}}(undef,n)
    desired_offset = d_in / 2 # user given distance between inner circles
    r = (d_out / 2) - (d_in / 2) - desired_offset
    i = 1
    while i <= n
        valid = true
        # println("i: ", i)
        angle = rand()*2*pi*sector
        dist = rand()*r
        # println("angle: ", angle, "with dist: ", dist)
        pos_t = Point2D(center.x+dist*cos(angle), center.y+dist*sin(angle))
        for j = 1:i
            pos_c = pos[j]
            dist_check = sqrt((pos_t.x-pos_c.x)^2 + (pos_t.y-pos_c.y)^2)
            if dist_check < (d_in + desired_offset)
                valid = false
                # println("penetration detected")
                break
            end
        end
        if valid == true
            # println("all clear")
            pos[i] = pos_t
            i += 1
        end
    end
    return pos
end

function model_circles_in_circle_rand(d_in, d_out, n, elemsize)
    # circles in circle - distribution: random
    m = create_empty_model()
    seed_in = round(Integer,(pi*d_in)/elemsize)
    seed_out = round(Integer,(pi*d_out)/elemsize)
    center = Point2D(0.0,0.0)
    pos = get_rand_pos_cirlce(d_in, d_out, center, n, sector = 1.0)
    for i = 1:n
        add!(m, circle(d_in, pos[i], seed = seed_in, dir = "neg"))
    end
    add!(m, circle(d_out, center, seed = seed_out, dir = "pos"))
    offset_model!(m)
    mi = make_model_immutable(m)
    println("model has ", mi.nelem, " elements")
    return mi
end

function model_circles_in_circle_rand_full(d_in, d_out, elemsize)
    # circles in circle - distribution: random case full
    m = create_empty_model()
    seed_in = round(Integer,(pi*d_in)/elemsize)
    seed_out = round(Integer,(pi*d_out)/elemsize)
    center = Point2D(0.0,0.0)
    pos = [ Point2D(-0.17160541430982149, 0.2336985360150676),
            Point2D(0.22451500551137007, -0.30768142248166275),
            Point2D(0.16828228357797986, -0.11171991410768657),
            Point2D(-0.17934496561024932, -0.7285037214126262),
            Point2D(-0.09085754774756338, 0.6457715735143245),
            Point2D(0.007325824302822567, 0.4179066847097554),
            Point2D(0.07716402068329058, -0.5337173965768568),
            Point2D(-0.2788863155366298, 0.398449173522567),
            Point2D(0.27973322582223453, -0.5352920087308007),
            Point2D(0.45515084453294763, 0.29132939345995706),
            Point2D(-0.12687522015645525, -0.25358199931056635),
            Point2D(-0.5898782348762814, 0.13012962748704132),
            Point2D(0.206434732538964, 0.4883359733387012),
            Point2D(0.25250771742586314, -0.7581101501890426),
            Point2D(0.05387174232420384, 0.15917369360119593),
            Point2D(0.41979923540554165, 0.6371896374415813),
            Point2D(-0.7000331990175072, -0.18264038174416705)]
    for i = 1:size(pos,1)
        add!(m, circle(d_in, pos[i], seed = seed_in, dir = "neg"))
    end
    add!(m, circle(d_out, center, seed = seed_out, dir = "pos"))
    offset_model!(m)
    mi = make_model_immutable(m)
    println("model has ", mi.nelem, " elements")
    return mi
end

function model_circles_in_circle_rand_half(d_in, d_out, elemsize)
    # circles in circle - distribution: random case half
    m = create_empty_model()
    seed_in = round(Integer,(pi*d_in)/elemsize)
    seed_out = round(Integer,(pi*d_out)/elemsize)
    center = Point2D(0.0,0.0)
    pos = [ Point2D(-0.05106910757240293, 0.5551081147435735),
            Point2D(0.6912293507261287, 0.1284491037094026),
            Point2D(0.3711110161319347, 0.1783226198209847),
            Point2D(-0.6655780206552504, 0.20649157369483673),
            Point2D(-0.14552889293587967, 0.1281009742606726),
            Point2D(-0.4507845542159105, 0.19789539410300944),
            Point2D(0.2400041066229391, 0.2845291942939924),
            Point2D(-0.05825139590019937, 0.32356766055237207),
            Point2D(0.6106953103224745, 0.3922045657836621),
            Point2D(0.2729193779508884, 0.6756029694480524),
            Point2D(-0.5609710546488887, 0.5657805033423454),
            Point2D(-0.26656002963029274, 0.3923815453043719),
            Point2D(0.14322739348953908, 0.15178785732278627),
            Point2D(0.39172907469276025, 0.4649650266278157),
            Point2D(0.4500353825397669, 0.02833073135506423),
            Point2D(0.1441893023975679, 0.5489509459934137),
            Point2D(-0.1152357517458235, 0.7092795091829979)]
    for i = 1:size(pos,1)
        add!(m, circle(d_in, pos[i], seed = seed_in, dir = "neg"))
    end
    add!(m, circle(d_out, center, seed = seed_out, dir = "pos"))
    offset_model!(m)
    mi = make_model_immutable(m)
    println("model has ", mi.nelem, " elements")
    return mi
end

function model_circles_in_circle_rand_quarter(d_in, d_out, elemsize)
    # circles in circle - distribution: random case quarter
    m = create_empty_model()
    seed_in = round(Integer,(pi*d_in)/elemsize)
    seed_out = round(Integer,(pi*d_out)/elemsize)
    center = Point2D(0.0,0.0)
    pos = [ Point2D(0.04653955138185156, 0.6005048683982765),
            Point2D(0.5058154728811148, 0.5154802821741737),
            Point2D(0.06354005754058685, 0.4232676434615091),
            Point2D(0.2474902955930964, 0.3449615973970197),
            Point2D(0.3706308034929659, 0.21933793063608178),
            Point2D(0.3290114648791275, 0.589939086034762),
            Point2D(0.23014554966526696, 0.0914868377927173),
            Point2D(0.47830449009160825, 0.3670302791427388),
            Point2D(0.1131074869381231, 0.27490908707385225),
            Point2D(0.5599208816828124, 0.11055418129300675),
            Point2D(0.6308260961526087, 0.27613890403113395),
            Point2D(0.42308992347522917, 0.016272533718979004),
            Point2D(0.09572990827015943, 0.7425117009982518),
            Point2D(0.7231774797046068, 0.12978353888817498),
            Point2D(0.027986933711252673, 0.15126801812809365),
            Point2D(0.18759948878965407, 0.5348799473499365),
            Point2D(0.2496608248774685, 0.7253070400849649)]
    for i = 1:size(pos,1)
        add!(m, circle(d_in, pos[i], seed = seed_in, dir = "neg"))
    end
    add!(m, circle(d_out, center, seed = seed_out, dir = "pos"))
    offset_model!(m)
    mi = make_model_immutable(m)
    println("model has ", mi.nelem, " elements")
    return mi
end

function model_line_to_line_with_obstacles(a, b, c, n, elemsize)
    # line to line with line obstacles
    # a - length of start and end line
    # b - distance between start and end line
    # c - distance between start and obstacle line
    # n - number of obstacle lines
    m = create_empty_model()
    seeda = round(Integer,a/elemsize)
    point1 = Point2D(0.0,0.0)
    point2 = Point2D(a,0.0)
    add!(m, edge(point1, point2, seed = seeda, dir = "pos"))
    point1 = Point2D(0.0,b)
    point2 = Point2D(a,b)
    add!(m, edge(point1, point2, seed = seeda, dir = "neg"))
    d = a / (2*n-1)
    seedd = round(Integer,d/elemsize)
    e = 0.0
    for i = 1:n
        point1 = Point2D(e,c)
        e += d
        point2 = Point2D(e,c)
        add!(m, edge(point1, point2, seed = seedd, dir = "neg"))
        e += d
    end
    offset_model!(m)
    mi = make_model_immutable(m)
    println("model has ", mi.nelem, " elements")
    return mi
end


text:
#script that estimates the dir bin gas network model on emid data and evaluates GAS
# forecasting performances

using Utilities,AReg,StaticNets,JLD,MLBase,StatsBase,DynNets
using PyCall; pygui(:qt); using PyPlot
#estimate and save for half the datase (after LTRO) or whole data?
halfPeriod = false

## Load dataj
fold_Path =  "/home/Domenico/Dropbox/Dynamic_Networks/data/emid_data/juliaFiles/"
loadFilePartialName = "Weekly_eMid_Data_from_"
halfPeriod? periodEndStr =  "2012_03_12_to_2015_02_27.jld": periodEndStr =  "2009_06_22_to_2015_02_27.jld"
@load(fold_Path*loadFilePartialName*periodEndStr, AeMidWeekly_T,banksIDs,inactiveBanks,
                YeMidWeekly_T,weekInd,datesONeMid ,degsIO_T,strIO_T)

#

## Load  piero's DATA
using MAT
pieroEmidEst = matread("/home/Domenico/Dropbox/Dynamic_Networks/data/emid_data/matlabFiles/data\ Piero/forecastEMID/fitness_timeseries.mat")
pieroEmidData = matread("/home/Domenico/Dropbox/Dynamic_Networks/data/emid_data/matlabFiles/data\ Piero/forecastEMID/eMID_weeklyaggregated_binaryDIRECT_postLTRO.mat")

#density plot
matY_T = YeMidWeekly_T[:,:,3:end]
 matA_T = matY_T.>0
    N = length(matA_T[1,:,1])
    T = size(matA_T)[3]
    Ttrain = 100#round(Int, T/2) #70 106 #
    threshVar = 0.00#5
    degsIO_T = [sumSq(matA_T,2);sumSq(matA_T,1)]
 allFitSS =  StaticNets.estimate( StaticNets.SnapSeqNetDirBin1(degsIO_T); identPost = false,identIter= true )

density_T = [sum(matA_T[:,:,t]) for t=1:T]./(N*(N-1))
sizeLab = 30
 close()
 subplot(2,1,1) ;plot(density_T)
 ylabel("Network Density",size = sizeLab)
 #xlabel("Time",size = sizeLab)
 grid()
 plt[:axvline](x=Ttrain,color = "r",linestyle = "--")
 S_T = [sum(matY_T[:,:,t]) for t=1:T]
 subplot(2,1,2);plot(S_T)
  ylabel("Total Volume",size = sizeLab)
  xlabel("Time",size = sizeLab)
  grid()
 plt[:axvline](x=Ttrain,color = "r",linestyle = "--")

#Estimate gas model on train sample
allFitConstTrain,~,~ =  StaticNets.estimate( StaticNets.ErgmDirBin1(meanSq(degsIO_T[:,1:Ttrain],2)) )
 modGasDirBin1_eMidTrain = DynNets.SdErgmDirBin1(degsIO_T[:,1:Ttrain],"FISHER-DIAG")
 estTargDirBin1_eMidTrain,~ = DynNets.estimateTarg(modGasDirBin1_eMidTrain;SSest = allFitSS )
 gasParEstOnTrain = estTargDirBin1_eMidTrain
 GasforeFit,~ = DynNets.score_driven_filter_or_dgp( DynNets.SdErgmDirBin1(degsIO_T),[gasParEstOnTrain[1];gasParEstOnTrain[2];gasParEstOnTrain[3]])
 gasforeFit = Float64.(GasforeFit)

#foreEval = DynNets.forecastEvalGasNetDirBin1(matA_T,gasParEstOnTrain,Ttrain;thVarConst = threshVar,mat2rem = remMat)
close()
 #select links among largest banks in training sample
 remMat =falses(N,N)#squeeze(prod(.!matA_T,3),3)#
 sizeLeg = 15
 tmpTm1 =  StaticNets.nowCastEvalFitNet(  gasforeFit[:,1:T] ,matA_T,Ttrain;mat2rem = remMat,shift = 0)
    legTex = ["GAS t-1 AUC = $(round(tmpTm1[3],3))"]
    # ROC for predictions from mean over Ttrain estimates
    tmpTm1 = StaticNets.nowCastEvalFitNet(   allFitSS[:,1:end] ,matA_T,Ttrain;mat2rem = remMat)
    legTex = [legTex; "t Par  AUC = $(round(tmpTm1[3],3))"]
    # tmpTm1 =  StaticNets.nowCastEvalFitNet(  gasforeFit[:,1:T] ,matA_T,Ttrain;mat2rem = remMat,shift = -1)
    # legTex = [legTex;"GAS t AUC = $(round(tmpTm1[3],3))"]
    shiftVal = 1
    tmpTm1 = StaticNets.nowCastEvalFitNet(   allFitSS[:,1:end] ,matA_T,Ttrain;mat2rem = remMat,shift = shiftVal)
    legTex = [legTex; "t-$(shiftVal) Par  AUC = $(round(tmpTm1[3],3))"]

    tmpTm1 = StaticNets.nowCastEvalFitNet( repmat( allFitConstTrain,1,T),matA_T,Ttrain;mat2rem = remMat)
    legTex = [legTex; "Const Par  AUC = $(round(tmpTm1[3],3))"]

    tmpTm1 = StaticNets.forecastEvalAR1Net( matA_T,Ttrain,mat2rem = remMat)
    legTex = [legTex; "AR1 SS = $(round(tmpTm1[3],3))"]
    legend(legTex,fontsize = sizeLeg)
    title("T train = $(Ttrain) of $(T)  ",size = sizeLab)

    xlabel("False Positive Rate",size = sizeLab)
    ylabel("False Negative Rate",size = sizeLab)

# Piero estimates of the fitnesses
N2,T = size(pieroEmidEst["thetaEST0"]);N =round(Int,N2/2);TtrainP = T - 50
 fitP1_T = -pieroEmidEst["thetaEST0x"]# [pieroEmidEst["thetaEST0"][1:N,:];pieroEmidEst["thetaEST0"][1:N,:]]
 fitP2_T = -pieroEmidEst["thetaEST0"]
 matAP_T = pieroEmidData["A"].>0
 shiftVal = 1
 mat2rem = falses(N,N)# squeeze(prod(.!matAP_T,3),3)
 tmpTm1 = StaticNets.nowCastEvalFitNet(fitP2_T,matAP_T,TtrainP,shift = shiftVal,mat2rem=mat2rem)
 legTex = ["theta0 t - $(shiftVal)   AUC = $(round(tmpTm1[3],3))"]
 tmpTm1 = StaticNets.nowCastEvalFitNet(fitP1_T,matAP_T,TtrainP,shift = shiftVal,mat2rem=mat2rem)
 legTex = [legTex ; "theta0X t - $(shiftVal)   AUC = $(round(tmpTm1[3],3))"]
 legend(legTex)

N = length(matAP_T[1,:,1])
    T = size(matAP_T)[3]
    Ttrain = 106#round(Int, T/2) #70 106 #
    threshVar = 0.00#5
    degsIOP_T = [sumSq(matAP_T,2);sumSq(matAP_T,1)]

 allFitSS =   StaticNets.estimate( StaticNets.SnapSeqNetDirBin1(degsIOP_T); identPost = false,identIter= true,targetErr = 1e-6 )
 allFitConstTrain,~,~ =  StaticNets.estimate( StaticNets.ErgmDirBin1(meanSq(degsIOP_T[:,1:Ttrain],2)) )
 modGasDirBin1_eMidTrain = DynNets.SdErgmDirBin1(degsIOP_T[:,1:Ttrain],"FISHER-DIAG")
 estTargDirBin1_eMidTrain,~ = DynNets.estimateTarg(modGasDirBin1_eMidTrain;SSest = allFitSS )
 gasParEstOnTrain = estTargDirBin1_eMidTrain
 GasforeFit,~ = DynNets.score_driven_filter_or_dgp( DynNets.SdErgmDirBin1(degsIOP_T),
                    [gasParEstOnTrain[1];gasParEstOnTrain[2];gasParEstOnTrain[3]])

 #Nconst,isConIn,isConOut = defineConstDegs(degsIOP_T[:,1:Ttrain];thVarConst = thVarConst )

 #foreEval = DynNets.forecastEvalGasNetDirBin1(matAP_T,gasParEstOnTrain,Ttrain;thVarConst = threshVar,mat2rem = remMat)
 close()
    #select links among largest banks in training sample
    remMat = mat2rem# falses(N,N)#;for i=1:N remMat[i,i] = true; end
    tmpGas =  StaticNets.nowCastEvalFitNet(   Float64.(GasforeFit) ,matAP_T,Ttrain;mat2rem = remMat,shift=0)
     legTex = ["GAS  AUC = $(round(tmpGas[3],3))"]
     legend(legTex)
     #tmpAR = StaticNets.forecastEvalAR1Net(matAP_T,Ttrain;thVarConst = threshVar,mat2rem = remMat)
     #legTex = [legTex ; "AR1  AUC = $(round(tmpAR[3],3))"]
     legend(legTex)
     title("T train = $(Ttrain) of $(T)  ")
       tmpTm1 = StaticNets.nowCastEvalFitNet(   allFitSS[:,1:end] ,matAP_T,Ttrain;mat2rem = remMat)
     legTex = [legTex; "t Par  AUC = $(round(tmpTm1[3],3))"]
     # tmpTm1 = StaticNets.nowCastEvalFitNet(  [ allFitSS[:,1]  allFitSS[:,1:end-1] ] ,matAP_T,Ttrain;mat2rem = remMat)
     # legTex = [legTex; "t-1 Par  AUC = $(round(tmpTm1[3],3))"]
     shiftVal = 1
     tmpTm1 = StaticNets.nowCastEvalFitNet(    allFitSS[:,1:end],matAP_T,Ttrain;mat2rem = remMat,shift=shiftVal)
     legTex = [legTex; "t-$(shiftVal) Par  AUC = $(round(tmpTm1[3],3))"]

     tmpTm1 = StaticNets.nowCastEvalFitNet( repmat( allFitConstTrain,1,T),matAP_T,Ttrain;mat2rem = remMat)
     legTex = [legTex; "Const Par  AUC = $(round(tmpTm1[3],3))"]

     tmpTm1 = StaticNets.forecastEvalAR1Net( matAP_T,Ttrain;mat2rem = remMat)
     legTex = [legTex; "AR1 SS = $(round(tmpTm1[3],3))"]
     legend(legTex)
     title("T train = $(Ttrain) of $(T)  ")
     grid()
     xlabel("False Positive Rate")
     ylabel("False Negative Rate")

using RCall

function Rroc(score,response)
    R"library('pROC')"
        @rput(score)
        @rput(response)
        R" plot(roc(response, score , direction='<'), print.auc=TRUE)
        par(new=TRUE)
        "
end

Rroc(tmpGas[5], tmpGas[4])
Rroc(tmpTm1[5], tmpTm1[4])

#

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
"""
Julia program to check if a number is a pronic number
Pronic Numbers are those numbers that can be arranged to form a rectangle.
It is a number which is the product of two consecutive integer
"""

function check_pronic(x)
    i = 0
    while(i <= isqrt(x))
        if(x == i * (i + 1))
            return true
        end
        i = i + 1
    end
    return false
end

print("Enter the number: ")
n = readline()
n = parse(Int, n)
res = check_pronic(abs(n))
if(res)
    print("The given number $n is a pronic number.")
else
    print("The given number $n is not a pronic number.")
end

"""
Time Complexity: O(log(n)), where 'n' is the given number
Space Complexity: O(1)

SAMPLE INPUT AND OUTPUT

SAMPLE 1

Enter the number: 154
The given number 154 is not a pronic number.

SAMPLE 2
Enter the number: 72
The given number 72 is a pronic number.
"""

max_issues_repo_name: adarshnjena/NeoAlgo

max_forks_repo_name: adarshnjena/NeoAlgo


text:
using ConsistencyResampling, Bootstrap
using Random

using ConsistencyResampling: bootstrap_direct, bootstrap_alias

Random.seed!(1234)

@testset "Simple example" begin
    predictions = reshape([1, 0], 2, 1)
    labels = [2]

    for f in (bootstrap, bootstrap_direct, bootstrap_alias)
        b = bootstrap(last, (predictions, labels), ConsistentSampling(20))

        @test data(b) == (predictions, labels)
        @test original(b) == (2,)
        @test straps(b) == (ones(Int, 20),)
    end
end


text:
###
# Here we implement inversion of banded matrices, beginning with Bidiagonal
###

function SemiseparableMatrix(Bi::ApplyMatrix{T,typeof(inv),<:Tuple{Bidiagonal}}) where T
    B = inv(Bi)
    n = size(B,1)
    ℓ = [1; cumprod(-B.ev./B.dv[1:end-1])]
    if B.uplo == 'L'
        SemiseparableMatrix(ApplyArray(*, inv.(B.dv) .* ℓ, inv.(ℓ)'), Zeros{T}(n,n), 0, 1)
    else
        SemiseparableMatrix(Zeros{T}(n,n), ApplyArray(*, inv.(ℓ), (inv.(B.dv) .* ℓ)'), 1, 0)
    end
end


text:
### -------------------------------------------------------------------------------
# supportPoint(solid, shape,r_abs,R_abs,e) - Return support point of shape shape at position r_abs/R_abs
#                                   in direction of unit vector e
#                                   function is called from mpr algorithm
# boundingBox!(solid, shape,r_abs,R_abs,AABB) - Update vertex positions of smallest AABB box that contains shape

"""
    r = supportPoint(solid, shape, r_abs, R_abs, e)

Return the absolute position vector `r` from world frame to *support point* of
solid `shape`, resolved in world frame in [m]. The support point is the point of solid `shape`
that is the most extreme in direction of unit vector `e`.

# Arguments
- `shape::Modia3D.AbstractGeometry`: Solid shape object.
- `r_abs::AbstractVector`: Absolute position vector of `shape` reference frame.
- `R_abs::AbstractMatrix`: Rotation matrix to rotate world frame in `shape` reference frame.
- `e::AbstractVector`: Unit vector pointing into the desired direction.
"""
@inline supportPoint_Box(shape::Box, r_abs::SVector{3,Float64}, R_abs::SMatrix{3,3,Float64,9}, e::SVector{3,Float64}, collisionSmoothingRadius::Float64) =
            r_abs + R_abs'*supportPoint_abs_Box(shape, R_abs*e, collisionSmoothingRadius) + e*collisionSmoothingRadius
#
@inline supportPoint_Cylinder( shape::Cylinder, r_abs::SVector{3,Float64}, R_abs::SMatrix{3,3,Float64,9}, e::SVector{3,Float64}, collisionSmoothingRadius::Float64) =
            r_abs + R_abs'*supportPoint_abs_Cylinder(shape, R_abs*e) + e*collisionSmoothingRadius

@inline supportPoint_Cone(shape::Cone, r_abs::SVector{3,Float64}, R_abs::SMatrix{3,3,Float64,9}, e::SVector{3,Float64}, collisionSmoothingRadius::Float64) =
            r_abs + R_abs'*supportPoint_abs_Cone(shape, R_abs*e) + e*collisionSmoothingRadius

@inline supportPoint_Capsule(shape::Capsule, r_abs::SVector{3,Float64}, R_abs::SMatrix{3,3,Float64,9}, e::SVector{3,Float64}, collisionSmoothingRadius::Float64) =
            r_abs + R_abs'*supportPoint_abs_Capsule(shape, R_abs*e) + e*collisionSmoothingRadius

@inline supportPoint_Beam(shape::Beam, r_abs::SVector{3,Float64}, R_abs::SMatrix{3,3,Float64,9}, e::SVector{3,Float64}, collisionSmoothingRadius::Float64) =
            r_abs + R_abs'*supportPoint_abs_Beam(shape, R_abs*e) + e*collisionSmoothingRadius

@inline supportPoint_Sphere(shape::Sphere, r_abs::SVector{3,Float64}, R_abs::SMatrix{3,3,Float64,9}, e::SVector{3,Float64}) =
            r_abs + (shape.diameter/2)*e

@inline supportPoint_Ellipsoid(shape::Ellipsoid, r_abs::SVector{3,Float64}, R_abs::SMatrix{3,3,Float64,9}, e::SVector{3,Float64}) =
            r_abs + R_abs'*supportPoint_abs_Ellipsoid(shape, R_abs*e)

@inline supportPoint_FileMesh(shape::FileMesh, r_abs::SVector{3,Float64}, R_abs::SMatrix{3,3,Float64,9}, e::SVector{3,Float64}) =
            r_abs + R_abs'*supportPoint_abs_FileMesh(shape, R_abs*e)

# supportPoint_abs(shape, e_abs) - Return support point of shape shape at the shape reference
#         frame of shape in direction of unit vector e_abs (resolved in shape reference frame = object3d frame)
# see e.g. [Gino v.d. Bergen, Collision Detection in Interactive 3D Environments, 2003]
# or [G. Snethen, Xenocollide: Complex collision made simple. In Game Programming Gems 7, 2008]
# or Andrea Neumayr took ideas for computing based on those books

# [Gino v.d. Bergen, p. 135]
@inline supportPoint_abs_Ellipsoid(shape::Ellipsoid, e_abs::SVector{3,Float64}) =
    @inbounds SVector( (shape.lengthX/2)^2*e_abs[1], (shape.lengthY/2)^2*e_abs[2], (shape.lengthZ/2)^2*e_abs[3] )/norm(SVector((shape.lengthX/2)*e_abs[1], (shape.lengthY/2)*e_abs[2], (shape.lengthZ/2)*e_abs[3]) )

# [Gino v.d. Bergen, p. 135]
@inline supportPoint_abs_Box(shape::Box, e_abs::SVector{3,Float64},  collisionSmoothingRadius::Float64) =
    @inbounds SVector( Basics.sign_eps(e_abs[1])*(shape.lengthX/2-collisionSmoothingRadius), Basics.sign_eps(e_abs[2])*(shape.lengthY/2-collisionSmoothingRadius), Basics.sign_eps(e_abs[3])*(shape.lengthZ/2-collisionSmoothingRadius) )

# [Gino v.d. Bergen, p. 136, XenoCollide, p. 168, 169]
@inline function supportPoint_abs_Cylinder(shape::Cylinder, e_abs::SVector{3,Float64})
    @inbounds begin
        if shape.axis == 1
            enorm = norm(SVector(e_abs[2], e_abs[3]))
            if enorm <= Modia3D.neps
                return Basics.sign_eps(e_abs[1])*SVector(shape.length/2, 0.0, 0.0)
            else
                return Basics.sign_eps(e_abs[1])*SVector(shape.length/2, 0.0, 0.0) + SVector(0.0, (shape.diameter/2)*e_abs[2], (shape.diameter/2)*e_abs[3])/enorm
            end
        elseif shape.axis == 2
            enorm = norm(SVector(e_abs[3], e_abs[1]))
            if enorm <= Modia3D.neps
                return Basics.sign_eps(e_abs[2])*SVector(0.0, shape.length/2, 0.0)
            else
                return Basics.sign_eps(e_abs[2])*SVector(0.0, shape.length/2, 0.0) + SVector((shape.diameter/2)*e_abs[1], 0.0, (shape.diameter/2)*e_abs[3])/enorm
            end
        else
            enorm = norm(SVector(e_abs[1], e_abs[2]))
            if enorm <= Modia3D.neps
                return Basics.sign_eps(e_abs[3])*SVector(0.0, 0.0, shape.length/2)
            else
                return Basics.sign_eps(e_abs[3])*SVector(0.0, 0.0, shape.length/2) + SVector((shape.diameter/2)*e_abs[1], (shape.diameter/2)*e_abs[2], 0.0)/enorm
            end
        end
    end
end

# G. Hippmann: Cylinder + Sphere as bottom and top
@inline function supportPoint_abs_Capsule(shape::Capsule, e_abs::SVector{3,Float64})
    @inbounds begin
        if shape.axis == 1
            return Basics.sign_eps(e_abs[1])*SVector(0.5*shape.length, 0.0, 0.0) + 0.5*shape.diameter*e_abs
        elseif shape.axis == 2
            return Basics.sign_eps(e_abs[2])*SVector(0.0, 0.5*shape.length, 0.0) + 0.5*shape.diameter*e_abs
        else
            return Basics.sign_eps(e_abs[3])*SVector(0.0, 0.0, 0.5*shape.length) + 0.5*shape.diameter*e_abs
        end
    end
end

# for cone: [Gino v.d. Bergen, p. 136]]
# for frustum of a cone: A. Neumayr, G. Hippmann
@inline function supportPoint_abs_Cone(shape::Cone, e_abs::SVector{3,Float64})
    @inbounds begin
        baseRadius = shape.diameter/2
        rightCone = shape.topDiameter == 0.0
        if rightCone
            sin_phi = baseRadius/sqrt(baseRadius^2 + shape.length^2)  # sin of base angle
        else
            topRadius = shape.topDiameter/2
            diffRadius = baseRadius - topRadius
            sin_phi = diffRadius/sqrt(diffRadius^2 + shape.length^2)  # sin of base angle
        end
        if shape.axis == 1
            value = e_abs[1] / norm(SVector(e_abs[1], e_abs[2], e_abs[3]))
            if value >= sin_phi
                if rightCone
                    return SVector(shape.length, 0.0, 0.0)  # apex is support point
                else  # frustum of a cone
                    enorm = norm(SVector(e_abs[2], e_abs[3]))
                    if enorm > Modia3D.neps
                        return SVector(shape.length, 0.0, 0.0) + SVector(0.0, topRadius*e_abs[2], topRadius*e_abs[3]) / enorm  # point on top circle is support point
                    else
                        return SVector(shape.length, 0.0, 0.0)  # top circle center is support point
                    end
                end
            else
                enorm = norm(SVector(e_abs[2], e_abs[3]))
                if enorm > Modia3D.neps
                    return SVector(0.0, baseRadius*e_abs[2], baseRadius*e_abs[3]) / enorm  # point on base circle is support point
                else
                    return SVector(0.0, 0.0, 0.0)  # base circle center is support point
                end
            end
        elseif shape.axis == 2
            value = e_abs[2] / norm(SVector(e_abs[1], e_abs[2], e_abs[3]))
            if value >= sin_phi
                if rightCone
                    return SVector(0.0, shape.length, 0.0)  # apex is support point
                else  # frustum of a cone
                    enorm = norm(SVector(e_abs[3], e_abs[1]))
                    if enorm > Modia3D.neps
                        return SVector(0.0, shape.length, 0.0) + SVector(topRadius*e_abs[1], 0.0, topRadius*e_abs[3]) / enorm  # point on top circle is support point
                    else
                        return SVector(0.0, shape.length, 0.0)  # top circle center is support point
                    end
                end
            else
                enorm = norm(SVector(e_abs[3], e_abs[1]))
                if enorm > Modia3D.neps
                    return SVector(baseRadius*e_abs[1], 0.0, baseRadius*e_abs[3]) / enorm  # point on base circle is support point
                else
                    return SVector(0.0, 0.0, 0.0)  # base circle center is support point
                end
            end
        else
            value = e_abs[3] / norm(SVector(e_abs[1], e_abs[2], e_abs[3]))
            if value >= sin_phi
                if rightCone
                    return SVector(0.0, 0.0, shape.length)  # apex is support point
                else  # frustum of a cone
                    enorm = norm(SVector(e_abs[1], e_abs[2]))
                    if enorm > Modia3D.neps
                        return SVector(0.0, 0.0, shape.length) + SVector(topRadius*e_abs[1], topRadius*e_abs[2], 0.0) / enorm  # point on top circle is support point
                    else
                        return SVector(0.0, 0.0, shape.length)  # top circle center is support point
                    end
                end
            else
                enorm = norm(SVector(e_abs[1], e_abs[2]))
                if enorm > Modia3D.neps
                    return SVector(baseRadius*e_abs[1], baseRadius*e_abs[2], 0.0) / enorm  # point on base circle is support point
                else
                    return SVector(0.0, 0.0, 0.0)  # base circle center is support point
                end
            end
        end
    end
end

# G. Hippmann: Outer half circles of beam
@inline function supportPoint_abs_Beam(shape::Beam, e_abs::SVector{3,Float64})
    @inbounds begin
        if shape.axis == 1
            enorm = norm(SVector(e_abs[1], e_abs[2]))
            if enorm <= Modia3D.neps
                return SVector(Basics.sign_eps(e_abs[1])*shape.length/2, 0.0, Basics.sign_eps(e_abs[3])*shape.thickness/2)
            else
                return SVector(Basics.sign_eps(e_abs[1])*shape.length/2, 0.0, Basics.sign_eps(e_abs[3])*shape.thickness/2) + SVector(shape.width/2*e_abs[1], shape.width/2*e_abs[2], 0.0)/enorm
            end
        elseif shape.axis == 2
            enorm = norm(SVector(e_abs[2], e_abs[3]))
            if enorm <= Modia3D.neps
                return SVector(Basics.sign_eps(e_abs[1])*shape.thickness/2, Basics.sign_eps(e_abs[2])*shape.length/2, 0.0)
            else
                return SVector(Basics.sign_eps(e_abs[1])*shape.thickness/2, Basics.sign_eps(e_abs[2])*shape.length/2, 0.0) + SVector(0.0, shape.width/2*e_abs[2], shape.width/2*e_abs[3])/enorm
            end
        else
            enorm = norm(SVector(e_abs[3], e_abs[1]))
            if enorm <= Modia3D.neps
                return SVector(0.0, Basics.sign_eps(e_abs[2])*shape.thickness/2, Basics.sign_eps(e_abs[3])*shape.length/2)
            else
                return SVector(0.0, Basics.sign_eps(e_abs[2])*shape.thickness/2, Basics.sign_eps(e_abs[3])*shape.length/2) + SVector(shape.width/2*e_abs[1], 0.0, shape.width/2*e_abs[3])/enorm
            end
        end
    end
end

# [Gino v.d. Bergen, p. 131]
@inline function supportPoint_abs_FileMesh(shape::FileMesh, e_abs::SVector{3,Float64})
    @inbounds begin
        max_value = Float64
        position = Float64
        for i in eachindex(shape.objPoints)
            if i == 1
            max_value = dot(shape.objPoints[i],e_abs)
            position = i
            else
            act_value = dot(shape.objPoints[i],e_abs)
            if act_value > max_value
                max_value = act_value
                position = i
        end; end; end
        return shape.objPoints[position]
    end
end

#------------------------------------------------------------------------------------
# boundingBox!(..) calculates the AABB (= Axis Aligned Bounding Box) for the geometries
# it also uses support points for finding the points which are furthest away in each x,y,z - direction
# it calles supportPoint_i
"""
    boundingBox!(shape, AABB, r_abs, R_abs; tight=true, scaleFactor=0.01)

Returns the *Axis Aligned Bounding Box* of solid `shape` in argument `AABB`.

# Arguments
- `shape::Modia3D.AbstractGeometry`: Solid shape object.
- `AABB::Modia3D.BoundingBox`: On return, update AABB with the actual *Axis Aligned Bounding Box* of solid `shape`.
- `r_abs::AbstractVector`: Absolute position vector of `shape` reference frame.
- `R_abs::AbstractMatrix`: Rotation matrix to rotate world frame in `shape` reference frame.
- `tight::Bool`: If true, return the tightest AABB. If false return an AABB that is
                scaleFactor bigger than the best fitting AABB
                (for example, scaleFactor=0.1 means that the returned AABB is 10 percent bigger
                than the best fitting AABB).
- `scaleFactor::Float64`: If `tight=false`, the returned AABB is `scaleFactor` bigger than the
                          best fitting AABB.
"""
#=
function boundingBox!(solid, shape::Modia3D.AbstractGeometry, AABB::Basics.BoundingBox, r_abs::SVector{3,Float64}, R_abs::SMatrix{3,3,Float64,9}; tight::Bool=true, scaleFactor::Float64=0.01)
    @inbounds begin
        xmin = supportPoint_i(solid, shape, r_abs[1], SVector(R_abs[:,1]), -1)
        xmax = supportPoint_i(solid, shape, r_abs[1], SVector(R_abs[:,1]), +1)
        ymin = supportPoint_i(solid, shape, r_abs[2], SVector(R_abs[:,2]), -1)
        ymax = supportPoint_i(solid, shape, r_abs[2], SVector(R_abs[:,2]), +1)
        zmin = supportPoint_i(solid, shape, r_abs[3], SVector(R_abs[:,3]), -1)
        zmax = supportPoint_i(solid, shape, r_abs[3], SVector(R_abs[:,3]), +1)

        if tight    # best fitting AABB which is possible
            AABB.x_min = xmin
            AABB.x_max = xmax
            AABB.y_min = ymin
            AABB.y_max = ymax
            AABB.z_min = zmin
            AABB.z_max = zmax
        else
            # scales each axis with scaleFactor 0.1 means 10 percent bigger than best fitting AABB and so on
            @assert(scaleFactor > 0.0)
            scale = max(abs(xmax - xmin)*scaleFactor/2,abs(ymax - ymin)*scaleFactor/2,abs(zmax - zmin)*scaleFactor/2)
            AABB.x_min = xmin - scale
            AABB.x_max = xmax + scale
            AABB.y_min = ymin - scale
            AABB.y_max = ymax + scale
            AABB.z_min = zmin - scale
            AABB.z_max = zmax + scale
        end
        return AABB
    end
 #=
  # other idea for AABB
  # scales each axis with scaleFactor 0.1 means 10 percent bigger than best fitting AABB and so on
  scaleX = abs(xmax - xmin)*scaleFactor/2
  scaleY = abs(ymax - ymin)*scaleFactor/2
  scaleZ = abs(zmax - zmin)*scaleFactor/2
  AABB.x_min = xmin - scaleX
  AABB.x_max = xmax + scaleX
  # so on

  # FileMesh
   longestEdge = max(shape.scaleFactor[1], shape.scaleFactor[2], shape.scaleFactor[3]) * shape.longestEdge * sqrt(2)
   x_min = -shape.longestEdge + r_rel[1]
   x_max =  shape.longestEdge + r_rel[1]
   y_min = -shape.longestEdge + r_rel[2]
   y_max =  shape.longestEdge + r_rel[2]
   z_min = -shape.longestEdge + r_rel[3]
   z_max = shape.longestEdge + r_rel[3]
 =#
end

function boundingBox!(solid, shape::Sphere, AABB::Basics.BoundingBox, r_abs::SVector{3,Float64}, R_abs::SMatrix{3,3,Float64,9}; tight::Bool=true, scaleFactor::Float64=0.01)
    @inbounds begin
        r = shape.diameter/2
        if tight    # best fitting AABB which is possible
            AABB.x_min = r_abs[1] - r
            AABB.x_max = r_abs[1] + r
            AABB.y_min = r_abs[2] - r
            AABB.y_max = r_abs[2] + r
            AABB.z_min = r_abs[3] - r
            AABB.z_max = r_abs[3] + r
        else
            # scales each axis with scaleFactor 0.1 means 10 percent bigger than best fitting AABB and so on
            @assert(scaleFactor > 0.0)
            scale = r*scaleFactor/2
            AABB.x_min = r_abs[1] - r - scale
            AABB.x_max = r_abs[1] + r + scale
            AABB.y_min = r_abs[2] - r - scale
            AABB.y_max = r_abs[2] + r + scale
            AABB.z_min = r_abs[3] - r - scale
            AABB.z_max = r_abs[3] + r + scale
        end
        return AABB
    end
end
=#
#------------------------------------------------------------------------------------
# supportPoint_i are utility functions for boundingBox!(..)
# it is based on supportPoint_abs (its nearly the same as supportPoint(..))
@inline supportPoint_i_Box(shape::Box, r_absi::Float64, R_absi::SVector{3,Float64}, isign::Int, collisionSmoothingRadius) =
               r_absi + R_absi'*supportPoint_abs_Box(shape, isign*R_absi, collisionSmoothingRadius) + isign*collisionSmoothingRadius

@inline supportPoint_i_Cylinder(shape::Cylinder, r_absi::Float64, R_absi::SVector{3,Float64}, isign::Int, collisionSmoothingRadius) =
               r_absi + R_absi'*supportPoint_abs_Cylinder(shape, isign*R_absi) + isign*collisionSmoothingRadius

@inline supportPoint_i_Cone(shape::Cone, r_absi::Float64, R_absi::SVector{3,Float64}, isign::Int, collisionSmoothingRadius) =
               r_absi + R_absi'*supportPoint_abs_Cone(shape, isign*R_absi) + isign*collisionSmoothingRadius

@inline supportPoint_i_Capsule(shape::Capsule, r_absi::Float64, R_absi::SVector{3,Float64}, isign::Int, collisionSmoothingRadius) =
               r_absi + R_absi'*supportPoint_abs_Capsule(shape, isign*R_absi) + isign*collisionSmoothingRadius
#
@inline supportPoint_i_Beam(shape::Beam, r_absi::Float64, R_absi::SVector{3,Float64}, isign::Int, collisionSmoothingRadius) =
               r_absi + R_absi'*supportPoint_abs_Beam(shape, isign*R_absi) + isign*collisionSmoothingRadius

@inline supportPoint_i_Ellipsoid(shape::Ellipsoid, r_absi::Float64, R_absi::SVector{3,Float64}, isign::Int) =
               r_absi + dot(R_absi, supportPoint_abs_Ellipsoid(shape, isign*R_absi))

@inline supportPoint_i_FileMesh(shape::FileMesh, r_absi::Float64, R_absi::SVector{3,Float64}, isign::Int) =
               r_absi + dot(R_absi, supportPoint_abs_FileMesh(shape, isign*R_absi))


text:
using JEMSS

println("Loading sim.")
sim = initSim(joinpath(@__DIR__, "input/sim_config.xml"))

println("Simulating.")
simulate!(sim)

println("Writing stats files to \"$(abspath(sim.outputPath))\"")
writeStatsFiles(sim)

# read some stats from file
println("Some statistics from stats_dict.csv output file:")
statsDict = readStatsDictFile(sim.outputFiles["statsDict"].path)
@show statsDict["ambs_avgDailyBusyDurationHours"]
@show statsDict["calls_all_fracResponsesInTime"]

# tables = readTablesFromFile(sim.outputFiles["ambulancesStats"].path)
# sim.outputFiles["callsStats"].path

# # Animate, if you like:
# animate!(sim)

nothing


text:
using SpecialFunctions

mutable struct MPCValues{T}
    # path information
    path            ::SplinePath

    # problem parameters
    n_modes         ::Int64
    N               ::Int64 # MPC horizon
    k_c             ::Int64 # consensus horizon
    dt              ::T # time step
    n               ::Int64 # state dimension (augmented to include path parameter s)
    m               ::Int64 # control dimension
    n_obs           ::Int64

    # initialization
    q0              ::Vector{T}
    u0              ::Vector{T}

    # objective function weights
    qc              ::T # contouring cost
    qcterm          ::T # terminal contouring cost
    ql              ::T # lag cost
    γ               ::T # progress reward
    rΔω             ::T # rate of change of yaw rate
    rΔa             ::T # rate of change of acceleration
    rΔvs            ::T # rate of change of path velocity
    R               ::Array{T, 2} # control penalty matrix
    L               ::Vector{T} # progress reward vector

    # tracking error Taylor expansion linear and quadratic matrices
    cs              ::Vector{Vector{T}}
    Γs              ::Vector{Array{T, 2}}
    c               ::Vector{T}
    Γ               ::Array{T, 2}

    # discrete, linearized, time-varying dynamics matrices
    # x_{k+1} = g_{k} + A_{q,u,k}*Δq_{q,u,k} + B_{q,u,k}*Δu_{q,u,k}
    gs              ::Vector{Vector{T}} # linearization point
    As              ::Vector{Array{T, 2}} # dynamics matrices
    Bs              ::Vector{Array{T, 2}} # input matrices

    # information for obstacle constraints
    ps              ::Vector{Vector{T}} # robot position vectors
    p_obs           ::Dict              # obstacle positions
    obs_on          ::Vector{Bool}      # obstacle constraint indicator array
    obs_horizon     ::Int64             # obstacle existence horizon
end

function MPCValues(path::SplinePath; n_modes = 1,
                                     N = 40,
                                     k_c = 4,
                                     dt = 0.25,
                                     n = 5,
                                     m = 3,
                                     n_obs = 1,
                                     q0 = [0.; 0.; 0.0; 1.0; 0.0],
                                     u0 = [0.; 0.; 0.],
                                     qc = 5.,
                                     qcterm = 5.,
                                     ql = 5.,
                                     γ = 0.01,
                                     rΔω = 0.01,
                                     rΔa = 0.01,
                                     rΔvs = 0.01,
                                     obs_horizon = 13)

    S_q = k_c + n_modes * (N - k_c)
    S_u = S_q - 1
    S_d = S_q - n_modes

    MPCValues(path, n_modes, N, k_c, dt, n, m, n_obs, q0, u0, qc, qcterm, ql, γ,
        rΔω, rΔa, rΔvs,
        diagm(repeat([rΔω, rΔa, rΔvs], S_u - n_modes)),     # R
        γ * dt * ones(S_u),                                 # L
        [zeros(3) for i in 1:S_q],                          # cs
        [zeros(3, 3) for i in 1:S_q],                       # Γs
        zeros(3 * S_q),                                     # c
        zeros(3 * S_q, 3 * S_q),                            # Γ
        [zeros(n) for i in 1:S_d],                          # gs
        [zeros(n, n) for i in 1:S_d],                       # As
        [zeros(n, m) for i in 1:S_d],                       # Bs
        [zeros(2) for i in 1:S_q],                          # ps
        Dict(),                                             # p_obs
        [false for i in 1:n_obs],                           # obs_on
        obs_horizon)                                        # obs_horizon
end

struct MPCVariables
    q       ::Array{Variable, 2}
    u       ::Array{Variable, 2}
end

mutable struct MPCParams
    qc      ::Parameter
    qcterm  ::Parameter
    ql      ::Parameter
    γ       ::Parameter
    R       ::Parameter
    L       ::Parameter
    dt      ::Parameter

    q0      ::Parameter
    u0      ::Parameter

    ω_max   ::Parameter
    ω_min   ::Parameter
    a_max   ::Parameter
    a_min   ::Parameter
    vs_max  ::Parameter
    vs_min  ::Parameter

    As
    Bs
    gs

    c       ::Parameter
    Γ       ::Parameter

    ps
    p_obs
    obs_on
end

"""
Produces an initial trajectory for the first MPC iteration.

NOTE: Check heading and wrap to π potential issues.
"""
function initial_guess(vals::MPCValues; v0=nothing)
    path = vals.path
    n_modes = vals.n_modes
    N = vals.N
    k_c = vals.k_c
    dt = vals.dt
    if v0 == nothing
        v0 = vals.q0[4]
    else
        v0 = v0
    end
    s0 = vals.q0[5]

    x_coefs = path.x_coefs
    y_coefs = path.y_coefs
    breaks = path.breaks

    s_guess = [range(s0, step=v0*dt, length=N);] #TODO: careful, if v0 is 0, this just just produces zeros
    spline_idces = find_spline_interval.(s_guess, (path,))

    X_guess = spline_x.(s_guess, (path,), spline_idces)
    Y_guess = spline_y.(s_guess, (path,), spline_idces)
    ψ_guess = heading.(s_guess, (path,), spline_idces)
    v_guess = v0 * ones(N)

    ω_guess = diff(ψ_guess) / dt
    a_guess = zeros(N-1)
    vs_guess = v0 * ones(N-1)

    qs = collect(hcat(X_guess, Y_guess, ψ_guess, v_guess, s_guess)')
    us = collect(hcat(ω_guess, a_guess, vs_guess)')

    qs = [qs repeat(qs[:, k_c+1:end], 1, n_modes - 1)]
    us = [us repeat(us[:, k_c:end], 1, n_modes - 1)]

    return qs, us
end

function construct_problem(dynamics::DynamicsModel, vals::MPCValues, scene, qs, us; verbose=false) where {T}
    # create model
    optimizer = OSQP.Optimizer(verbose=verbose)
    model = Model(optimizer)

    # convenience definitions
    n_modes = vals.n_modes
    N = vals.N
    k_c = vals.k_c
    n = vals.n
    m = vals.m

    # MPC parameters
    qc = Parameter(() -> vals.qc, model)
    qcterm = Parameter(() -> vals.qcterm, model)
    ql = Parameter(() -> vals.ql, model)
    γ = Parameter(() -> vals.γ, model)
    R = Parameter(() -> vals.R, model)
    L = Parameter(() -> vals.L, model)
    dt = Parameter(() -> vals.dt, model)
    q0 = Parameter(() -> vals.q0, model)
    u0 = Parameter(() -> vals.u0, model)

    # dynamics parameters
    ω_max = Parameter(() -> dynamics.u_limits.ω_max, model)
    ω_min = Parameter(() -> dynamics.u_limits.ω_min, model)
    a_max = Parameter(() -> dynamics.u_limits.a_max, model)
    a_min = Parameter(() -> dynamics.u_limits.a_min, model)
    vs_max = Parameter(() -> dynamics.u_limits.vs_max, model)
    vs_min = Parameter(() -> dynamics.u_limits.vs_min, model)

    # decision variables
    S_q = k_c + n_modes * (N - k_c)
    S_u = S_q - 1
    q = [Variable(model) for i in 1:n, t in 1:S_q]        # (x, y, θ, V) + s
    u = [Variable(model) for i in 1:m, t in 1:S_u]        # (ω, a) + vs

    # initial/final state constraints
    state_constraints!(model, q, q0, vals)

    # dynamics constraints
    control_limits = [ω_max, ω_min, a_max, a_min, vs_max, vs_min]
    control_constraints!(model, q, u, control_limits, vals)
    As, Bs, gs = dynamics_constraints!(model, q, u, qs, us, vals)

    # collision avoidance constraints
    ps, p_obs, obs_on = obstacle_constraints!(model, q, qs, vals, scene)

    # contour/lag cost matrices
    tracking_cost_matrices!(qs, vals)
    c = Parameter(() -> vals.c, model)
    Γ = Parameter(() -> vals.Γ, model)

    # set up for control effort term
    Δu = diff(u, dims = 2)[:]
    # consensus portion (for transitions 1 to 2, 2 to 3, ..., k_c - 1 to k_c)
    Δu = diff(u[:, 1:k_c], dims = 2)[:]
    # transition and continuation to parallel horizons, k_c to subsequent time
    # step for each mode
    for j = 1:n_modes
        start_idx = k_c + (j - 1) * (N - k_c) + 1 # parallel horizon start
        end_idx = start_idx + N - k_c - 2 # parallel horizon end

        append!(Δu, u[:, start_idx] - u[:, k_c])
        append!(Δu, diff(u[:, start_idx:end_idx], dims = 2)[:])
    end

    # assemble objective function
    obj = @expression transpose(c) * q[[1,2,5],:][:]  +
                      transpose(q[[1,2,5],:][:]) * Γ * q[[1,2,5],:][:] +
                      transpose(Δu) * R * Δu -
                      transpose(L) * u[m, :] +
                      transpose((u[:, 1] - u0)) * diagm([vals.rΔω, vals.rΔa, vals.rΔvs]) * (u[:, 1] - u0)

    @objective(model, Minimize, obj)
    params = MPCParams(qc, qcterm, ql, γ, R, L, dt, q0, u0, ω_max, ω_min, a_max, a_min,
        vs_max, vs_min, As, Bs, gs, c, Γ, ps, p_obs, obs_on)
    # params = 0
    variables = MPCVariables(q, u)

    # model
    model, params, variables, vals
end

function relinearize!(vals::MPCValues, qs, us)
    path = vals.path
    n_modes = vals.n_modes
    N = vals.N
    k_c = vals.k_c
    n = vals.n
    m = vals.m
    dt = vals.dt
    qc = vals.qc
    qcterm = vals.qcterm
    ql = vals.ql

    S_q = k_c + n_modes * (N - k_c)
    S_d = S_q - n_modes

    # dynamic constraints
    for k in 1:S_d
        Ak, Bk = linearize_dynamics(discrete_dynamics, qs[1:n-1, k], us[1:m-1, k], vals.dt)
        vals.As[k] = Ak
        vals.Bs[k] = Bk
        vals.gs[k] = discrete_dynamics(qs[1:n-1, k], us[1:m-1, k], vals.dt) - Ak * qs[1:n-1, k] - Bk * us[1:m-1, k]
    end

    # reproject to find arc length
    for k in 1:S_q
        qs[5, k] = find_best_s(qs[:, k], path)
    end

    # contour/lag cost matrices
    for k in 1:S_q
        P0 = [qs[1, k], qs[2, k], qs[5, k]] #X, Y, s of the initial guess
        spline_idx = find_spline_interval(P0[3], vals.path)
        if k == N
            ck = tracking_linear_term(P0, vals.qcterm, vals.ql, vals.path, spline_idx)
            Γk = tracking_quadratic_term(P0, vals.qcterm, vals.ql, vals.path, spline_idx) # TODO: fix function signature
        else
            ck = tracking_linear_term(P0, vals.qc, vals.ql, vals.path, spline_idx)
            Γk = tracking_quadratic_term(P0, vals.qc, vals.ql, vals.path, spline_idx) # TODO: fix function signature
        end
        vals.cs[k] = ck
        vals.Γs[k] = Γk
    end

    vals.c = vcat(vals.cs...)
    vals.Γ = Array(blockdiag(sparse.(vals.Γs)...))

    # update positions for obstacle constraints
    for k in 1:S_q
        vals.ps[k] = qs[1:2, k]
    end
end

function update_problem!(vals::MPCValues, qs, us)
    path = vals.path
    n_modes = vals.n_modes
    N = vals.N
    k_c = vals.k_c
    n = vals.n
    m = vals.m
    dt = vals.dt
    qc = vals.qc
    qcterm = vals.qcterm
    ql = vals.ql

    S_q = k_c + n_modes * (N - k_c)

    #update initial state and prev control
    vals.q0 = qs[:, 2]
    vals.q0[5] = find_best_s(vals.q0, path)
    vals.u0 = us[:, 1]

    # shift states and controls over consensus portion
    qs[:, 1:k_c] = qs[:, 2:k_c+1]
    us[:, 1:k_c-1] = us[:, 2:k_c]

    # shift dynamics and tracking error matrices over consensus portion
    vals.As[1:k_c] = vals.As[2:k_c+1]
    vals.Bs[1:k_c] = vals.Bs[2:k_c+1]
    vals.gs[1:k_c] = vals.gs[2:k_c+1]

    vals.cs[1:k_c] = vals.cs[2:k_c+1]
    vals.Γs[1:k_c] = vals.Γs[2:k_c+1]

    for j = 1:n_modes
        horizon_start = k_c + (j - 1) * (N - k_c) + 1
        horizon_end = horizon_start + N - k_c - 1

        # shift states and controls for each horizon
        qs[:, horizon_start:horizon_end-1] = qs[:, horizon_start+1:horizon_end]
        us[:, horizon_start:horizon_end-2] = us[:, horizon_start+1:horizon_end-1]

        # simulate for new last state for each horizon
        q_last = [discrete_dynamics(qs[1:n-1, horizon_end], us[1:m-1, horizon_end-1], dt);
                  simulate_path_dynamics(qs[end, horizon_end], us[horizon_end-1])]
        qs[:, horizon_end] = q_last

        # shift dynamics matrices for each horizon
        start_idx = horizon_start
        end_idx = horizon_start + N - k_c - 2 - (j - 1)

        vals.As[start_idx:end_idx-1] = vals.As[start_idx+1:end_idx]
        vals.Bs[start_idx:end_idx-1] = vals.Bs[start_idx+1:end_idx]
        vals.gs[start_idx:end_idx-1] = vals.gs[start_idx+1:end_idx]

        # fill in last entry of dynamics matrices for each horizon
        q_last = qs[:, horizon_end-2]
        u_last = us[:, horizon_end-2]
        vals.As[end_idx], vals.Bs[end_idx] = linearize_dynamics(discrete_dynamics, qs[1:n-1, horizon_end-1], us[1:m-1, horizon_end-1], dt)
        vals.gs[end_idx] = discrete_dynamics(qs[1:n-1, horizon_end-1], us[1:m-1, horizon_end-1], dt) - vals.As[end_idx] * qs[1:n-1, horizon_end-1] - vals.Bs[end_idx] * us[1:m-1, horizon_end-1]

        # shift tracking cost matrices for each horizon
        vals.cs[horizon_start:horizon_end-1] = vals.cs[horizon_start+1:horizon_end]
        vals.Γs[horizon_start:horizon_end-1] = vals.Γs[horizon_start+1:horizon_end]

        # correct second to last (stl) so that the weight is qc rather than qcterm
        P_stl = [qs[1, horizon_end-1], qs[2, horizon_end-1], qs[5, horizon_end-1]] #X, Y, s of last state of previous solution
        spline_idx = find_spline_interval(P_stl[3], path)
        c_stl = tracking_linear_term(P_stl, qc, ql, path, spline_idx)
        Γ_stl = tracking_quadratic_term(P_stl, qc, ql, path, spline_idx)
        vals.cs[horizon_end-1] = c_stl
        vals.Γs[horizon_end-1] = Γ_stl

        # update last entry
        P_last = [q_last[1], q_last[2], q_last[5]] #X, Y, s of new last state appended
        spline_idx = find_spline_interval(P_last[3], path)
        c_last = tracking_linear_term(P_last, qcterm, ql, path, spline_idx)
        Γ_last = tracking_quadratic_term(P_last, qcterm, ql, path, spline_idx)
        vals.cs[horizon_end] = c_last
        vals.Γs[horizon_end] = Γ_last
    end

    vals.c = vcat(vals.cs...)
    vals.Γ = Array(blockdiag(sparse.(vals.Γs)...))

    # update positions for obstacle constraints
    for k in 1:S_q
        vals.ps[k] = qs[1:2, k]
    end

    return qs, us
end

max_stars_repo_name: StanfordASL/MATS

max_issues_repo_name: StanfordASL/MATS

max_forks_repo_name: StanfordASL/MATS