text:
using RoME, Distributions

# [OPTIONAL] add more julia processes to speed up inference
using Distributed
nprocs() < 3 ? addprocs(4-nprocs()) : nothing

@everywhere using RoME

# number of particles
N = 100

# start with an empty factor graph object
fg = initfg()

# Add the first pose :x0
addVariable!(fg, :x0, Pose2, N=N)
# Add a few more poses
for i in 1:5
  addVariable!(fg, Symbol("x$(i)"), Pose2, N=N)
end
# Add landmarks
for i in 1:2
  addVariable!(fg, Symbol("l$(i)"), Point2, N=N)
end
# Add at a fixed location Prior to pin :x0 to a starting location (0,0,0)
addFactor!(fg, [:x0], PriorPose2( MvNormal([0; 0; 0], Matrix(Diagonal([0.3;0.3;0.3].^2)) )))

# Step 0
ppr = Pose2Point2Range(MvNormal([7.323699045815659], Matrix(Diagonal([0.3].^2))))
addFactor!(fg, [:x0; :l1], ppr )

# Step 1
i = 1
pp = Pose2Pose2(MvNormal([9.82039106655709;-0.040893643507463134;0.7851856362659602], Matrix(Diagonal([0.3;0.3;0.05].^2))))
ppr = Pose2Point2Range(MvNormal([6.783355788017825], Matrix(Diagonal([0.3].^2))))
addFactor!(fg, [Symbol("x$(i-1)"); Symbol("x$(i)")], pp )
addFactor!(fg, [Symbol("x$(i)"); Symbol("l1")], ppr )

# Step 2
i = i + 1
pp = Pose2Pose2(MvNormal([9.824301116341609;-0.1827443713759503;0.7586281983933238], Matrix(Diagonal([0.3;0.3;0.05].^2))))
ppr = Pose2Point2Range(MvNormal([6.768864478545091], Matrix(Diagonal([0.3].^2))))
addFactor!(fg, [Symbol("x$(i-1)"); Symbol("x$(i)")], pp )
addFactor!(fg, [Symbol("x$(i)"); Symbol("l2")], ppr )

# Step 3
i = i + 1
pp = Pose2Pose2(MvNormal([9.776502885351334; -0.010587078502017132; 1.5591793408311467], Matrix(Diagonal([0.3;0.3;0.05].^2))))
ppr = Pose2Point2Range(MvNormal([7.401417053438512], Matrix(Diagonal([0.3].^2))))
addFactor!(fg, [Symbol("x$(i-1)"); Symbol("x$(i)")], pp )
addFactor!(fg, [Symbol("x$(i)"); Symbol("l2")], ppr )

# Step 4
i = i + 1
pp = Pose2Pose2(MvNormal([9.644657137571507; 0.5847494762836476; 0.7422440549101994], Matrix(Diagonal([0.3;0.3;0.05].^2))))
ppr = Pose2Point2Range(MvNormal([7.331883435735532], Matrix(Diagonal([0.3].^2))))
addFactor!(fg, [Symbol("x$(i-1)"); Symbol("x$(i)")], pp )
addFactor!(fg, [Symbol("x$(i)"); Symbol("l2")], ppr )

# Step 5
i = i + 1
pp = Pose2Pose2(MvNormal([9.725096752125593; 0.6094800276622434; 0.7750400402527422], Matrix(Diagonal([0.3;0.3;0.05].^2))))
ppr = Pose2Point2Range(MvNormal([6.476782344345081], Matrix(Diagonal([0.3].^2))))
addFactor!(fg, [Symbol("x$(i-1)"); Symbol("x$(i)")], pp )
addFactor!(fg, [Symbol("x$(i)"); Symbol("l1")], ppr )

using Cairo, RoMEPlotting
Gadfly.set_default_plot_size(35cm, 20cm)

# draw initialized
plotSLAM2D(fg, contour=false)

for i in 1:10
# solve factor graph
tree,_,_ = solveTree!(fg, storeOld=true);

# after solve
plotSLAM2D(fg, contour=false)

pl1 = plotSLAM2DPoses(fg, to=0, drawhist=true);
X0 = IncrementalInference.getBelief(fg, :x0) |> getPoints;
# pl2 = Gadfly.layer(x=X0[1,:],y=X0[2,:], Geom.hexbin);
# you can also do
# push!(pl1.layers, pl2[1])
# pl1 # to show
# pl1 |> PDF("/tmp/test.pdf", 12cm,8cm)
  # plH = hstack(pl1, pl2)

# convert to file
p1H |> PNG("/home/chad/Pictures/julia_x0_step0.png")

#drawPoses(fg)
# If you have landmarks, you can instead call
plotSLAM2D(fg)

# # Draw the KDE for x0
# plotKDE(fg, :x0)
# # Draw the KDE's for x0 and x1
# plotKDE(fg, [:x1])
# plotKDE(fg, [:x2])
# plotKDE(fg, [:x3])
# plotKDE(fg, [:x4])
# plotKDE(fg, [:x5])
# plotKDE(fg, [:l1])
# plotKDE(fg, [:l2])

max_stars_repo_name: MarineRoboticsGroup/NF-iSAM

max_issues_repo_name: MarineRoboticsGroup/NF-iSAM

max_forks_repo_name: MarineRoboticsGroup/NF-iSAM


text:
function test_solve()
    prob = LightDark1D()
    sim = MCVISimulator()

    solver = MCVISolver(sim, nothing, 1, 100, 8, 500, 1000, 5000, 50, LightDark1DLowerBound(sim.rng), LightDark1DUpperBound(sim.rng))
    println("Solving...")
    policy = @time solve(solver, prob)
    println("...Solved")
    up = updater(policy)
    reward = simulate(sim, prob, policy, up, up.root)
    println("Reward:", reward)
    return true
end

function test_dummy_heuristics()
    prob = LightDark1D()
end

ext: jl


text:
using CUDAnative, CUDAdrv
using Base.Test

# NOTE: all kernel function definitions are prefixed with @eval to force toplevel definition,
#       avoiding boxing as seen in https://github.com/JuliaLang/julia/issues/18077#issuecomment-255215304

# NOTE: based on test/pkg.jl::capture_stdout, but doesn't discard exceptions
macro grab_output(ex)
    quote
        let fname = tempname()
            try
                ret = nothing
                open(fname, "w") do fout
                    redirect_stdout(fout) do
                        ret = $(esc(ex))
                    end
                end
                ret, readstring(fname)
            finally
                rm(fname, force=true)
            end
        end
    end
end

# Run some code on-device, returning captured standard output
macro on_device(exprs)
    @gensym kernel_fn
    quote
        let
            @eval function $kernel_fn()
                $(esc(exprs))

                return nothing
            end

            @cuda (1,1) $kernel_fn()
            synchronize()
        end
    end
end

@testset "CUDAnative" begin

@test devcount() > 0
if devcount() > 0
    include("base.jl")
    include("codegen.jl")

    # pick most recent device (based on compute capability)
    global dev = nothing
    for i in 0:devcount()-1
        newdev = CuDevice(i)
        if dev == nothing || capability(newdev) > capability(dev)
            dev = newdev
        end
    end
    info("Testing using device $(name(dev))")
    global ctx = CuContext(dev, CUDAdrv.SCHED_BLOCKING_SYNC)

    if capability(dev) < v"2.0"
        warn("native execution not supported on SM < 2.0")
    else
        include("execution.jl")
        include("array.jl")
        include("intrinsics.jl")

        include("examples.jl")
    end
end

end


text:
module RayCastings

export Point

struct Point{T}
    x::T
    y::T
end
Base.show(io::IO, p::Point) = print(io, "($(p.x), $(p.y))")

const Edge = Tuple{Point{T}, Point{T}} where T
Base.show(io::IO, e::Edge) = print(io, "$(e[1]) ∘-∘ $(e[2])")

function rayintersectseg(p::Point{T}, edge::Edge{T}) where T
    a, b = edge
    if a.y > b.y
        a, b = b, a
    end
    if p.y ∈ (a.y, b.y)
        p = Point(p.x, p.y + eps(p.y))
    end

    rst = false
    if (p.y > b.y || p.y < a.y) || (p.x > max(a.x, b.x))
        return false
    end

    if p.x < min(a.x, b.x)
        rst = true
    else
        mred = (b.y - a.y) / (b.x - a.x)
        mblu = (p.y - a.y) / (p.x - a.x)
        rst = mblu ≥ mred
    end

    return rst
end

isinside(poly::Vector{Tuple{Point{T}, Point{T}}}, p::Point{T}) where T =
    isodd(count(edge -> rayintersectseg(p, edge), poly))

connect(a::Point{T}, b::Point{T}...) where T =
    [(a, b) for (a, b) in zip(vcat(a, b...), vcat(b..., a))]

end  # module RayCastings

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@enum(PrimitiveType,
points,
lines,
linesStrip,
triangles,
trianglesStrip,
triangleFan,
quads)

immutable Vertex
    position::Vector2f
    color::Color
    texCoords::Vector2f
end

type VertexArray <: Drawable
    ptr::Ptr{Void}

    function VertexArray(ptr::Ptr{Void})
        v = new(ptr)
        finalizer(v, destroy)
        v
    end
end

function VertexArray()
    VertexArray(ccall((:sfVertexArray_create, libcsfml_graphics), Ptr{Void}, ()))
end

function copy(arr::VertexArray)
    return VertexArray(ccall((:sfVertexArray_copy, libcsfml_graphics), Ptr{Void}, (Ptr{Void},), arr.ptr))
end

function destroy(arr::VertexArray)
    ccall((:sfVertexArray_destroy, libcsfml_graphics), Void, (Ptr{Void},), arr.ptr)
end

function get_vertexcount(arr::VertexArray)
    return Int(ccall((:sfVertexArray_getVertexCount, libcsfml_graphics), UInt32, (Ptr{Void},), arr.ptr))
end

function get_vertex(arr::VertexArray, index::Integer)
    ptr = ccall((:sfVertexArray_getVertex, libcsfml_graphics), Ptr{Vertex}, (Ptr{Void}, Csize_t,), arr.ptr, index)
    return unsafe_load(ptr)
end

function clear(arr::VertexArray)
    ccall((:sfVertexArray_clear, libcsfml_graphics), Void, (Ptr{Void},), arr.ptr)
end

function resize(arr::VertexArray, vertex_count::Integer)
    ccall((:sfVertexArray_resize, libcsfml_graphics), Void, (Ptr{Void}, UInt32,), arr.ptr, vertex_count)
end

function append(arr::VertexArray, vertex::Vertex)
    ccall((:sfVertexArray_append, libcsfml_graphics), Void, (Ptr{Void}, Vertex,), arr.ptr, vertex)
end

function set_primitive_type(arr::VertexArray, primitive_type::PrimitiveType)
    ccall((:sfVertexArray_setPrimitiveType, libcsfml_graphics), Void, (Ptr{Void}, Int32,), arr.ptr, Int32(primitive_type))
end

function get_primitive_type(arr::VertexArray)
    return PrimitiveType(ccall((:sfVertexArray_getPrimitiveType, libcsfml_graphics), Int32, (Ptr{Void},), arr.ptr))
end

function getbounds(arr::VertexArray)
    return ccall((:sfVertexArray_getBounds, libcsfml_graphics), FloatRect, (Ptr{Void},), arr.ptr)
end


text:
export Particle, Particles

struct Particle

    x :: Float64
    v :: ComplexF64
    e :: ComplexF64
    b :: Float64

end

mutable struct Particles

    nbpart :: Int64

    x :: Array{Float64,2}
    v :: Array{Float64,2}
    e :: Array{Float64,2}
    b :: Vector{Float64}
    t :: Vector{Float64}
    w :: Float64

    function Particles( nbpart :: Int64, w :: Float64 )

        x = zeros(Float64, (2,nbpart))
        v = zeros(Float64, (2,nbpart))
        e = zeros(Float64, (2,nbpart))
        b = zeros(Float64, nbpart)
        t = zeros(Float64, nbpart)

        new( nbpart, x, v, e, b, t , w )

    end

end

include("read_particles.jl")
include("plasma.jl")
include("compute_rho.jl")
include("interpolation.jl")


text:
module JSONConverter

export convert_JSON

using ..FlexPlan
const _FP = FlexPlan
import ..FlexPlan: _LOGGER

import JSON
import Memento

include("base.jl")
include("nw.jl")
include("time_series.jl")
include("singular_data.jl")

end


text:
@testset "AlgAssAbsOrd" begin
  Qx, x = FlintQQ["x"]

  @testset "Quaternion Algebras" begin
    function sum_of_two_squares(a::fmpz)
      for i=1:Int(root(a,2))
        for j=1:Int(root(a,2))
          if a==i^2+j^2
            return true
          end
        end
      end
      return false
    end

    for b in Hecke.squarefree_up_to(100)[2:end]
      K, a = NumberField(x^2-b, check = false, cached = false)
      O = maximal_order(K);
      cocval = Array{nf_elem, 2}(undef, 2, 2)
      G = NfToNfMor[Hecke.NfToNfMor(K,K,a),Hecke.NfToNfMor(K,K,-a)]
      cocval[1,1] = K(1)
      cocval[1,2] = K(1)
      cocval[2,1] = K(1)
      cocval[2,2] = K(-1)
      A = Hecke.CrossedProductAlgebra(K,G,cocval)
      if Hecke.issplit(A)
        A1 = Hecke.CrossedProductAlgebra(O, G, cocval)
        O1 = Order(A1, basis(A1))
        d = discriminant(O1)
        fac = factor(d)
        for p in keys(fac.fac)
          On = Hecke.pmaximal_overorder(O1, Int(p))
          @test valuation(discriminant(On), p) == 0
        end
        @test sum_of_two_squares(fmpz(b))
      else
        @test !sum_of_two_squares(fmpz(b))
      end
    end

    A = Hecke.quaternion_algebra2(4,36)
    @test Hecke.issplit(A)
    A=Hecke.quaternion_algebra2(-1,-1)
    O= Order(A, [A[i] for i=1:4])
    @test Hecke.schur_index_at_real_plc(O)==2

  end

  @testset "Crossed Product Order" begin

    K, a = NumberField(x^4-4*x^2+1)
    O = maximal_order(K)
    Autos = Array{NfToNfMor, 1}(undef, 4)
    Autos[1] = NfToNfMor(K, K, a)
    Autos[2] = NfToNfMor(K, K, -a)
    Autos[3] = NfToNfMor(K, K, a^3 - 4*a)
    Autos[4] = NfToNfMor(K, K, -a^3 + 4*a)
    MatCoc = [0 0 0 0; 0 1 0 1; 0 1 1 0; 0 0 1 1]
    Coc = Array{nf_elem, 2}(undef, 4, 4)
    for i = 1:4
      for j = 1:4
        Coc[i, j] = K(-1)^MatCoc[i, j]
      end
    end
    A = Hecke.CrossedProductAlgebra(O, Autos, Coc)
    @test Hecke.issplit(A)
    O1 = Order(A, basis(A))
    d = discriminant(O1)
    fac1 = factor(discriminant(O))
    fac2 = factor(d)
    @test Set(collect(keys(fac1.fac))) == Set(collect(keys(fac2.fac)))
    for p in keys(fac1.fac)
      O3 = Hecke.pmaximal_overorder(O1, Int(p))
      @test valuation(discriminant(O3), p) == 0
    end
  end

  @testset "Any order" begin
    A = AlgAss(x^2 - fmpq(1, 5))

    O = any_order(A)

    # Test whether O is indeed an order
    @test one(A) in O

    for b in basis(O, copy = false)
      @test isintegral(elem_in_algebra(b, copy = false))
      for c in basis(O, copy = false)
        @test elem_in_algebra(b*c, copy = false) in O
      end
    end

    # Some more miscellaneous tests
    b = rand(O, 10)
    @test elem_in_algebra(b) in O

    b = A(fmpq[ fmpq(1), fmpq(1, 3) ])
    @test denominator(b, O)*b in O
  end

  @testset "Maximal Order" begin
    A = AlgAss(x^2 + 10x + 7)
    AA = deepcopy(A) # avoid caching
    O1 = maximal_order(A)
    O2 = Hecke.maximal_order_via_decomposition(AA)
    @test discriminant(O1) == discriminant(O2)

    QG = group_algebra(FlintQQ, small_group(2, 1))
    O = maximal_order(QG)
    @test isone(abs(discriminant(O)))

    M = matrix_algebra(FlintQQ, 3)
    O = maximal_order(M)
    @test isone(abs(discriminant(O)))
  end
end


text:
# Autogenerated wrapper script for SCALAPACK32_jll for x86_64-apple-darwin-libgfortran3
export libscalapack32

using OpenBLAS32_jll
using MPICH_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("SCALAPACK32")
JLLWrappers.@declare_library_product(libscalapack32, "@rpath/libscalapack32.dylib")
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS32_jll, MPICH_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libscalapack32,
        "lib/libscalapack32.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export HyperPlanesIntersection, LinesHull, affinehull, linespace

# It is easy to go from H-rep of affine space to V-rep of affine space by computing the kernel of a matrix using RowEchelon
# However, it is really worth it since Base.in for an HRepElement in HyperPlanesIntersection and Base.in for an VRepElement in LinesHull are false already.
function remproj(x::RepElement{N, <:Integer}, l::RepElement{N, <:Integer}) where N
    # generates large numbers but keeps the type integer
    x * dot(coord(l), coord(l)) - l * dot(coord(x), coord(l))
end
function remproj(x, l)
    ll = dot(coord(l), coord(l))
    if isapproxzero(ll)
        x
    else
        simplify(x - l * ((dot(coord(x), coord(l)) / ll)))
    end
end

# H-representation

"""
    hrep(hyperplanes::HyperPlaneIt)

Creates an affine space from the list of hyperplanes `hyperplanes`.

### Examples
```julia
hrep([HyperPlane([0, 1, 0], 1), HyperPlane([0, 0, 1], 0)])
```
creates the 1-dimensional affine subspace containing all the points ``(x_1, 0, 0)``, i.e. the ``x_1``-axis.

```julia
hrep([HyperPlane([1, 1], 1), HyperPlane([1, 0], 0)])
```
creates the 0-dimensional affine subspace only containing the point ``(0, 1)``.
"""
hrep(hyperplanes::HyperPlaneIt) = HyperPlanesIntersection(hyperplanes)

# Representation of an affine space as the intersection of hyperplanes.
# Also called affine set, affine manifold, affine variety, linear variety or flat.
# An affine space L satisfies:
# λx + (1-λ)y ∈ L, ∀x, y ∈ L, ∀ λ ∈ R
# Note that λ is not rhyperplaneuired to be between 0 and 1 as in convex sets.
struct HyperPlanesIntersection{N, T, AT} <: HAffineSpace{N, T}
    # HyperPlanes whose intersection is the affine space
    hyperplanes::Vector{HyperPlane{N, T, AT}}
    function HyperPlanesIntersection{N, T, AT}(hps::HyperPlaneIt{N, T}) where {N, T, AT}
        new{N, T, AT}(lazy_collect(hps))
    end
end
hvectortype(L::Type{HyperPlanesIntersection{N, T, AT}}) where {N, T, AT} = AT
similar_type(PT::Type{<:HyperPlanesIntersection}, d::FullDim{N}, ::Type{T}) where {N, T} = HyperPlanesIntersection{N, T, similar_type(hvectortype(PT), d, T)}

HyperPlanesIntersection{N, T, AT}() where {N, T, AT} = HyperPlanesIntersection{N, T, AT}(HyperPlane{N, T, AT}[])
HyperPlanesIntersection(it::ElemIt{HyperPlane{N, T, AT}}) where {N, T, AT} = HyperPlanesIntersection{N, T, AT}(it)

Base.intersect!(L::HyperPlanesIntersection{N}, h::HyperPlane{N}) where N = push!(L.hyperplanes, h)

@vecrepelem HyperPlanesIntersection HyperPlane hyperplanes

hreptype(::Type{<:HyperPlanesIntersection{N, T, AT}}) where {N, T, AT} = Intersection{N, T, AT}

# Returns an HyperPlanesIntersection representing the affine hull of p.
# The affine hull is defined as
# {λx + (1-λ)y | x, y ∈ p, λ ∈ R}
function affinehull(h::HRep, current=false)
    if !current
        detecthlinearity!(h)
    end
    HyperPlanesIntersection(hyperplanes(h))
end

function remproj(h::HRepElement, L::HyperPlanesIntersection)
    for hp in hyperplanes(L)
        h = remproj(h, hp)
    end
    h
end
function Base.in(h::HRepElement, L::HyperPlanesIntersection)
    h = remproj(h, L)
    isapproxzero(h)
end

function removeduplicates(L::HyperPlanesIntersection{N, T, AT}) where {N, T, AT}
    H = HyperPlanesIntersection{N, T, AT}()
    for h in hyperplanes(L)
        if !(h in H)
            intersect!(H, h)
        end
    end
    H
end

# V-representation
struct VEmptySpace{N, T, AT} <: VLinearSpace{N, T} end
emptyspace(v::VRep{N, T}) where {N, T} = VEmptySpace{N, T, vvectortype(typeof(v))}()

Base.in(v::VRepElement, L::VEmptySpace) = isapproxzero(v)

"""
    vrep(lines::LineIt)

Creates an affine space from the list of lines `lines`.

### Examples
```julia
vrep([Line([1, 0, 0]), Line([0, 1, 0])])
```
creates the 2-dimensional affine subspace containing all the points ``(x_1, x_2, 0)``, i.e. the ``x_1````x_2``-plane.
"""
vrep(lines::LineIt) = LinesHull(lines)

# Representation of an affine space containing the origin by the minkowsky sum of lines
struct LinesHull{N, T, AT} <: VLinearSpace{N, T}
    lines::Vector{Line{N, T, AT}}
    function LinesHull{N, T, AT}(lines::LineIt{N, T}) where {N, T, AT}
        new{N, T, AT}(lazy_collect(lines))
    end
end
vvectortype(::Type{LinesHull{N, T, AT}}) where {N, T, AT} = AT
similar_type(PT::Type{<:LinesHull}, d::FullDim{N}, ::Type{T}) where {N, T} = LinesHull{N, T, similar_type(vvectortype(PT), d, T)}

LinesHull{N, T, AT}() where {N, T, AT} = LinesHull{N, T, AT}(Line{N, T, AT}[])
LinesHull(it::ElemIt{Line{N, T, AT}}) where {N, T, AT} = LinesHull{N, T, AT}(it)

convexhull!(L::LinesHull{N}, l::Line{N}) where {N} = push!(L.lines, l)

@vecrepelem LinesHull Line lines

conetype(::Type{LinesHull{N, T, AT}}) where {N, T, AT} = RaysHull{N, T, AT}
vreptype(::Type{LinesHull{N, T, AT}}) where {N, T, AT} = Hull{N, T, AT}

# Returns a LinesHull representing the following set (TODO does it have a name?, does someone has a reference talking about it ?)
# {x | ⟨a, x⟩ = 0 ∀ a such that (α, β) is a valid hyperplane for p}
function linespace(v::VRep, current=false)
    if !current
        detectvlinearity!(v)
    end
    LinesHull(lines(v))
end

function remproj(v::VRepElement, L::LinesHull)
    for l in lines(L)
        v = remproj(v, l)
    end
    v
end
function Base.in(v::VRepElement, L::LinesHull)
    v = remproj(v, L)
    isapproxzero(v)
end

function removeduplicates(L::LinesHull{N, T, AT}) where {N, T, AT}
    V = LinesHull{N, T, AT}()
    for l in lines(L)
        if !(l in H)
            convexhull!(H, h)
        end
    end
    H
end


text:
function movements!(r,side::HorizonSide)
    while isborder(r,side) ==  false
        move!(r,side)
        putmarker!(r)
    end
end

function movements_mark(r,side::HorizonSide)
    while ismarker(r) == true
        move!(r,side)
    end
end

function krest(r)

    for side in [Nord, Ost, Sud, West]
        movements!(r,side)
        movements_mark(r,HorizonSide((Int(side)+2)%4))
    end

    putmarker!(r)
    show(r)
end

krest(r)
show(r)

max_stars_repo_name: AntonDudDm/mirea-progs

max_issues_repo_name: AntonDudDm/mirea-progs

max_forks_repo_name: AntonDudDm/mirea-progs


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

# grid neighborhood for gradients
function gridneighbors(img, i::CartesianIndex, window::Int)
    dim = Size(img)
    minid(d) = max(1,i[d]-window)
    maxid(d) = min(dim[d],i[d]+window)
    idx = Tuple([minid(d):maxid(d) for d in 1:length(dim)])
    CartesianIndices(idx)
end

# nearest grid data id to some other sample data
function grid2hd_ids(pdata,pdomain)
  hd   = [coordinates(centroid(pdata, i)) for i in 1:nelements(pdata)]
  grid = [coordinates(centroid(pdomain, i)) for i in 1:nelements(pdomain)]

  tree = KDTree(grid)
  idxs, dists = nn(tree, hd)

  [i for i in idxs]
end

# nearest grid data mahalanobis matrix to some other sample data
function grid2hd_qmat(pdata,pdomain,localaniso)
  idxs = grid2hd_ids(pdata,pdomain)
  [qmat(rotation(localaniso,i),magnitude(localaniso,i)) for i in idxs]
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -3.892086739821988 <= q <= 1.133112254487742)

                     add_NL_constraint(m, :(log(1 + exp(-0.9215225812660202 + -0.8659264200816277*log(1 + exp(0.9754315885766114 + 0.9182647481159649*$(x[1]) + 0.2963572993024828*$(x[2]))) + -0.7655479829179535*log(1 + exp(0.7516798050016731 + -0.4187259277447679*$(x[1]) + -0.991762240929563*$(x[2]))))) + log(1 + exp(0.9152212659843197 + -0.09622429015475253*log(1 + exp(0.9754315885766114 + 0.9182647481159649*$(x[1]) + 0.2963572993024828*$(x[2]))) + 0.18727818139193309*log(1 + exp(0.7516798050016731 + -0.4187259277447679*$(x[1]) + -0.991762240929563*$(x[2]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Underscores
using StaticArrays

const Num = Int32

mutable struct Circle
    start::Int
    list::Array{Num,1}
end

function Circle(list)
    re = Array{Num}(undef, length(list))
    for i in 2:length(list)
        re[list[i-1]] = list[i]
    end
    re[list[end]] = list[1]
    Circle(list[1], re)
end

function parse_file(f)::Circle
    @_ collect(f) |> map(parse(Num, _), __) |> Circle
end

Base.length(circle::Circle) = length(circle.list)

Base.copy(circle::Circle) = Circle(circle.start, copy(circle.list))

Base.:(==)(lhs::Circle, rhs::Circle) = lhs.start == rhs.start && lhs.list == rhs.list

const REMOVE = 3

function extract_after(circle::Circle)
    re = MVector{REMOVE,Num}(0, 0, 0)
    which = circle.start
    for i in 1:REMOVE
        which = circle.list[which]
        re[i] = which
    end
    re
end

function insert_after!(into::Circle, what, after)
    into.list[into.start] = into.list[what[end]]
    into.list[what[end]] = into.list[after]
    into.list[after] = what[1]
end

function step!(circle::Circle)
    pickup = extract_after(circle)
    next = circle.list[pickup[end]]
    searchfor = mod1(circle.start - 1, length(circle))
    while searchfor in pickup
        searchfor = mod1(searchfor - 1, length(circle))
    end
    insert_after!(circle, pickup, searchfor)
    circle.start = next
end

function ex1(circle::Circle, rounds=10)
    circle = copy(circle)
    for i in 1:rounds
        step!(circle)
    end
    re = []
    itr = circle.list[1]
    while isempty(re) || itr != 1
        push!(re, itr)
        itr = circle.list[itr]
    end
    join(re)
end

const MILLION = 1_000_000

function ex2(problem::Array{Num,1}, rounds=10_000_000)
    big = collect(1:MILLION)
    big[1:length(problem)] .= problem
    circle = Circle(big)
    for i in 1:rounds
        step!(circle)
    end
    fst = circle.list[1]
    snd = circle.list[fst]
    widemul(fst, snd)
end

problem = parse_file("364289715")
println("Result of game one: ", ex1(problem))
println("Product of the huge game: ", ex2(Num[3,6,4,2,8,9,7,1,5]))

using Test

@testset "Crab Cups" begin
    example = [3,8,9,1,2,5,4,6,7]
    example_circle = Circle(example)

    @testset "circle" begin
        @test example_circle.start == 3
        @test example_circle.list[example_circle.start] == 8

        @test extract_after(example_circle).data == (8,9,1)

        @testset "insert_after!" begin
            to_change = copy(example_circle)
            after_step1 = Circle([2,8,9,1,5,4,6,7,3]) # 3 (2) 8  9  1  5  4  6  7
            after_step2 = Circle([5,4,6,7,8,9,1,3,2]) # 3  2 (5) 4  6  7  8  9  1
            after_step3 = Circle([8,9,1,3,4,6,7,2,5]) # 7  2  5 (8) 9  1  3  4  6
            after_step4 = Circle([4,6,7,9,1,3,2,5,8]) # 3  2  5  8 (4) 6  7  9  1

            set = to_change.list[1]
            insert_after!(to_change, [8,9,1], 2)
            to_change.start = set
            @test to_change == after_step1

            set = to_change.list[1]
            insert_after!(to_change, [8,9,1], 7)
            to_change.start = set
            @test to_change == after_step2

            set = to_change.list[7]
            insert_after!(to_change, [4,6,7], 3)
            to_change.start = set
            @test to_change == after_step3

            set = to_change.list[3]
            insert_after!(to_change, [9,1,3], 7)
            to_change.start = set
            @test to_change == after_step4
        end
    end

    @testset "example 1" begin
        example = parse_file("389125467")

        @test ex1(example, 10) == "92658374"
    end

    @testset "results" begin
        @test ex1(problem, 100) == "98645732"
        @test ex2(Num[3,6,4,2,8,9,7,1,5]) == 689500518476
    end
end

max_stars_repo_name: CmdQ/AoC


text:
using LinearAlgebra
using SparseArrays
using KernelAbstractions
using BenchmarkTools

using ExaPF
import ExaPF: LinearSolvers

const LS = LinearSolvers

function build_nlp(datafile, device)
    print("Load data\t")
    polar = @time PolarForm(datafile, device)

    constraints = Function[
        ExaPF.voltage_magnitude_constraints,
        ExaPF.active_power_constraints,
        ExaPF.reactive_power_constraints,
    ]
    print("Constructor\t")
    powerflow_solver = NewtonRaphson(tol=1e-10)
    nlp = @time ExaPF.ReducedSpaceEvaluator(polar; constraints=constraints,
                                            powerflow_solver=powerflow_solver)
    return nlp
end

function run_reduced_evaluator(nlp, u; device=CPU())
    # Update nlp to stay on manifold
    print("Update   \t")
    @time ExaPF.update!(nlp, u)
    # Compute objective
    print("Objective\t")
    c = @btime ExaPF.objective($nlp, $u)
    # Compute gradient of objective
    g = similar(u)
    fill!(g, 0)
    print("Gradient \t")
    @btime ExaPF.gradient!($nlp, $g, $u)

    # Constraint
    ## Evaluation of the constraints
    cons = similar(nlp.g_min)
    fill!(cons, 0)
    print("Constrt \t")
    @btime ExaPF.constraint!($nlp, $cons, $u)

    print("Jac-prod \t")
    jv = copy(g) ; fill!(jv, 0)
    v = copy(cons) ; fill!(v, 1)
    @btime ExaPF.jtprod!($nlp, $jv, $u, $v)
    hv = similar(u) ; fill!(hv, 0)
    v = similar(u) ; fill!(v, 0)
    v[1] = 1
    print("Hessprod \t")
    @btime ExaPF.hessprod!($nlp, $hv, $u, $v)
    y = similar(cons) ; fill!(y, 1.0)
    w = similar(cons) ; fill!(w, 1.0)
    print("HLagPen-prod \t")
    @btime ExaPF.hessian_lagrangian_penalty_prod!($nlp, $hv, $u, $y, 1.0, $v, $w)
    return
end

function run_slack(nlp, u; device=CPU())
    slack = ExaPF.SlackEvaluator(nlp)
    w = ExaPF.initial(slack)
    print("Update   \t")
    @time ExaPF.update!(slack, w)
    print("Objective\t")
    c = @time ExaPF.objective(slack, w)
    g = similar(w)
    fill!(g, 0)
    print("Gradient \t")
    @time ExaPF.gradient!(slack, g, w)
    cons = similar(nlp.g_min)
    fill!(cons, 0)
    print("Constrt \t")
    @time ExaPF.constraint!(slack, cons, w)

    hv = similar(w) ; fill!(hv, 0)
    v = similar(w) ;  fill!(v, 1)
    print("Hessprod \t")
    @time ExaPF.hessprod!(slack, hv, w, v)
    y = similar(cons) ; fill!(y, 1.0)
    z = similar(cons) ; fill!(w, 1.0)
    print("HLagPen-prod \t")
    @time ExaPF.hessian_lagrangian_penalty_prod!(slack, hv, w, y, 1.0, v, z)
    return
end

function run_penalty(nlp, u; device=CPU())
    pen = ExaPF.AugLagEvaluator(nlp, u)
    print("Update   \t")
    @time ExaPF.update!(pen, u)
    print("Objective\t")
    c = @time ExaPF.objective(pen, u)
    g = similar(u)
    fill!(g, 0)
    print("Gradient \t")
    @time ExaPF.gradient!(pen, g, u)
    hv = similar(u) ; fill!(hv, 0)
    v = similar(u) ;  fill!(v, 1)
    print("Hessprod \t")
    @time ExaPF.hessprod!(pen, hv, u, v)
    return
end

function run_slackaug(nlp, u; device=CPU())
    nlp2 = ExaPF.SlackEvaluator(nlp)
    w = ExaPF.initial(nlp2)
    pen = ExaPF.AugLagEvaluator(nlp2, w)
    print("Update   \t")
    @time ExaPF.update!(pen, w)
    print("Objective\t")
    c = @time ExaPF.objective(pen, w)
    g = similar(w)
    fill!(g, 0)
    print("Gradient \t")
    @time ExaPF.gradient!(pen, g, w)
    hv = similar(w) ; fill!(hv, 0)
    v = similar(w) ;  fill!(v, 1)
    print("Hessprod \t")
    @time ExaPF.hessprod!(pen, hv, w, v)
    return
end

function run_proxal(nlp, u; device=CPU())
    pen = ExaPF.ProxALEvaluator(nlp, ExaPF.Normal)
    w = ExaPF.initial(pen)
    print("Update   \t")
    @time ExaPF.update!(pen, w)
    print("Objective\t")
    c = @time ExaPF.objective(pen, w)
    g = similar(w)
    fill!(g, 0)
    print("Gradient \t")
    @time ExaPF.gradient!(pen, g, w)
    return
end

function run_proxaug(nlp, u; device=CPU())
    pen = ExaPF.ProxALEvaluator(nlp, ExaPF.Normal)
    w = ExaPF.initial(pen)
    aug = ExaPF.AugLagEvaluator(pen, w)
    print("Update   \t")
    @time ExaPF.update!(aug, w)
    print("Objective\t")
    c = @time ExaPF.objective(aug, w)
    g = similar(w)
    fill!(g, 0)
    print("Gradient \t")
    @time ExaPF.gradient!(aug, g, w)
    return
end

datafile = joinpath(dirname(@__FILE__), "..", "data", "case300.m")
device = CPU()
nlp = build_nlp(datafile, CPU())
u = ExaPF.initial(nlp)

@info("ReducedSpaceEvaluator")
run_reduced_evaluator(nlp, u, device=CPU())
ExaPF.reset!(nlp)
# @info("SlackEvaluator")
# run_slack(nlp, u, device=CPU())
# ExaPF.reset!(nlp)
# @info("AugLagEvaluator")
# run_penalty(nlp, u, device=CPU())
# ExaPF.reset!(nlp)
# @info("SlackAugLagEvaluator")
# run_slackaug(nlp, u, device=CPU())
# ExaPF.reset!(nlp)


text:
using triangle_area
using Base.Test

# write your own tests here
area=heron(3,4,5)
println("변의 길이가 3, 4, 5인 삼각형의 넓이는 ",area)
@test area==6


text:
using Test
using UArrayInts

@testset "array funs" begin

    tmp1 = [0xff, 0xfe]
    tmp2 = UArrayInt(tmp1)
    @test tmp1[1] == tmp2[1]
    @test tmp1[end] == tmp2[end]
    @test length(tmp1) == length(tmp2)

    tmp2[1] = 0x00
    @test tmp2[1] == 0x00

end

@testset "add_carry" begin
    function test_add_carry(x::T, y::T) where T
        tmp = add_carry(x, y)
        tmp[1] + (T(1) << (sizeof(T) * 8)) * tmp[2] == x + y
    end
    @test test_add_carry(0xff, 0xff)
    @test test_add_carry(0xffff, 0xffff)
    @test test_add_carry(0xffffffff, 0xffffffff)
    @test test_add_carry(0xffffffffffffffff, 0xffffffffffffffff)
    @test test_add_carry(0xffffffffffffffffffffffffffffffff, 0xffffffffffffffffffffffffffffffff)
    @test test_add_carry(0x0f, 0x0f)
    @test test_add_carry(0x000f, 0x000f)
    @test test_add_carry(0x0000000f, 0x0000000f)
    @test test_add_carry(0x000000000000000f, 0x000000000000000f)
    @test test_add_carry(0x0000000000000000000000000000000f, 0x0000000000000000000000000000000f)
    @test UArrayInts.to_unsigned(UArrayInts.add!(([0xff, 0xff]), ([0xff, 0xff]), 1)) == 0xffff + 0xff00
end

@testset "add" begin
    function inner_test_add(x, y)
        x2 = deepcopy(x)
        add!(x, y)
        tmp1 = to_unsigned(x)
        tmp2 = (to_unsigned(x2) + to_unsigned(y)) % typeof(tmp1)
        tmp1 == tmp2
    end
    function test_add(::Type{T}) where T
        a = UArrayInt
        z = zero(T)
        o = one(T)
        m = typemax(T)

        @test inner_test_add([m], [m])
        sizeof(T) < 16 && @test inner_test_add([m, m], [m])
        sizeof(T) < 8  && @test inner_test_add([m, m], [m, m, m])
        @test inner_test_add([o], [o])
        sizeof(T) < 16 && @test inner_test_add([z, o], [o])
        sizeof(T) < 8  && @test inner_test_add([z, o], [z, z, o])

        @test inner_test_add(a([m]), a([m]))
        sizeof(T) < 16 && @test inner_test_add(a([m, m]), a([m]))
        sizeof(T) < 8  && @test inner_test_add(a([m, m]), a([m, m, m]))
        @test inner_test_add(a([o]), a([o]))
        sizeof(T) < 16 && @test inner_test_add(a([z, o]), a([o]))
        sizeof(T) < 8  && @test inner_test_add(a([z, o]), a([z, z, o]))

        return nothing
    end

    test_add(UInt8)
    test_add(UInt16)
    test_add(UInt32)
    test_add(UInt64)

end

@testset "mul_carry" begin function test_mul_carry(x::T, y::T) where T
        tmp = mul_carry(x, y)
        tmp[1] + (T(1) << (sizeof(T) * 8)) * tmp[2] == x * y
    end
    @test test_mul_carry(0xff, 0xff)
    @test test_mul_carry(0xffff, 0xffff)
    @test test_mul_carry(0xffffffff, 0xffffffff)
    @test test_mul_carry(0xffffffffffffffff, 0xffffffffffffffff)
    @test_throws ErrorException test_mul_carry(0xffffffffffffffffffffffffffffffff, 0xffffffffffffffffffffffffffffffff)
    @test test_mul_carry(0x0f, 0x0f)
    @test test_mul_carry(0x000f, 0x000f)
    @test test_mul_carry(0x0000000f, 0x0000000f)
    @test test_mul_carry(0x000000000000000f, 0x000000000000000f)
    @test_throws ErrorException test_mul_carry(0x0000000000000000000000000000000f, 0x0000000000000000000000000000000f)
end

@testset "mul" begin
    function inner_test_mul(x, y)
        x2 = deepcopy(x)
        mul!(x2, x, y)
        tmp1 = to_unsigned(x2)
        tmp2 = (to_unsigned(x) * to_unsigned(y)) % typeof(tmp1)
        tmp1 == tmp2
    end
    function test_mul(::Type{T}) where T
        a = UArrayInt
        z = zero(T)
        o = one(T)
        m = typemax(T)

        @test inner_test_mul([m], [m])
        sizeof(T) < 16 && @test inner_test_mul([m, m], [m])
        sizeof(T) < 8  && @test inner_test_mul([m, m], [m, m, m])
        @test inner_test_mul([o], [o])
        sizeof(T) < 16 && @test inner_test_mul([z, o], [o])
        sizeof(T) < 8  && @test inner_test_mul([z, o], [z, z, o])

        @test inner_test_mul(a([m]), a([m]))
        sizeof(T) < 16 && @test inner_test_mul(a([m, m]), a([m]))
        sizeof(T) < 8  && @test inner_test_mul(a([m, m]), a([m, m, m]))
        @test inner_test_mul(a([o]), a([o]))
        sizeof(T) < 16 && @test inner_test_mul(a([z, o]), a([o]))
        sizeof(T) < 8  && @test inner_test_mul(a([z, o]), a([z, z, o]))

        return nothing
    end

    test_mul(UInt8)
    test_mul(UInt16)
    test_mul(UInt32)
    test_mul(UInt64)

end

@testset "to_unsigned" begin
    @test to_unsigned([0x04, 0x03, 0x02, 0x01], Val(4)) == 0x04030201
    @test to_unsigned([0x0004, 0x0003, 0x0002, 0x0001], Val(4)) == 0x0004000300020001
    tmp1 = to_unsigned([0x00000004, 0x00000003, 0x00000002, 0x00000001], Val(4))
    tmp2 = 0x00000004000000030000000200000001
    @test tmp1 == tmp2
    @test to_unsigned([0x03, 0x02, 0x01], Val(3)) == 0x00030201
end


text:
## Main Finite Difference Model ##
type FiniteDifferenceModel{T}
  evolver::T
  stoppingTimes::Vector{Float64}

  FiniteDifferenceModel{T}(evolver::T, stoppingTimes::Vector{Float64}) = new(evolver, sort(unique(stoppingTimes)))
end

function FiniteDifferenceModel(L::TridiagonalOperator, BCs::Vector{BoundaryCondition}, evolverFunc::Function)
  evolver = evolverFunc(L, FdmBoundaryConditionSet(BCs))
  return FiniteDifferenceModel{typeof(evolver)}(evolver, Vector{Float64}())
end

function FiniteDifferenceModel(L::Vector{TridiagonalOperator}, BCs::Matrix{BoundaryCondition}, evolverFunc::Function)
  evolver = ParallelEvolver(L, BCs, evolverFunc)
  return FiniteDifferenceModel{ParallelEvolver}(evolver, Vector{Float64}())
end

function rollback_impl!(model::FiniteDifferenceModel, a::Vector, from::Float64, to::Float64, steps::Int, condition::StepCondition)
  dt = (from - to) / steps
  t = from
  set_step!(model.evolver, dt)

  if length(model.stoppingTimes) > 0 && model.stoppingTimes[end] == from
    apply_to!(condition, a, from)
  end

  for i = 1:steps
    _now = t
    _next = t - dt
    hit = false
    for j = length(model.stoppingTimes):-1:1
      if _next <= model.stoppingTimes[j] && model.stoppingTimes[j] < _now
        # a stopping time was hit
        hit = true

        # perform a small step to stoppingTimes[j]
        set_step!(model.evolver, _now - model.stoppingTimes[j])
        step!(model.evolver, a, _now)

        apply_to!(condition, a, model.stoppingTimes[j])

        # and continue the cycle
        _now = model.stoppingTimes[j]
      end
    end

    # if we did hit
    if hit
      # we might have to make a small step to complete the big one
      if _now > _next
        set_step!(model.evolver, _now - _next)
        step!(model.evolver, a, _now)
        apply_to!(condition, a, _next)
      end

      # and in any case, we have to reset the evolver to the default step
      set_step!(model.evolver, dt)
    else
      # if we didn't, the evolver is already set to the default step, which is ok
      step!(model.evolver, a, _now)
      apply_to!(condition, a, _next)
    end
    t -= dt
  end

  return model, a
end

rollback!(model::FiniteDifferenceModel, a::Vector, from::Float64, to::Float64, steps::Int, condition::StepCondition = FdmNullStepCondition()) = rollback_impl!(model, a, from, to, steps, condition)


text:
x = linspace(0, 0.2, 10)
y = zeros(x)
z = zeros(x)
r = [x'; y'; z']
a = randn(10).^2  # squared to make sure the're all positive

# Different types to check promotion
h = 1.02 * ones(Float32, length(a))
g = 1.04 * ones(Float16, length(a))

weird_zoop = Scatterer(r, a, h, g)
werid_zoop = Scatterer(r, a, 1.02, 1.04)
werid_zoop = Scatterer(r, a, 1.02, 1)

# Make sure all models load
krill1 = Models.krill_mcgeehee
krill2 = Models.krill_conti
cope = Models.calanoid_copepod
sandeel = Models.sandeel
daphnia = Models.daphnia

@test_approx_eq(length(rotate(krill1, roll=90, tilt=30, yaw=12)), length(krill1))
@assert length(resize(krill1, 0.03)) == 0.03


text:
module JavisGraphs

import Base
#=
This funny mix of using and imports is required due to a coding style followed in JavisGraphs.
* When 'using' something from an external package, never make use of qualifiers like `Javis.` or `Graph.` for simplicity and to make sure there is no naming conflict introduced by JavisGraphs itself. Thus method calls in JavisGraphs.jl are distinguished by argument lists solely and not namespaces. (As of this date, Javis.jl and Graphs.jl are single moduled packages i.e. no intra-naming conflicts exists, and they have no inter-naming conflicts either.)
- Hence unexported types/functions like Javis.get_frames, Javis.Frames has to be used with `using Javis: Frames, ...`
- Rest of the exported names can be brought in simply with `using Javis`
* Only when writing methods for exported Javis/Graphs functions use the namespace qualifiers for clarity.
* ref - https://docs.julialang.org/en/v1/manual/modules/#using-and-import-with-specific-identifiers,-and-adding-methods
=#
using Javis: isempty, svgwh, get_frames, get_latex_svg, Frames, AbstractObject, LaTeXString
using Javis
import Javis: text, setopacity, setline, setdash, setfont, sethue
using Graphs
import Graphs: weights, add_vertex!, rem_vertex!, add_edge!, rem_edge!
using LaTeXStrings
using GraphPlot
using LinearAlgebra
using OrderedCollections
using SparseArrays: SparseMatrixCSC

abstract type AbstractJavisGraphElement end
abstract type AbstractJavisGraph <: AbstractJavisGraphElement end
abstract type AbstractGraphElement <: AbstractJavisGraphElement end
abstract type AbstractGraphVertex <: AbstractGraphElement end
abstract type AbstractGraphEdge <: AbstractGraphElement end

include("utils.jl")

include("structs/WeightedGraph.jl")
include("structs/ReferenceGraph.jl")
include("structs/Graph.jl")
include("structs/GraphVertex.jl")
include("structs/GraphEdge.jl")

include("vertex_styles.jl")
include("vertex_animations.jl")
include("edge_styles.jl")
include("edge_animations.jl")

include("animations.jl")
include("overrides.jl")

# Export each function from Javis
for func in names(Javis; imported = true)
    eval(Meta.parse("import Javis." * string(func)))
    eval(Expr(:export, func))
end

export JGraph, WeightedGraph, ReferenceGraph, GraphVertex, GraphEdge

export edges, vertices, register_style_opts, add_styles, get_draw

export vertex_shape, vertex_text_style, vertex_text, vertex_fill, vertex_border

export edge_shape, edge_style, edge_arrow, edge_text

# Export each function from Luxor
for func in names(Graphs; imported = true)
    eval(Meta.parse("import Graphs." * string(func)))
    eval(Expr(:export, func))
end

export convert, @add_styles, @register_style_opts, graph_from_matrix

end # module


text:
if basename(pwd()) == "aoc"
    cd("2017/6")
end

loadblocks(filename::AbstractString) = parse.(UInt8, split(readline(filename)))

part1(filename::AbstractString) = part1(loadblocks(filename))

function part1(blocks)
    previous = copy(blocks)
    seen = Set{typeof(blocks)}()
    steps = 0
    while !(previous in seen)
        @debug join(previous, ' ')
        steps += 1
        push!(seen, previous)
        maxval, maxind = findmax(blocks)
        blocks[maxind] = zero(blocks[maxind])
        rest = maxind+1:min(maxind + maxval, lastindex(blocks))
        blocks[rest] .+= oneunit(maxval)
        maxval -= length(rest)
        if !iszero(maxval)
            d, r = divrem(maxval, oftype(maxval, length(blocks)))
            blocks .+= d
            blocks[begin:begin+r-1] .+= oneunit(maxval)
        end
        previous = copy(blocks)
    end
    @debug join(previous, ' ')
    steps, blocks
end

@assert part1("example.txt")[1] == 5

part2(blocks) = part1(part1(blocks)[2])[1]

@assert part2("example.txt")[1] == 4
part2("input.txt")[1]

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
include("cc_maximin_load_parameters_chain.jl")
include("soyster_load_parameters_chain.jl")
include("de_jong_load_parameters_chain.jl")
include("robust_load_parameters_chain.jl")


text:
include("../games/ospiel_ttt/main.jl")
# export OSpielTictactoe
Examples.games["ospiel_ttt"] = OSpielTictactoe.GameSpec()
Examples.experiments["ospiel_ttt"] = OSpielTictactoe.Training.experiment


text:
module inference

using ergm.sampler, ergm.models, ergm.stats, ergm.optim
using Statistics
import StatsBase
export mcmc_mle, mcmc_mle_from_stats

function mcmc_mle(
    observations, model :: ExponentialFamily, optimizer,
    estimation_steps, burn_in, sample_interval)
    observation_stats = hcat([get_stats(model.stats, o) for o ∈ observations]...)
    target_Es = mean(observation_stats, dims=2)[:, 1]
    # seed sampler with observations
    mcmc_mle_from_stats(target_Es, observations, model, optimizer, estimation_steps, burn_in, sample_interval)
end

function mcmc_mle_from_stats(
    target_Es, sampler_seeds, model :: ExponentialFamily, optimizer,
    estimation_steps, burn_in, sample_interval)

    m = length(get_params(model))
    θs = []
    Ls = []
    Es = []

    while !done(optimizer)
        G0 = StatsBase.sample(sampler_seeds)
        s = ParallelGibbsSampler(
            G0,
            model,
            burn_in,
            sample_interval,
            Threads.nthreads()
        )
        _, ss = sample(s, estimation_steps)
        push!(Es, mean(ss, dims=2))
        θ = optim_step(optimizer, ss, target_Es)
        update_params(model, θ)
        push!(θs, θ)
        g = mean(ss, dims=2) - target_Es
        L = sum(g .^ 2)
        push!(Ls, L)
    end

    hcat(θs...), Ls, target_Es, Es
end

function ee(
    observations, model :: ExponentialFamily,
    fitting_steps, estimation_steps, learning_rate)
end

end


text:
# SPDX-License-Identifier: X11
# 2020-08-30
# Forbidden List (300pt)

function main()
  X, N = parse.(Int, split(readline(), " "))

  if N == 0
    println(X)
    return
  end

  P = parse.(Int, split(readline(), " "))
  for i ∈ 0:cld(N, 2)
    if X - i ∉ P
      println(X - i)
      return
    elseif X + i ∉ P
      println(X + i)
      return
    end
  end
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
export InvMcmcMixtureMove

struct InvMcmcMixtureMove{N,M<:NTuple{N,InvMcmcMove}} <: InvMcmcMove
    moves::M
    p::NTuple{N,Float64}
    p_cusum::NTuple{N,Float64}
end

InvMcmcMixtureMove(moves, p) = InvMcmcMixtureMove(moves, p, cumsum(p))

acceptance_prob(x::InvMcmcMixtureMove) = [acceptance_prob(move) for move in x.moves]
function reset_counts!(x::InvMcmcMixtureMove)
    for move in x.moves
        reset_counts!(move)
    end
end

function Base.show(
    io::IO, x::InvMcmcMixtureMove
    )

    title = "InvMcmcMixtureMove ($(length(x.moves)) components)"
    println(io, title)
    for (p,move) in zip(x.p,x.moves)
        println(io, "   * $p - ", "$move")
    end
end

# We specialise the accept_reject function
function accept_reject!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    pointers::InteractionSequence{Int},
    mix_move::InvMcmcMixtureMove,
    model::T
    ) where {T<:Union{SIS,SIM}}

    # Sample move
    p = mix_move.p_cusum  # Mixture proportions
    u = rand()      # Random unif(0,1)
    i = findfirst(x->x>u, p)

    # Select ith move
    move = mix_move.moves[i]

    # Do accept-reject for the move
    accept_reject!(S_curr, S_prop, pointers, move, model)

end


text:
module ConstraintModels

using CBLS
using Constraints
using Dictionaries
using JuMP
using LocalSearchSolvers

const LS = LocalSearchSolvers

import LocalSearchSolvers: Options

export chemical_equilibrium
export golomb
export qap
export magic_square
export mincut
export n_queens
export scheduling
export sudoku

include("assignment.jl")
include("chemical_equilibrium.jl")
include("cut.jl")
include("golomb.jl")
include("magic_square.jl")
include("n_queens.jl")
include("scheduling.jl")
include("sudoku.jl")

end


text:
"""
CxxWrap wrapper module for the HelFEM shared library.
"""
module helfem
    using CxxWrap, HelFEM_jll
    @wrapmodule(libhelfem)
    function __init__()
        @initcxx
        verbose(false)
    end

    function ArmaVector(v::AbstractVector)
        av = ArmaVector(length(v))
        for (i, x) in enumerate(v)
            at!(av, i - 1, convert(Float64, x))
        end
        return av
    end
    Base.size(v::ArmaVector) = (Int(nrows(v)),)
    function Base.size(v::ArmaVector, dim::Integer)
        dim <= 0 && throw(ArgumentError("dimension $dim out of range"))
        dim == 1 ? Int(nrows(v)) : 1
    end
    Base.collect(v::ArmaVector) = [at(v, i - 1) for i = 1:size(v, 1)]

    function ArmaMatrix(M::AbstractMatrix)
        aM = ArmaMatrix(size(M, 1), size(M, 2))
        idxs = CartesianIndices(M)
        for (k, x) in enumerate(M)
            i, j = Tuple(idxs[k])
            at!(aM, i - 1, j - 1, convert(Float64, x))
        end
        return aM
    end
    Base.size(m::ArmaMatrix) = (Int(nrows(m)), Int(ncols(m)))
    function Base.size(m::ArmaMatrix, dim::Integer)
        dim <= 0 && throw(ArgumentError("dimension $dim out of range"))
        dim == 1 ? Int(nrows(m)) : dim == 2 ? Int(ncols(m)) : 1
    end
    Base.collect(m::ArmaMatrix) = [at(m, i - 1, j - 1) for i = 1:size(m, 1), j = 1:size(m, 2)]
end


text:
# !INPUT/OUTPUT PARAMETERS:
#   m : order of derivative (in,integer)
#   n : number of points (in,integer)
#   x : abscissa array (in,real(n))
#   f : function array (in,real(n))
#   g : (anti-)derivative of f (out,real(n))
# !DESCRIPTION:
#   Given function $f$ defined on a set of points $x_i$ then if $m\ge 0$ this
#   routine computes the $m$th derivative of $f$ at each point. If $m=-1$ the
#   anti-derivative of $f$ given by
#   $$ g(x_i)=\int_{x_1}^{x_i} f(x)\,dx $$
#   is calculated. Both derivatives and integrals are computed by first fitting
#   the function to a clamped cubic spline.
function fderiv!(m, n, x, f, g)

    @assert n > 0

    # automatic arrays
    cf = zeros(Float64,3,n)

    # high accuracy integration/differentiation from spline interpolation
    spline!(n, x, f, cf)
    if m == -1
        ss =0.0
        g[1] = 0.0
        for i in 1:n-1
          dx = x[i+1] - x[i]
          ss = ss + dx*( f[i] + dx*( 0.50*cf[1,i] + dx*( cf[2,i]/3 + dx*0.25*cf[3,i] ) ) )
          g[i+1] = ss
        end
    elseif m == 1
        for i in 1:n
            g[i] = cf[1,i]
        end
    elseif m == 2
        for i in 1:n
            g[i] = 2.0*cf[2,i]
        end
    elseif m == 3
        for i in 1:n
            g[i] = 6.0*cf[3,i]
        end
    else
        for i in 1:n
            g[i] = 0.0
        end
    end
    return
end


text:
using ConservationLawsDiffEq
using OrdinaryDiffEq
using LinearAlgebra
u_exact(x,t) = sin(4*pi*(x-t))
# Define the flux and flux Jacobian
Jf(u) = 1
f(u) = u

# Initial condition (using integral cell averages)
f0(x) = sin(4*pi*x)

#Setup problem for a given N (number of cells/control volumenes on a uniform mesh)
#and given final time (Tend) with periodic boundary conditions
function get_problem(N, scheme; Tend = 2.0, CFL = 0.5)

  mesh = Uniform1DFVMesh(N, [0.0, 1.0])

  f_h = getSemiDiscretization(f,scheme,mesh,[Periodic()]; Df = Jf, use_threads = false,numvars = 1)

  #Compute discrete initial data
  u0 = getInitialState(mesh,f0,use_threads = true)

  #Setup ODE problem for a time interval = [0.0,1.0]
  ode_prob = ODEProblem(f_h,u0,(0.0,Tend))

  #Setup callback in order to fix CFL constant value
  cb = getCFLCallback(f_h, CFL)

  #Estimate an initial dt
  dt = update_dt!(u0, f_h, CFL)
  return ode_prob,mesh,cb, dt
end

mesh_ncells = [40,80,160,320]
t1 = get_conv_order_table(LaxWendroffScheme(),solve, get_problem, u_exact, mesh_ncells, Euler())

using PrettyTables
header = ["M","e_tot","order"]
pretty_table(t1.data, header)


text:
#!/usr/bin/env julia --project
#=
# This experiment file establishes the initial conditions, boundary conditions,
# source terms and simulation parameters (domain size + resolution) for the
# BOMEX LES case. The set of parameters presented in the `master` branch copy
# include those that have passed offline tests at the full simulation time of
# 6 hours. Suggested offline tests included plotting horizontal-domain averages
# of key properties (see AtmosDiagnostics). The timestepper configuration is in
# `src/Driver/solver_configs.jl` while the `AtmosModel` defaults can be found in
# `src/Atmos/Model/AtmosModel.jl` and `src/Driver/driver_configs.jl`
#
# To simulate the full 6 hour experiment, change `timeend` to (3600*6) and type in
#
# julia --project experiments/AtmosLES/bomex.jl
#
# See `src/Driver/driver_configs.jl` for additional flags (e.g. VTK, diagnostics,
# update-interval, output directory settings)
#
# Upcoming changes:
# 1) Atomic sources
# 2) Improved boundary conditions
# 3) Collapsed experiment design
# 4) Updates to generally keep this in sync with master

[Siebesma2003](@cite)
=#

using ArgParse
using Distributions
using DocStringExtensions
using LinearAlgebra
using Printf
using StaticArrays
using Test

using ClimateMachine
using ClimateMachine.Atmos
using ClimateMachine.Orientations
using ClimateMachine.ConfigTypes
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.Diagnostics
using ClimateMachine.GenericCallbacks
using ClimateMachine.Mesh.Filters
using ClimateMachine.Mesh.Grids
using ClimateMachine.ODESolvers
using Thermodynamics
using ClimateMachine.TurbulenceClosures
using ClimateMachine.TurbulenceConvection
using ClimateMachine.VariableTemplates
using ClimateMachine.BalanceLaws
import ClimateMachine.BalanceLaws: prognostic_vars

using CLIMAParameters
using CLIMAParameters.Planet: e_int_v0, grav, day
using CLIMAParameters.Atmos.Microphysics

struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

import ClimateMachine.BalanceLaws: source, prognostic_vars
using ClimateMachine.Atmos: altitude, recover_thermo_state
using UnPack

"""
  Bomex Geostrophic Forcing (Source)
"""
struct BomexGeostrophic{FT} <: TendencyDef{Source}
    "Coriolis parameter [s⁻¹]"
    f_coriolis::FT
    "Eastward geostrophic velocity `[m/s]` (Base)"
    u_geostrophic::FT
    "Eastward geostrophic velocity `[m/s]` (Slope)"
    u_slope::FT
    "Northward geostrophic velocity `[m/s]`"
    v_geostrophic::FT
end

prognostic_vars(::BomexGeostrophic) = (Momentum(),)

function source(::Momentum, s::BomexGeostrophic, m, args)
    @unpack state, aux = args
    @unpack f_coriolis, u_geostrophic, u_slope, v_geostrophic = s

    z = altitude(m, aux)
    # Note z dependence of eastward geostrophic velocity
    u_geo = SVector(u_geostrophic + u_slope * z, v_geostrophic, 0)
    ẑ = vertical_unit_vector(m, aux)
    fkvector = f_coriolis * ẑ
    # Accumulate sources
    return -fkvector × (state.ρu .- state.ρ * u_geo)
end

"""
  Bomex Sponge (Source)
"""
struct BomexSponge{FT} <: TendencyDef{Source}
    "Maximum domain altitude (m)"
    z_max::FT
    "Altitude at with sponge starts (m)"
    z_sponge::FT
    "Sponge Strength 0 ⩽ α_max ⩽ 1"
    α_max::FT
    "Sponge exponent"
    γ::FT
    "Eastward geostrophic velocity `[m/s]` (Base)"
    u_geostrophic::FT
    "Eastward geostrophic velocity `[m/s]` (Slope)"
    u_slope::FT
    "Northward geostrophic velocity `[m/s]`"
    v_geostrophic::FT
end

prognostic_vars(::BomexSponge) = (Momentum(),)

function source(::Momentum, s::BomexSponge, m, args)
    @unpack state, aux = args
    @unpack z_max, z_sponge, α_max, γ = s
    @unpack u_geostrophic, u_slope, v_geostrophic = s

    z = altitude(m, aux)
    u_geo = SVector(u_geostrophic + u_slope * z, v_geostrophic, 0)
    ẑ = vertical_unit_vector(m, aux)
    # Accumulate sources
    if z_sponge <= z
        r = (z - z_sponge) / (z_max - z_sponge)
        β_sponge = α_max * sinpi(r / 2)^s.γ
        return -β_sponge * (state.ρu .- state.ρ * u_geo)
    else
        FT = eltype(state)
        return SVector{3, FT}(0, 0, 0)
    end
end

"""
    BomexTendencies (Source)

Moisture, Temperature and Subsidence tendencies
"""
struct BomexTendencies{FT} <: TendencyDef{Source}
    "Advection tendency in total moisture `[s⁻¹]`"
    ∂qt∂t_peak::FT
    "Lower extent of piecewise profile (moisture term) `[m]`"
    zl_moisture::FT
    "Upper extent of piecewise profile (moisture term) `[m]`"
    zh_moisture::FT
    "Cooling rate `[K/s]`"
    ∂θ∂t_peak::FT
    "Lower extent of piecewise profile (subsidence term) `[m]`"
    zl_sub::FT
    "Upper extent of piecewise profile (subsidence term) `[m]`"
    zh_sub::FT
    "Subsidence peak velocity"
    w_sub::FT
    "Max height in domain"
    z_max::FT
end

prognostic_vars(::BomexTendencies) = (Mass(), Energy(), TotalMoisture())

function compute_bomex_tend_params(s, m, args)
    @unpack state, aux = args
    FT = eltype(state)
    ρ = state.ρ
    z = altitude(m, aux)

    # Moisture tendencey (sink term)
    # Temperature tendency (Radiative cooling)
    # Large scale subsidence
    # Unpack struct
    @unpack zl_moisture, zh_moisture, z_max, zl_sub, zh_sub = s
    @unpack w_sub, ∂qt∂t_peak, ∂θ∂t_peak = s
    zl_temperature = zl_sub
    k̂ = vertical_unit_vector(m, aux)

    # Piecewise term for moisture tendency
    linscale_moisture = (z - zl_moisture) / (zh_moisture - zl_moisture)
    if z <= zl_moisture
        ρ∂qt∂t = ρ * ∂qt∂t_peak
    elseif zl_moisture < z <= zh_moisture
        ρ∂qt∂t = ρ * (∂qt∂t_peak - ∂qt∂t_peak * linscale_moisture)
    else
        ρ∂qt∂t = -zero(FT)
    end

    # Piecewise term for internal energy tendency
    linscale_temp = (z - zl_sub) / (z_max - zl_sub)
    if z <= zl_sub
        ρ∂θ∂t = ρ * ∂θ∂t_peak
    elseif zl_temperature < z <= z_max
        ρ∂θ∂t = ρ * (∂θ∂t_peak - ∂θ∂t_peak * linscale_temp)
    else
        ρ∂θ∂t = -zero(FT)
    end

    # Piecewise terms for subsidence
    linscale_sub = (z - zl_sub) / (zh_sub - zl_sub)
    w_s = -zero(FT)
    if z <= zl_sub
        w_s = -zero(FT) + z * (w_sub) / (zl_sub)
    elseif zl_sub < z <= zh_sub
        w_s = w_sub - (w_sub) * linscale_sub
    else
        w_s = -zero(FT)
    end
    return (w_s = w_s, ρ∂qt∂t = ρ∂qt∂t, ρ∂θ∂t = ρ∂θ∂t, k̂ = k̂)
end

function source(::Mass, s::BomexTendencies, m, args)
    @unpack state, diffusive = args
    params = compute_bomex_tend_params(s, m, args)
    @unpack ρ∂qt∂t, w_s, k̂ = params
    ρ = state.ρ
    return ρ∂qt∂t - state.ρ * w_s * dot(k̂, diffusive.moisture.∇q_tot)
end
function source(::Energy, s::BomexTendencies, m, args)
    @unpack state, diffusive = args
    @unpack ts = args.precomputed
    params = compute_bomex_tend_params(s, m, args)
    FT = eltype(state)
    @unpack ρ∂qt∂t, ρ∂θ∂t, w_s, k̂ = params
    cvm = cv_m(ts)
    param_set = parameter_set(m)
    _e_int_v0 = FT(e_int_v0(param_set))
    term1 = cvm * ρ∂θ∂t * exner(ts) + _e_int_v0 * ρ∂qt∂t
    term2 = state.ρ * w_s * dot(k̂, diffusive.energy.∇h_tot)
    return term1 - term2
end
function source(::TotalMoisture, s::BomexTendencies, m, args)
    @unpack state, diffusive = args
    params = compute_bomex_tend_params(s, m, args)
    @unpack ρ∂qt∂t, w_s, k̂ = params
    return ρ∂qt∂t - state.ρ * w_s * dot(k̂, diffusive.moisture.∇q_tot)
end

add_perturbations!(state, localgeo) = nothing

"""
  Initial Condition for BOMEX LES
"""
function init_bomex!(problem, bl, state, aux, localgeo, t)
    (x, y, z) = localgeo.coord
    # This experiment runs the BOMEX LES Configuration
    # (Shallow cumulus cloud regime)
    # x,y,z imply eastward, northward and altitude coordinates in `[m]`

    # Problem floating point precision
    FT = eltype(state)
    param_set = parameter_set(bl)

    P_sfc::FT = 1.015e5 # Surface air pressure
    qg::FT = 22.45e-3 # Total moisture at surface
    q_pt_sfc = PhasePartition(qg) # Surface moisture partitioning
    Rm_sfc = gas_constant_air(param_set, q_pt_sfc) # Moist gas constant
    θ_liq_sfc = FT(299.1) # Prescribed θ_liq at surface
    T_sfc = FT(300.4) # Surface temperature
    _grav = FT(grav(param_set))

    # Initialise speeds [u = Eastward, v = Northward, w = Vertical]
    u::FT = 0
    v::FT = 0
    w::FT = 0

    # Prescribed altitudes for piece-wise profile construction
    zl1::FT = 520
    zl2::FT = 1480
    zl3::FT = 2000
    zl4::FT = 3000

    # Assign piecewise quantities to θ_liq and q_tot
    θ_liq::FT = 0
    q_tot::FT = 0

    # Piecewise functions for potential temperature and total moisture
    if FT(0) <= z <= zl1
        # Well mixed layer
        θ_liq = 298.7
        q_tot = 17.0 + (z / zl1) * (16.3 - 17.0)
    elseif z > zl1 && z <= zl2
        # Conditionally unstable layer
        θ_liq = 298.7 + (z - zl1) * (302.4 - 298.7) / (zl2 - zl1)
        q_tot = 16.3 + (z - zl1) * (10.7 - 16.3) / (zl2 - zl1)
    elseif z > zl2 && z <= zl3
        # Absolutely stable inversion
        θ_liq = 302.4 + (z - zl2) * (308.2 - 302.4) / (zl3 - zl2)
        q_tot = 10.7 + (z - zl2) * (4.2 - 10.7) / (zl3 - zl2)
    else
        θ_liq = 308.2 + (z - zl3) * (311.85 - 308.2) / (zl4 - zl3)
        q_tot = 4.2 + (z - zl3) * (3.0 - 4.2) / (zl4 - zl3)
    end

    # Set velocity profiles - piecewise profile for u
    zlv::FT = 700
    if z <= zlv
        u = -8.75
    else
        u = -8.75 + (z - zlv) * (-4.61 + 8.75) / (zl4 - zlv)
    end

    # Convert total specific humidity to kg/kg
    q_tot /= 1000
    # Scale height based on surface parameters
    H = Rm_sfc * T_sfc / _grav
    # Pressure based on scale height
    P = P_sfc * exp(-z / H)

    # Establish thermodynamic state and moist phase partitioning
    ts = PhaseEquil_pθq(param_set, P, θ_liq, q_tot)
    T = air_temperature(ts)
    ρ = air_density(ts)

    # Compute momentum contributions
    ρu = ρ * u
    ρv = ρ * v
    ρw = ρ * w

    # Compute energy contributions
    e_kin = FT(1 // 2) * (u^2 + v^2 + w^2)
    e_pot = _grav * z
    ρe_tot = ρ * total_energy(e_kin, e_pot, ts)

    # Assign initial conditions for prognostic state variables
    state.ρ = ρ
    state.ρu = SVector(ρu, ρv, ρw)
    state.energy.ρe = ρe_tot
    state.moisture.ρq_tot = ρ * q_tot
    if moisture_model(bl) isa NonEquilMoist
        state.moisture.ρq_liq = FT(0)
        state.moisture.ρq_ice = FT(0)
    end

    add_perturbations!(state, localgeo)
    init_state_prognostic!(turbconv_model(bl), bl, state, aux, localgeo, t)
end

function bomex_model(
    ::Type{FT},
    config_type,
    zmax,
    surface_flux;
    turbconv = NoTurbConv(),
    moisture_model = "equilibrium",
) where {FT}

    ics = init_bomex!     # Initial conditions

    C_smag = FT(0.23)     # Smagorinsky coefficient

    u_star = FT(0.28)     # Friction velocity
    C_drag = FT(0.0011)   # Bulk transfer coefficient

    T_sfc = FT(300.4)     # Surface temperature `[K]`
    q_sfc = FT(22.45e-3)  # Surface specific humiity `[kg/kg]`
    LHF = FT(147.2)       # Latent heat flux `[W/m²]`
    SHF = FT(9.5)         # Sensible heat flux `[W/m²]`
    moisture_flux = LHF / latent_heat_vapor(param_set, T_sfc)

    ∂qt∂t_peak = FT(-1.2e-8)  # Moisture tendency (energy source)
    zl_moisture = FT(300)     # Low altitude limit for piecewise function (moisture source)
    zh_moisture = FT(500)     # High altitude limit for piecewise function (moisture source)
    ∂θ∂t_peak = FT(-2 / FT(day(param_set)))  # Potential temperature tendency (energy source)

    z_sponge = FT(2400)     # Start of sponge layer
    α_max = FT(0.75)        # Strength of sponge layer (timescale)
    γ = 2              # Strength of sponge layer (exponent)

    u_geostrophic = FT(-10)        # Eastward relaxation speed
    u_slope = FT(1.8e-3)     # Slope of altitude-dependent relaxation speed
    v_geostrophic = FT(0)          # Northward relaxation speed

    zl_sub = FT(1500)         # Low altitude for piecewise function (subsidence source)
    zh_sub = FT(2100)         # High altitude for piecewise function (subsidence source)
    w_sub = FT(-0.65e-2)     # Subsidence velocity peak value

    f_coriolis = FT(0.376e-4) # Coriolis parameter

    # Assemble source components
    source_default = (
        Gravity(),
        BomexTendencies{FT}(
            ∂qt∂t_peak,
            zl_moisture,
            zh_moisture,
            ∂θ∂t_peak,
            zl_sub,
            zh_sub,
            w_sub,
            zmax,
        ),
        BomexSponge{FT}(
            zmax,
            z_sponge,
            α_max,
            γ,
            u_geostrophic,
            u_slope,
            v_geostrophic,
        ),
        BomexGeostrophic{FT}(f_coriolis, u_geostrophic, u_slope, v_geostrophic),
        turbconv_sources(turbconv)...,
    )
    if moisture_model == "equilibrium"
        source = source_default
        moisture = EquilMoist(; maxiter = 5, tolerance = FT(0.1))
    elseif moisture_model == "nonequilibrium"
        source = (source_default..., CreateClouds())
        moisture = NonEquilMoist()
    else
        @warn @sprintf(
            """
%s: unrecognized moisture_model, using the defaults""",
            moisture_model,
        )
        source = source_default
        moisture = EquilMoist(; maxiter = 5, tolerance = FT(0.1))
    end

    # Set up problem initial and boundary conditions
    if surface_flux == "prescribed"
        energy_bc = PrescribedEnergyFlux((state, aux, t) -> LHF + SHF)
        moisture_bc = PrescribedMoistureFlux((state, aux, t) -> moisture_flux)
    elseif surface_flux == "bulk"
        energy_bc = BulkFormulaEnergy(
            (bl, state, aux, t, normPu_int) -> C_drag,
            (bl, state, aux, t) -> (T_sfc, q_sfc),
        )
        moisture_bc = BulkFormulaMoisture(
            (state, aux, t, normPu_int) -> C_drag,
            (state, aux, t) -> q_sfc,
        )
    else
        @warn @sprintf(
            """
%s: unrecognized surface flux; using 'prescribed'""",
            surface_flux,
        )
    end

    # Assemble model components
    physics = AtmosPhysics{FT}(
        param_set;
        turbulence = SmagorinskyLilly{FT}(C_smag),
        moisture = moisture,
        turbconv = turbconv,
    )

    problem = AtmosProblem(
        boundaryconditions = (
            AtmosBC(
                physics;
                momentum = Impenetrable(DragLaw(
                    # normPu_int is the internal horizontal speed
                    # P represents the projection onto the horizontal
                    (state, aux, t, normPu_int) -> (u_star / normPu_int)^2,
                )),
                energy = energy_bc,
                moisture = moisture_bc,
                turbconv = turbconv_bcs(turbconv)[1],
            ),
            AtmosBC(physics; turbconv = turbconv_bcs(turbconv)[2]),
        ),
        init_state_prognostic = ics,
    )

    # Assemble model components
    model =
        AtmosModel{FT}(config_type, physics; problem = problem, source = source)

    return model
end


text:
# include("admm_test.jl")
# @time test_train_logistic()
# @time test_train_logistic_optimized()
# @time test_train_logistic(size=1000*1000)
# @time test_train_logistic_optimized(size=1000*1000)

# function test_admm_without_autodiff_smallinput(;gd_iter = 3,dbg = false)
#     admm_test_matrix1 = rand([(Distributions.Bernoulli(0.7), 100 => 50, 3),(Distributions.Gaussian(3, 1), 100 => 50, 3)])
#     admm_test_matrix_missing1 = sample(BernoulliModel(), x = admm_test_matrix1, rate = 0.8)
#     @time admm_test_matrix_output_1 = complete(A = admm_test_matrix_missing1, maxiter = 200, use_autodiff = false, gd_iter = gd_iter, debug_mode = dbg)
#     gaussian_acc = accuracyImputedContinuousPart(truth = admm_test_matrix1, completedMatrix = admm_test_matrix_output_1)
#     bernoulli_acc = accuracyImputedBinaryPart(truth = admm_test_matrix1, completedMatrix = admm_test_matrix_output_1)
#     @printf("gaussian acc: %f\n", gaussian_acc)
#     @printf("bernoulli acc: %f\n",bernoulli_acc)
# end

@testset "$(format("ADMM Algorithm: Small Input[Gaussian + Bernoulli]"))" begin
  let
    truth_matrix = rand([(FixedRankMatrix(Distributions.Gaussian(0, 1), rank = 5), 200, 100),
                         (FixedRankMatrix(Distributions.Bernoulli(0.5), rank = 5), 200, 100)])
    sample_model = provide(Sampler{BernoulliModel}, rate = 0.8)
    input_matrix = sample_model.draw(truth_matrix)

  end
end

# test_admm_with_autodiff_smallinput(gd_iter=3,dbg=true)
# test_admm_without_autodiff_smallinput(gd_iter=3)
# test_admm_without_autodiff_largeinput(gd_iter=3,dbg=true)
#@time POISSON_SMALL_TEST_SET_LOOSE()


text:
"""Specification for a "time-local" pulse parametrization.

The parametrization is given as a collection of three functions:

* ``ϵ(u(t))``
* ``u(ϵ(t))``
* ``∂ϵ/∂u`` as a function of ``u(t)``.
"""
struct PulseParametrization
    epsilon_of_u :: Function
    u_of_epsilon :: Function
    de_du_derivative :: Function
end

"""Parametrization where ``ϵ(t) ≡ u(t)``."""
NoParametrization() = PulseParametrization(u->u, ϵ->ϵ, u->1.0)

"""Parametrization ϵ(t) = u²(t), enforcing pulse values ≥ 0."""
SquareParametrization() = PulseParametrization(
    u -> begin # epsilon_of_u
        ϵ = u^2
    end,
    ϵ -> begin # u_of_epsilon
        a = max(ϵ, 0.0)
        u = √a
    end,
    u -> begin # de_du_derivative
        ∂ϵ╱∂u = 2u
    end
)

"""Parametrization with a tanh function that enforces ϵ_min < ϵ(t) < ϵ_max.
"""
function TanhParametrization(ϵ_min, ϵ_max)

    Δ = ϵ_max - ϵ_min
    Σ = ϵ_max + ϵ_min
    εₚ = eps(1.0)  # 2⋅10⁻¹⁶ (machine precision)
    @assert ϵ_max > ϵ_min

    PulseParametrization(
        u -> begin # epsilon_of_u
            ϵ = tanh(u) * Δ / 2 + Σ / 2
        end,
        ϵ -> begin # u_of_epsilon
            a = clamp(2ϵ/Δ -  Σ/Δ, -1+εₚ, 1-εₚ)
            u = atanh(a)  # -18.4 < u < 18.4
        end,
        u -> begin # de_du_derivative
            ∂ϵ╱∂u = (Δ / 2) * sech(u)^2
        end
    )

end

"""Parametrization with a tanh² function that enforces ``0 ≤ ϵ(t) < ϵ_max``.
"""
function TanhSqParametrization(ϵ_max)

    εₚ = eps(1.0)  # 2⋅10⁻¹⁶ (machine precision)
    @assert ϵ_max > 0

    PulseParametrization(
        u -> begin # epsilon_of_u
            ϵ = ϵ_max * tanh(u)^2
        end,
        ϵ -> begin # u_of_epsilon
            a = clamp(ϵ / ϵ_max, 0, 1-εₚ)
            u = atanh(√a)
        end,
        u -> begin # de_du_derivative
            ∂ϵ╱∂u = 2ϵ_max * tanh(u) * sech(u)^2
        end
    )

end

"""Parametrization with a Logistic function that enforces ϵ_min < ϵ(t) < ϵ_max.
"""
function LogisticParametrization(ϵ_min, ϵ_max; k=1.0)

    Δ = ϵ_max - ϵ_min
    ε₀ = eps(0.0)  # 5⋅10⁻³²⁴
    @assert ϵ_max > ϵ_min

    PulseParametrization(
        u -> begin # epsilon_of_u
            ϵ = Δ / (1 + exp(-k * u)) + ϵ_min
        end,
        ϵ -> begin # u_of_epsilon
            ϵ′ = ϵ - ϵ_min
            a = max(ϵ′/(Δ - ϵ′), ε₀)
            u = log(a) / k
        end,
        u -> begin # de_du_derivative
            e⁻ᵏᵘ = exp(-k * u)
            ∂ϵ╱∂u = Δ * k * e⁻ᵏᵘ / (1 + e⁻ᵏᵘ)^2
        end
    )

end

"""Parametrization with a Logistic-Square function that enforces 0 ≤ ϵ(t) < ϵ_max.
"""
function LogisticSqParametrization(ϵ_max; k=1.0)

    ε₀ = eps(0.0)  # 5⋅10⁻³²⁴
    @assert ϵ_max > 0

    PulseParametrization(
        u -> begin # epsilon_of_u
            ϵ = ϵ_max * (2 / (1 + exp(-k * u)) - 1)^2
        end,
        ϵ -> begin # u_of_epsilon
            ρ = clamp(ϵ / ϵ_max, 0.0, 1.0)
            a = clamp((2 / (√ρ + 1)) - 1, ε₀, 1.0)
            u = -log(a) / k
        end,
        u -> begin # de_du_derivative
            eᵏᵘ = exp(k * u)
            ∂ϵ╱∂u = 4k * ϵ_max * eᵏᵘ * (eᵏᵘ  - 1) / (eᵏᵘ + 1)^3
        end
    )

end


text:
export ModuleProvider, dispose

# TODO: use this class implicitly? ie. when passing Modules to functions expecting a
#       ModuleProvider. Especially because dispose(mp) destroys the underlying mod.

@checked struct ModuleProvider
    ref::API.LLVMModuleProviderRef
end
reftype(::Type{ModuleProvider}) = API.LLVMModuleProviderRef

ModuleProvider(mod::Module) =
    ModuleProvider(API.LLVMCreateModuleProviderForExistingModule(ref(mod)))

function ModuleProvider(f::Core.Function, args...)
    mp = ModuleProvider(args...)
    try
        f(mp)
    finally
        dispose(mp)
    end
end

# NOTE: this destroys the underlying module
dispose(mp::ModuleProvider) = API.LLVMDisposeModuleProvider(ref(mp))


text:
using SparseArrays

struct InitialParams{L}
    centers::Matrix{Float64}
    covs::Array{CovMat{L},1};
    assignment::Vector{Int64};

    InitialParams(centers::Matrix{<:Real}, covs::Float64, assignment::Vector{Int64}) where N =
        new{N}(deepcopy(centers), [CovMat{N}(diagm(0 => (ones(N) .* covs))) for i in 1:size(centers, 1)], deepcopy(assignment))

    InitialParams(centers::Array{T, 2} where T <: Real, covs::Array{<:CovMat{N}, 1}, assignment::Array{Int64,1}) where N =
        new{N}(deepcopy(centers), deepcopy(covs), deepcopy(assignment))
end


text:
using StochasticDiffEq, Test, Random, DiffEqDevTools
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
using DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_2Dlinear
Random.seed!(100)
prob = prob_sde_2Dlinear

## Solve and plot
println("Solve and Plot")
sol = solve(prob,EM(),dt=1/2^(3))
sol = solve(prob,LambaEM(),dt=1/2^(3))
sol = solve(prob,LambaEulerHeun(),dt=1/2^(3))
sol = solve(prob,RKMil(),dt=1/2^(3))
sol = solve(prob,RKMilCommute(),dt=1/2^(3))
sol = solve(prob,RKMilGeneral(),dt=1/2^(3))
sol = solve(prob,SRI(),dt=1/2^(3))
sol = solve(prob,SRIW1(),dt=1/2^(3))
sol = solve(prob,SRIW2(),dt=1/2^(3))
sol = solve(prob,SOSRI(),dt=1/2^(3))
sol = solve(prob,SOSRI2(),dt=1/2^(3))
sol = solve(prob,ImplicitEM(),dt=1/2^(3))
sol = solve(prob,ImplicitEM(nlsolve=StochasticDiffEq.NLFunctional()),dt=1/2^(3),adaptive=false)
sol = solve(prob,ImplicitEM(autodiff=false),dt=1/2^(3))
sol = solve(prob,ImplicitRKMil(),dt=1/2^(3))

sol = solve(prob,SRIW1(),dt=1/2^(3),save_everystep=false)
#sol = solve(prob,SRIW1(),dt=1/2^(3),progress=true,progress_steps=1)

sol = solve(prob,SRIW1(),dt=1/2^(3),seed=1)
sol2 = solve(prob,SRI(error_terms=2),dt=1/2^(3),seed=1,delta=1/6)
sol.t ≈ sol2.t

sol = solve(prob,SOSRI(),dt=1/2^(3),seed=1)
sol2 = solve(prob,SRI(tableau=StochasticDiffEq.constructSRIOpt1()),dt=1/2^(3),seed=1)
sol.t ≈ sol2.t

## Convergence Testing
println("Convergence Test on 2D Linear")
dts = (1/2) .^ (7:-1:4) #14->7 good plot

sim = test_convergence(dts,prob,EM(),trajectories=1000)
@test abs(sim.𝒪est[:l2]-.5) < 0.1

sim = test_convergence(dts,prob,LambaEM(),trajectories=1000)
@test abs(sim.𝒪est[:l2]-.5) < 0.1

sim = test_convergence(dts,prob,ImplicitEM(),trajectories=300)
@test abs(sim.𝒪est[:l2]-.5) < 0.1

sim = test_convergence(dts,prob,ImplicitEM(theta=1),trajectories=100)
@test abs(sim.𝒪est[:l2]-.5) < 0.1

sim = test_convergence(dts,prob,ImplicitEM(symplectic=true),trajectories=500)
@test abs(sim.𝒪est[:l2]-.5) < 0.1

sim = test_convergence(dts,prob,ImplicitEM(symplectic=true,autodiff=false),trajectories=100)
@test abs(sim.𝒪est[:l2]-.5) < 0.1

sim = test_convergence(dts,prob,ISSEM(),trajectories=500)
@test abs(sim.𝒪est[:l2]-.5) < 0.1

sim = test_convergence(dts,prob,ImplicitRKMil(),trajectories=100)
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim = test_convergence(dts,prob,ImplicitRKMil(theta=1),trajectories=100)
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim = test_convergence(dts,prob,ImplicitRKMil(theta=1,autodiff=false),trajectories=200)
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim = test_convergence(dts,prob,ImplicitRKMil(symplectic=true),trajectories=150)
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim = test_convergence(dts,prob,ImplicitRKMil(symplectic=true,autodiff=false),trajectories=100)
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim = test_convergence(dts,prob,ImplicitRKMil(nlsolve=StochasticDiffEq.NLFunctional()),trajectories=100)
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim2 = test_convergence(dts,prob,RKMil(),trajectories=1000)
@test abs(sim2.𝒪est[:l∞]-1) < 0.2

sim2 = test_convergence(dts,prob,RKMilCommute(),trajectories=1000)
@test abs(sim2.𝒪est[:l∞]-1) < 0.2

sim2 = test_convergence(dts,prob,RKMilGeneral(),trajectories=1000)
@test abs(sim2.𝒪est[:l∞]-1) < 0.2

print(".")

sim2 = test_convergence(dts,prob,WangLi3SMil_A(),trajectories=100)
@test abs(sim2.𝒪est[:l∞]-1) < 0.2

sim2 = test_convergence(dts,prob,WangLi3SMil_B(),trajectories=100)
@test abs(sim2.𝒪est[:l∞]-1) < 0.2

sim2 = test_convergence(dts,prob,WangLi3SMil_C(),trajectories=100)
@test abs(sim2.𝒪est[:l∞]-1) < 0.2

sim2 = test_convergence(dts,prob,WangLi3SMil_D(),trajectories=100)
@test abs(sim2.𝒪est[:l∞]-1) < 0.2

sim2 = test_convergence(dts,prob,WangLi3SMil_E(),trajectories=100)
@test abs(sim2.𝒪est[:l∞]-1) < 0.2

sim2 = test_convergence(dts,prob,WangLi3SMil_F(),trajectories=100)
@test abs(sim2.𝒪est[:l∞]-1) < 0.2

print(".")

Random.seed!(100)
eigen_est = (integrator) -> integrator.eigen_est = 10.0
for Alg in [SROCK1, SROCK2], alg in [Alg(), Alg(eigen_est=eigen_est)]
  local sim2
  sim2 = test_convergence(dts,prob,alg,trajectories=150)
  @test abs(sim2.𝒪est[:l∞]-1) < 0.2
  @show sim2.𝒪est[:l∞]
end

sim2 = test_convergence(dts,prob,SROCKEM(strong_order_1=false),trajectories=100)
@test abs(sim2.𝒪est[:l∞]-0.5) < 0.2

sim2 = test_convergence(dts,prob,SROCKEM(),trajectories=100)
@test abs(sim2.𝒪est[:l∞]-1) < 0.2

sim2 = test_convergence(dts,prob,SROCKEM(eigen_est=eigen_est),trajectories=100)
@test abs(sim2.𝒪est[:l∞]-1) < 0.2

sim2 = test_convergence(dts,prob,SKSROCK(),trajectories=500)
@test abs(sim2.𝒪est[:l∞]-0.5) < 0.2
sim2 = test_convergence(dts,prob,SKSROCK(eigen_est=eigen_est),trajectories=500)
@test abs(sim2.𝒪est[:l∞]-0.5) < 0.2

sim2 = test_convergence(dts,prob,SROCKC2(),trajectories=100)
@test abs(sim2.𝒪est[:l∞]-1) < 0.2
sim2 = test_convergence(dts,prob,SROCKC2(eigen_est=eigen_est),trajectories=100)
@test abs(sim2.𝒪est[:l∞]-1) < 0.2

# #omitting tests for incomplete methods
# sim = test_convergence(dts,prob,TangXiaoSROCK2(version_num=1),trajectories=100)
# @test abs(sim.𝒪est[:l∞]- 1) < 0.2
#
# sim = test_convergence(dts,prob,TangXiaoSROCK2(version_num=2),trajectories=100)
# @test abs(sim.𝒪est[:l∞]- 1) < 0.2
#
# sim = test_convergence(dts,prob,TangXiaoSROCK2(version_num=3),trajectories=100)
# @test abs(sim.𝒪est[:l∞]- 1) < 0.2
#
# sim = test_convergence(dts,prob,TangXiaoSROCK2(version_num=4),trajectories=100)
# @test abs(sim.𝒪est[:l∞]- 1) < 0.2
#
# sim = test_convergence(dts,prob,TangXiaoSROCK2(version_num=5),trajectories=100)
# @test abs(sim.𝒪est[:l∞]- 1) < 0.2

print(".")

sim3 = test_convergence(dts,prob,SRI(),trajectories=100)
@test abs(sim3.𝒪est[:final]-1.5) < 0.3

sim4 = test_convergence(dts,prob,SRIW1(),trajectories=100)
@test abs(sim4.𝒪est[:final]-1.5) < 0.3

sim4 = test_convergence(dts,prob,SRIW2(),trajectories=100)
@test abs(sim4.𝒪est[:final]-1.5) < 0.3

sim4 = test_convergence(dts,prob,SOSRI(),trajectories=100)
@test abs(sim4.𝒪est[:final]-1.5) < 0.3

sim4 = test_convergence(dts,prob,SOSRI2(),trajectories=100)
@test abs(sim4.𝒪est[:final]-1.5) < 0.3

# 2D oop
f_oop(u,p,t) = u
g_oop(u,p,t) = u
prob = SDEProblem(f_oop, g_oop, ones(2, 2), (0., 1.))
@test_nowarn solve(prob, ImplicitEM())

max_stars_repo_name: JuliaDiffEq/StochasticDiffEq


text:
ENV["PYCALL_JL_RUNTIME_PYTHON"] = Sys.which("python")
using Distributed
# using 6 cores
if 6 - Distributed.nprocs() > 0
    Distributed.addprocs(6 - Distributed.nprocs());
end
@everywhere using SACBP
using SACBP
using LinearAlgebra
using Plots
pyplot()
using ProgressMeter
using Printf
using Random
using JLD2
using FileIO

include(joinpath(@__DIR__, "../baseline_implementations/p_control.jl"))
include(joinpath(@__DIR__, "../baseline_implementations/mcts.jl"))
include(joinpath(@__DIR__, "../baseline_implementations/belief_ilqg.jl"))
include(joinpath(@__DIR__, "../baseline_implementations/tlqg.jl"))

function get_initial_true_state()
    px0_true = 1.5;    # True px(0)
    py0_true = 1.;     # True py(0)
    θ0_true = 11*pi/6.;# True theta(0) in (0,2*pi)
    vx0_true = 0.;     # True vx(0)
    vy0_true = 0.;     # True vy(0)
    ω0_true = 0.;      # True ω(0)
    m_true = 2.0;      # True mass
    a_true = 0.5;      # Height of the Plate.
    b_true = 1.0;      # Width of the Plate.
    J_true = 1/3*m_true*(a_true^2 + b_true^2); # Moment of Inertia
    rx_true = 0.25;    # Arm Length rx
    ry_true = 0.25;    # Arm Length ry
    μ_true = 4.0;      # Linear Friction Coefficient
    x_init = PhysManipulate2D(0.,[px0_true,py0_true],θ0_true,[vx0_true,vy0_true],ω0_true,m_true,J_true,[rx_true,ry_true],μ_true);
    return x_init, a_true, b_true
end

function get_initial_belief_state()
    μ_init = [4.,4.,pi/4,0.1,-0.1,pi/10,6.,2.,0.3,0.1,7.0];
    Σ_init = diagm([10.,10.,(pi/2)^2,2.,2.,(pi/4)^2,1.,1.0,1.0,1.0,5.0]);
    b_init = BelMvNormal(0.,μ_init,Σ_init);
    return b_init
end

function visualizeSimEnv(px::Real,py::Real,θ::Real,vx::Real,vy::Real,ω::Real,
                         a::Real,b::Real,rx::Real,ry::Real,dto=0.0;
                         plot_velocity=false)
    plt = plot();
    Vertices = [[-b/2,a/2],[-b/2,-a/2],[b/2,-a/2],[b/2,a/2]];
    xVertices = map(v -> v[1], Vertices);
    yVertices = map(v -> v[2], Vertices);
    VRotated = map(v -> [cos(θ) -sin(θ); sin(θ) cos(θ)]*v + [px, py], Vertices);
    xRotated = map(v -> v[1], VRotated);
    yRotated = map(v -> v[2], VRotated);
    Robot = [cos(θ) -sin(θ); sin(θ) cos(θ)]*[rx, ry] + [px, py];
    plot!(xRotated,yRotated,seriestype=:shape,color=:cyan,linecolor=:darkcyan,alpha=0.4,aspect_ratio=1.0,label="Object",xlabel="x [m]",ylabel="y [m]");
    scatter!((px,py),color=:darkcyan,markersize=5.0,label="C.M.")
    scatter!((Robot[1],Robot[2]),color=:red,markersize=10.0,label="Robot")
    if plot_velocity
        plot!([px,px+dto*vx],[py,py+dto*vy],line=:arrow,color=:green,width=2.0,label="Lin. Vel.")
        plot!([px + min(a,b)/4*cos(pi/2+ii) for ii = range(0.,stop=5*dto*ω,length=10)],
              [py + min(a,b)/4*sin(pi/2+ii) for ii = range(0.,stop=5*dto*ω,length=10)],line=:arrow,color=:orange,width=2.0,label="Ang. Vel.");
    end
    plot!(xVertices,yVertices,seriestype=:shape,color=:black,linecolor=:black,alpha=0.1,label="Goal",size=(500,400))
    xlims!((-1.,2.3))
    ylims!((-1.,1.5))
    return plt
end

function simulate(b_init::BelMvNormal{<:Real},
                  x_init::PhysManipulate2D{<:Real},
                  a_true::Real, b_true::Real,
                  controlMode::String, rng_num::Int64,
                  simModel::SimulateManipulate2D;
                  numSamples=10, # for sacbp only
                  u_param_min::Vector{Float64}=[-3.0, -3.0, -3.0],
                  u_param_max::Vector{Float64}=[3.0, 3.0, 3.0],
                  posGain::Real,
                  rotGain::Real,
                  Th::Float64=2.0, # planning horizon [s]
                  tf::Float64=20.0, # simulation horizon [s]
                  mcts_param::Union{Nothing, MCTSParams}=nothing,
                  ilqgPolicy::Union{Nothing, Dict{String,Any}}=nothing,
                  tlqg::Union{Nothing, TLQG_Planner_Manipulate2D}=nothing,
                  animate=true, plotfig=true, verbose=true)
    rng = MersenneTwister(rng_num)
    @assert controlMode in ["sacbp", "mcts", "ilqg", "tlqg", "tlqg_offline"];
    if animate
        # Animation Object.
        anim = Animation();
    end
    # Only for MCTS policies.
    if controlMode == "mcts"
        @assert !isnothing(mcts_param) "MCTS Parameters not given"
        bmdp = BelMDP(simModel, u_param_min, u_param_max, mcts_param.numActions)
        mctsSolverPControl = MCTS.DPWSolver(depth=Int64(Th/simModel.dto),
                                            n_iterations=mcts_param.numSamples,
                                            exploration_constant=mcts_param.expConst,
                                            k_action=mcts_param.kAct,
                                            alpha_action=mcts_param.αAct,
                                            k_state=mcts_param.kState,
                                            alpha_state=mcts_param.αState,
                                            rng=rng,
                                            estimate_value=RolloutEstimator(PControl(u_param_min, u_param_max, posGain, rotGain)));
        mctsPolicyPControl = solve(mctsSolverPControl,bmdp);
    end

    contTransModel = simModel.continuousTransModel
    discTransModel = simModel.discreteTransModel
    dtc, dto = simModel.dtc, simModel.dto
    Q, R, Cs, Cu = simModel.Q, simModel.R, simModel.Cs, simModel.Cu

    # Controller Compilation.
    U_default = [MControl2D(round(dtc*t,digits=5),zeros(3)) for t = 0:Int64(Th/dtc)-1];
    if controlMode == "sacbp" # SACBP (with PControl as nominal)
        @time U_default = pControlUpdate(simModel,b_init,U_default,u_param_min,u_param_max,posGain,rotGain);
        @time U_default = sacControlUpdate(simModel,b_init,U_default,u_param_min,u_param_max,numSamples,rng)[1];
    elseif controlMode == "mcts" # MCTS (with PControl as nominal)
        @time U_default = mctsControlUpdate(mctsPolicyPControl,b_init,U_default,dto)
    elseif controlMode == "ilqg" # Belief iLQG (initialized with PControl)
        @assert !isnothing(ilqgPolicy) "iLQG Policy file not given"
        @time U_default = ilqgControlUpdate(ilqgPolicy,b_init,U_default,u_param_min,u_param_max,dto,dto)
    elseif controlMode == "tlqg"
        @assert !isnothing(tlqg) "T-LQG Planner object not given"
        @time ~,U_default = tlqgControlUpdate(tlqg,b_init,U_default,verbose=verbose,online=true)
    elseif controlMode == "tlqg_offline"
        # offline == planner can take much longer time to compute NLP iterations (until convergence).
        # Only for collecting computation time statistics
        @assert !isnothing(tlqg) "T-LQG Planner object not given"
        @time ~,U_default = tlqgControlUpdate(tlqg,b_init,U_default,verbose=verbose,online=false)
    end

    println("Controller Compiled. Starting Simulation...")

    # Simulation.
    y_history = Vector{Float64}[];
    tcalc_true_history = Float64[];
    x_history = [x_init];
    b_history = [b_init];
    u_history = MControl2D{Float64}[];
    U_pool = copy(U_default);

    observeModel = ObserveModel_Manipulate2D();
    numSteps = Int64(tf/dtc);
    nextSampleTime = dto;
    @showprogress 1 for t = 1:numSteps
        if round(t*dtc,digits=5) > nextSampleTime
        # Belief & Control Update.
            x = x_history[end];
            if animate
                plt = visualizeSimEnv(x.pos[1],x.pos[2],x.θ,x.vel[1],x.vel[2],x.ω,a_true,b_true,x.r[1],x.r[2],dto,plot_velocity=true)
                time_str = @sprintf "%3.2f" x.t
                title!("Time: $(time_str) [s]")
                Plots.frame(anim);
            end
            u = u_history[end];
            # Observation.
            yVec = observe(observeModel,x,MControl2D(x.t,vec(u)),R,rng);
            push!(y_history,yVec);
            # Belief Update.
            b = pop!(b_history);
            b_new = trans(discTransModel,b,MControl2D(b.t,vec(u)),yVec,R);
            push!(b_history,b_new);
            if verbose
                println((b_new.t,"Belief Updated."));
            end
            # Control Update
            if controlMode == "sacbp";  # SACBP (with PControl)
                tcalc_true = @elapsed U_pool = pControlUpdate(simModel,b_new,U_pool,u_param_min,u_param_max,posGain,rotGain);
                tcalc_true += @elapsed U_pool,act_time_init,act_time_final,tcalc_simulated = sacControlUpdate(simModel,b_new,U_pool,u_param_min,u_param_max,numSamples,rng,offline=true);
                push!(tcalc_true_history,tcalc_true);
                if verbose
                    println((b.t,"Control Updated.","time_init=$(act_time_init)","time_final=$(act_time_final)"));
                end
            elseif controlMode == "mcts" # MCTS (with PControl)
                tcalc_true = @elapsed U_pool = mctsControlUpdate(mctsPolicyPControl,b_new,U_pool,dto);
                push!(tcalc_true_history,tcalc_true);
                if verbose
                    println((b.t,"Control Updated."));
                end
            elseif controlMode == "ilqg" # Belief iLQG (with PControl)
                tcalc_true = @elapsed U_pool = ilqgControlUpdate(ilqgPolicy,b_new,U_pool,u_param_min,u_param_max,dto,dto);
                push!(tcalc_true_history,tcalc_true);
                if verbose
                    println((b.t,"Control Updated."));
                end
            elseif controlMode == "tlqg" # T-LQG (with PControl)
                #println(b_new)
                tcalc_true, U_pool = tlqgControlUpdate(tlqg,b_new,U_pool,verbose=verbose,online=true)
                if !isnothing(tcalc_true)
                    push!(tcalc_true_history,tcalc_true);
                end
                if verbose
                    println((b.t,"Control Updated."));
                end
            elseif controlMode == "tlqg_offline" # T-LQG (with PControl), offline
                tcalc_true, U_pool = tlqgControlUpdate(tlqg,b_new,U_pool,verbose=verbose,online=false)
                if !isnothing(tcalc_true)
                    push!(tcalc_true_history,tcalc_true);
                end
                if verbose
                    println((b.t,"Control Updated."));
                end
            end
            # Control Execution.
            u = popfirst!(U_pool);
            push!(u_history,u);
            x_new = trans(TransModel_Manipulate2D(),x,u,dtc);
            push!(x_history,x_new);
            b_new = trans(contTransModel,b_new,MControl2D(b.t,vec(u)),dtc,Q);
            push!(b_history,b_new);

            # Append nominal control.
            u_nominal = MControl2D(round(b.t+Th,digits=5),zeros(3));
            push!(U_pool,u_nominal);
            nextSampleTime = round(nextSampleTime+dto,digits=5);
        else
            # No Control Update.
            b = b_history[end];
            # Control Execution.
            x = x_history[end];
            if animate && round(x.t/0.02, digits=5) == round(x.t/0.02)
                plt = visualizeSimEnv(x.pos[1],x.pos[2],x.θ,x.vel[1],x.vel[2],x.ω,a_true,b_true,x.r[1],x.r[2],dto,plot_velocity=true)
                time_str = @sprintf "%3.2f" x.t
                title!("Time: $(time_str) [s]")
                Plots.frame(anim);
            end
            u = popfirst!(U_pool);
            push!(u_history,u);
            x_new = trans(TransModel_Manipulate2D(),x,u,dtc);
            push!(x_history,x_new);
            b_new = trans(contTransModel,b,u,dtc,Q);
            push!(b_history,b_new);

            # Append nominal control.
            u_nominal = MControl2D(round(b.t+Th,digits=5),zeros(3));
            push!(U_pool,u_nominal);
        end
    end
    # Observation.
    x = x_history[end];
    if animate
        plt = visualizeSimEnv(x.pos[1],x.pos[2],x.θ,x.vel[1],x.vel[2],x.ω,a_true,b_true,x.r[1],x.r[2],dto,plot_velocity=true)
        time_str = @sprintf "%3.2f" x.t
        title!("Time: $(time_str) [s]")
        Plots.frame(anim);
    end
    u = u_history[end];
    # Observation.
    yVec = observe(observeModel,x,MControl2D(x.t,vec(u)),R,rng);
    push!(y_history,yVec);
    # Belief Update.
    b = pop!(b_history);
    b_new = trans(discTransModel,b,MControl2D(x.t,vec(u)),yVec,R);
    push!(b_history,b_new);
    if verbose
        println((b_new.t,"Belief Updated."));
    end
    log_dir = normpath(joinpath(@__DIR__, "../data",
                                "manipulation_with_unknown_parameters"));
    if animate
        gif(anim,joinpath(log_dir, "$(controlMode)_tf_$(tf)_rng_$(rng_num).gif"),fps=50)
        #gif(anim,joinpath(log_dir, "$(controlMode)_tf_$(tf)_rng_$(rng_num).gif"),fps=5)
    end
    return y_history,tcalc_true_history,b_history,x_history,u_history,U_pool
end

function simulate_main(controlMode::String; rng_num::Int64,
                       animate=true, plotfig=true, verbose=true)
    # PControl Parameters
    posGain = 1.0;
    rotGain = 0.5;

    # SAC Controller Parameters.
    dtc = 0.01;        # Euler Approximation Interval.
    dto = 0.2;         # Observation & Control Update Interval.
    Th = 2.0;          # Planning Horizon.
    dtexec = 0.2*dto;  # Control Execution Duration.
    Q = diagm([0.05, 0.05, 1.0*pi/180, 0.05, 0.05, 1.0*pi/180., 0., 0., 0., 0., 0.]); # state: [px, py, θ, vx, vy, ω, m, J, rx, ry, μ]
    R = diagm([0.1, 0.1, 5.0*pi/180, 0.1, 0.1, 5.0*pi/180., 0.2, 0.2, 10.0*pi/180.]) # observation: [px, py, θ, vx, vy, ω, ax, ay(linear acceleration), α(angular acceleration)]
    Cs = diagm([20.0, 20.0, 20.0, 15.0, 15.0, 15.0, 0., 0., 0., 0., 0.]); # State Cost Coefficient.
    Cu = diagm([0.1, 0.1, 0.1]);  # Control Cost Coefficient.
    u_param_min = [-3.,-3.,-3.];  # Lower Limit for PosControl
    u_param_max = [3.,3.,3.];     # Upper Limit for PosControl
    contTransModel = CPredictManipulate2D();
    discTransModel = DUpdateManipulate2D();
    costModel = CostManipulate2D();
    simModel = SimulateManipulate2D(contTransModel,discTransModel,costModel,dtc,dto,dtexec,Q,R,Cs,Cu);
    numSamples = 10;     # Number of Monte Carlo Samples for SACBP.

    # MCTS Parameters
    numActions = 7;             # Number of Actions per Coordinate.
    numSamples_mcts = 240;
    expConst = 100.0            # Exploration Constant.
    kAct = 8.0;                 # Action k.
    αAct = 0.2;                 # Action α.
    kState = 8.0;               # State k.
    αState = 0.2;               # State α.
    mcts_param = MCTSParams(numActions, numSamples_mcts, expConst, kAct, αAct, kState, αState);

    # iLQG Policy
    ilqg_policy_path = normpath(joinpath(@__DIR__, "../baseline_implementations",
                                         "ilqg_pcontrol_policy_20sec_2.jld2"));
    ilqgPolicy = FileIO.load(ilqg_policy_path);

    # Simulation Parameters
    tf = 20.;    # Simulation Time.

    # T-LQG Planner
    Th_tlqg = Th;
    kl_threshold = 25.0;
    tlqg = TLQG_Planner_Manipulate2D(Th_tlqg,dto,Q,R,Cs,Cu,u_param_min,u_param_max,posGain,rotGain,kl_threshold)

    x_init, a_true, b_true = get_initial_true_state();
    if plotfig
        frame = visualizeSimEnv(x_init.pos[1],x_init.pos[2],x_init.θ,
                                x_init.vel[1],x_init.vel[2],x_init.ω,
                                a_true, b_true, x_init.r[1], x_init.r[2])
        init_fig_path = normpath(joinpath(@__DIR__, "../data",
                                          "manipulation_with_unknown_parameters",
                                          "manipulation_initial.pdf"))
        savefig(init_fig_path);
    end

    b_init = get_initial_belief_state();

    y_history,tcalc_true_history,b_history,x_history,u_history,U_pool =
        simulate(b_init, x_init, a_true, b_true, controlMode, rng_num,
                 simModel, numSamples=numSamples,
                 u_param_min=u_param_min, u_param_max=u_param_max,
                 posGain=posGain, rotGain=rotGain,
                 Th=Th, tf=tf, mcts_param=mcts_param,
                 ilqgPolicy=ilqgPolicy, tlqg=tlqg,
                 animate=animate, plotfig=plotfig, verbose=verbose);
    println("Average control computation time: $(round(mean(tcalc_true_history), digits=3)) [s]")
    return y_history,tcalc_true_history,b_history,x_history,u_history,U_pool, Cs, Cu
end

if abspath(PROGRAM_FILE) == @__FILE__
    println("1. Simulation with SACBP Controller")
    simulate_main("sacbp", rng_num=202, animate=true, plotfig=true, verbose=false);
    println("2. Simulation with MCTS-DPW Controller")
    simulate_main("mcts", rng_num=202, animate=true, plotfig=true, verbose=false);
    println("3. Simulation with Belief iLQG Controller")
    simulate_main("ilqg", rng_num=202, animate=true, plotfig=true, verbose=false);
    println("4. Simulation with T-LQG Controller")
    simulate_main("tlqg", rng_num=202, animate=true, plotfig=true, verbose=false);
end


text:
include("./utils.jl")
include("../src/unionfind.jl")
using .UnionFindAlgo

PROJECT_PATH = joinpath(@__DIR__, "..")
tinyuf = loadufdata(joinpath(PROJECT_PATH, "data/tinyUF.txt"))
mediumuf = loadufdata(joinpath(PROJECT_PATH, "data/mediumUF.txt"))

@testset "Union find" begin
    uf = UnionFind(tinyuf)
    @test uf.ncomp == 2
    @test isconnected(uf, 3, 4) == false
    @test isconnected(uf, 4, 5) == true
    @test isconnected(uf, 9, 10) == true

    uf = UnionFind(mediumuf)
    @test uf.ncomp == 3
end

max_stars_repo_name: Jetafull/AlgorithmsInJulia

max_issues_repo_name: Jetafull/algorithms-in-julia

max_forks_repo_name: Jetafull/algorithms-in-julia


text:
function save_radau_table_to_file(n_rule::Int64)
    @eval begin  # https://stackoverflow.com/questions/26071317/declaring-top-level-variables-in-julia-using-metaprogramming/26071597#26071597
        n_rule = $n_rule
        A, b, c, bi, b_hat = radau_butcher_table_plus(2 * n_rule - 1)
        inv_A = inv(A)
        lambda, T = big_eigen(inv_A)
        inv_T = inv(T)
        f_64_type(v::VecOrMat{BigFloat}) = Float64
        f_64_type(v::VecOrMat{Complex{BigFloat}}) = Complex{Float64}
        symbols_save = (:A, :inv_A, :b, :c, :lambda, :T, :inv_T, :bi, :b_hat)
        symbols_type = f_64_type.(eval.(symbols_save))

        ### Table Data
        radau_table_path = get_radau_table_path()
        save_path = joinpath(radau_table_path, string(n_rule) * "_rule")
        for (k, var) = enumerate(symbols_save)
            file_path = joinpath(save_path, string(var) * ".txt")
            open(file_path, "w") do io
                d_type_k = symbols_type[k]
                writedlm(io, d_type_k.(eval(var)))
            end
        end
    end

    ### Table Types
    save_path_types = joinpath(radau_table_path, "type.txt")
    open(save_path_types, "w") do io
        writedlm(io, string.(symbols_type))
    end

    ### Table Symbols
    save_path_names = joinpath(radau_table_path, "symbols.txt")
    open(save_path_names, "w") do io
        writedlm(io, string.(symbols_save))
    end
end


text:
module Resparse

include("SvdLowrank.jl")

end

max_stars_repo_name: SamanthaWangdl/ReSparse

max_issues_repo_name: SamanthaWangdl/ReSparse

max_forks_repo_name: SamanthaWangdl/ReSparse