text:
export velocities, velocities!
export velocities_gen, velocities_gen!

function velocities(
    x::ArrayPartition{F, T},
    p::Union{SampledModel, IntegratedModel},
    t
) where {
    F,
    T <: Tuple{Vararg{AbstractVector{<:Real}}}
}
    dx = similar(x)
    velocities!(dx, x, p, t)
    dx
end

function velocities!(
    dx::ArrayPartition{F, T},
    x::ArrayPartition{F, T},
    p::SampledModel{N, TVs, TWprimes, Tmobilities},
    t
) where {
    F,
    T <: Tuple{Vararg{AbstractVector{<:Real}}},
    N, TVs, TWprimes, Tmobilities
}
    dens = pwc_densities(x.x)
    for spec in 1:N
        dx.x[spec] .= p.Vs[spec].(t, x.x[spec])
        for other in 1:N, i in eachindex(x.x[spec])
            dx.x[spec][i] += sampled_interaction(t, x.x[spec][i], p.Wprimes[spec][other], x.x[other])
        end
        # for other in 1:N
        #     dx.x[spec] .+= sampled_interaction.(t, x.x[spec]; Wprime=p.Wprimes[spec][other], particles=x.x[other])
        # end
        d = dens[spec]
        for i in eachindex(dx.x[spec])
            if dx.x[spec][i] < 0
                mob = p.mobilities[spec](d[:, 1, i]...)
            else
                mob = p.mobilities[spec](d[:, 2, i]...)
            end
            dx.x[spec][i] *= mob
        end
    end
end

function velocities!(
    dx::ArrayPartition{F, T},
    x::ArrayPartition{F, T},
    p::IntegratedModel{N, TVs, TWs, Tmobilities},
    t
) where {
    F,
    T <: Tuple{Vararg{AbstractVector{<:Real}}},
    N, TVs, TWs, Tmobilities
}
    dens = pwc_densities(x.x)
    dens_diff = similar(x)
    for s in 1:N
        dens_diff.x[s] .= dens[s][s, 2, :] .- dens[s][s, 1, :]
    end
    for spec in 1:N
        dx.x[spec] .= p.Vs[spec].(t, x.x[spec])
        for other in 1:N, i in eachindex(x.x[spec])
            dx.x[spec][i] += integrated_interaction(t, x.x[spec][i], p.Ws[spec][other], x.x[other], dens_diff.x[other])
        end
        d = dens[spec]
        for i in eachindex(dx.x[spec])
            if dx.x[spec][i] < 0
                mob = p.mobilities[spec](d[:, 1, i]...)
            else
                mob = p.mobilities[spec](d[:, 2, i]...)
            end
            dx.x[spec][i] *= mob
        end
    end
end

function velocities_gen(
    x::ArrayPartition{F, T},
    p::Union{SampledModel, IntegratedModel},
    t
) where {
    F,
    T <: Tuple{Vararg{AbstractVector{<:Real}}}
}
    dx = similar(x)
    velocities_gen!(dx, x, p, t)
    dx
end

@generated function velocities_gen!(
    dx::ArrayPartition{F, T},
    x::ArrayPartition{F, T},
    p::SampledModel{N, TVs, TWprimes, Tmobilities},
    t
) where {
    F,
    T <: Tuple{Vararg{AbstractVector{<:Real}}},
    N, TVs, TWprimes, Tmobilities
}
quote
    dens = pwc_densities(x.x)
    $((quote
        dx.x[$spec] .= p.Vs[$spec].(t, x.x[$spec])
        $((quote
            for i in eachindex(x.x[$spec])
                dx.x[$spec][i] += sampled_interaction(t, x.x[$spec][i], p.Wprimes[$spec][$other], x.x[$other])
            end
        end for other in 1:N)...)
        d = dens[$spec]
        for i in eachindex(dx.x[$spec])
            if dx.x[$spec][i] < 0
                mob = p.mobilities[$spec]($((:(d[$j, 1, i]) for j in 1:N)...))
            else
                mob = p.mobilities[$spec]($((:(d[$j, 2, i]) for j in 1:N)...))
            end
            dx.x[$spec][i] *= mob
        end
    end for spec in 1:N)...)
end
end

@generated function velocities_gen!(
    dx::ArrayPartition{F, T},
    x::ArrayPartition{F, T},
    p::IntegratedModel{N, TVs, TWs, Tmobilities},
    t
) where {
    F,
    T <: Tuple{Vararg{AbstractVector{<:Real}}},
    N, TVs, TWs, Tmobilities
}
quote
    dens = pwc_densities(x.x)
    dens_diff = similar(x)
    for s in 1:N
        dens_diff.x[s] .= dens[s][s, 2, :] .- dens[s][s, 1, :]
    end
    # $((:(dens_diff.x[$s] .= dens[$s][$s, 1, :] .- dens[$s][$s, 2, :]) for s in 1:N)...)
    $((quote
        dx.x[$spec] .= p.Vs[$spec].(t, x.x[$spec])
        $((quote
            for i in eachindex(x.x[$spec])
                dx.x[$spec][i] += integrated_interaction(t, x.x[$spec][i], p.Ws[$spec][$other], x.x[$other], dens_diff.x[$other])
            end
        end for other in 1:N)...)
        d = dens[$spec]
        for i in eachindex(dx.x[$spec])
            if dx.x[$spec][i] < 0
                mob = p.mobilities[$spec]($((:(d[$j, 1, i]) for j in 1:N)...))
            else
                mob = p.mobilities[$spec]($((:(d[$j, 2, i]) for j in 1:N)...))
            end
            dx.x[$spec][i] *= mob
        end
    end for spec in 1:N)...)
end
end

"""
    make_velocity(V::Function, Wprime::Function, mobility::Function)

Creates a function `velocity(dx, x, p, t)` that computes the velocity of the particles
under the influence of an external velocity `V`, a mutual interaction `Wprime` and the
congestion given by `mobility`.

See also [`make_velocities`](@ref).
"""
function make_velocity(V::Function, Wprime::Function, mobility::Function)
    function velocity(dx, x, p, t)
        R = pwc_density(x)
        len = length(x)
        for i in 1:len
            v::eltype(x) = 0
            for j in 1:i-1
                v += Wprime(x[j] - x[i])
            end
            for j in i+1:len
                v += Wprime(x[j] - x[i])
            end
            v /= len - 1
            v += V(x[i])
            if v < 0
                mob = mobility(R[i])
            else
                mob = mobility(R[i+1])
            end
            dx[i] = v * mob
        end
    end
    velocity
end

"""
    make_velocities((V₁, ...), ((W′₁₁, ...), ...), (mob₁, ...))

Creates a function `velocities(dx, x, p, t)` that computes the velocity of the particles
under the influence of the external velocities `Vᵢ`, mutual interactions `W′ᵢⱼ` and
congestions given by `mobᵢ`.

See also [`make_velocity`](@ref).
"""
function make_velocities(
        Vs::Tuple{Vararg{Any,N}},
        Wprimes::Tuple{Vararg{Tuple{Vararg{Any,N}},N}},
        mobilities::Tuple{Vararg{Any,N}}
        ) where N
    function velocities(
            dx::ArrayPartition{F, T},
            x::ArrayPartition{F, T},
            p,
            t
            ) where F where T<:Tuple{Vararg{AbstractVector{<:Real}}}
        dens = pwc_densities(x.x)
        for spec in 1:N
            for i in 1:length(x.x[spec])
                v::F = Vs[spec](t, x.x[spec][i])
                for other in 1:N
                    v += sampled_interaction(t, x.x[spec][i], Wprimes[spec][other], x.x[other])
                end
                if v < 0
                    mob = mobilities[spec](dens[spec][:, 1, i]...)
                else
                    mob = mobilities[spec](dens[spec][:, 2, i]...)
                end
                dx.x[spec][i] = v * mob
            end
        end
    end
    velocities
end

export abstract_velocities, abstract_velocities!
export abstract_velocities_gen, abstract_velocities_gen!

function apply_velocity!(dx, x, p::AbstractModel, t, s::Integer)
    xs = species(x, s)
    dxs = species(dx, s)
    dxs .= external_velocity(p, s).(t, xs)
end

function apply_all_interactions!(dx, x, p::AbstractModel, t, s::Integer, dens_diff)
    for o in eachspecies(p)
        apply_interaction!(dx, x, p, t, s, o, dens_diff.x[o])
    end
end

function apply_interaction!(dx, x, p::Union{SampledModel,DiffusiveSampledModel}, t, s::Integer, o::Integer, dens_diff)
    xs = species(x, s)
    dxs = species(dx, s)
    for i in eachindex(xs)
        dxs[i] += sampled_interaction(t, xs[i], interaction(p, s, o), species(x, o))
    end
end

function apply_interaction!(dx, x, p::Union{IntegratedModel,DiffusiveIntegratedModel}, t, s::Integer, o::Integer, dens_diff)
    xs = species(x, s)
    dxs = species(dx, s)
    for i in eachindex(xs)
        dxs[i] += integrated_interaction(t, xs[i], interaction(p, s, o), species(x, o), dens_diff)
    end
end

function apply_interaction!(dx, x, p::Union{HyperbolicModel,ParabolicModel}, t, s::Integer, o::Integer, dens_diff)
    xs = species(x, s)
    dxs = species(dx, s)
    for i in eachindex(xs)
        dxs[i] += compute_interaction(t, xs[i], interaction(p, s, o), species(x, o), dens_diff)
    end
end

function apply_diffusion!(dx, x, p::AbstractModel, t, s::Integer, dens)
    diffuse!(species(dx, s), species(x, s), dens, diffusion(p, s))
end

function abstract_velocities(
    x::ArrayPartition{F, T},
    p::AbstractModel,
    t
) where {
    F,
    T <: Tuple{Vararg{AbstractVector{<:Real}}}
}
    dx = similar(x)
    abstract_velocities!(dx, x, p, t)
    dx
end

function abstract_velocities!(
    dx::ArrayPartition{F, T},
    x::ArrayPartition{F, T},
    p::AbstractModel,
    t
) where {
    F,
    T <: Tuple{Vararg{AbstractVector{<:Real}}},
}
    dens = pwc_densities(x.x)
    dens_diff = similar(x)
    for s in eachspecies(p)
        dens_diff.x[s] .= dens[s][s, 2, :] .- dens[s][s, 1, :]
    end
    for s in eachspecies(p)
        xs = species(x, s)
        dxs = species(dx, s)
        apply_velocity!(dx, x, p, t, s)
        apply_all_interactions!(dx, x, p, t, s, dens_diff)
        d = dens[s]
        mob = mobility(p, s)
        for i in eachindex(dxs)
            if dxs[i] < 0
                m = mob(d[:, 1, i]...)
            else
                m = mob(d[:, 2, i]...)
            end
            dxs[i] *= m
        end
        apply_diffusion!(dx, x, p, t, s, pwc_density(xs))
    end
end

function abstract_velocities!(
    dx::ArrayPartition{F, T},
    x::ArrayPartition{F, T},
    p::HyperbolicModel,
    t
) where {
    F,
    T <: Tuple{Vararg{AbstractVector{<:Real}}},
}
    dens = pwc_densities(x.x)
    dens_diff = similar(x)
    for s in eachspecies(p)
        dens_diff.x[s] .= dens[s][s, 2, :] .- dens[s][s, 1, :]
    end
    for s in eachspecies(p)
        xs = species(x, s)
        dxs = species(dx, s)
        apply_velocity!(dx, x, p, t, s)
        apply_all_interactions!(dx, x, p, t, s, dens_diff)
        d = dens[s]
        mob = mobility(p, s)
        for i in eachindex(dxs)
            if dxs[i] < 0
                m = mob(d[:, 1, i]...)
            else
                m = mob(d[:, 2, i]...)
            end
            dxs[i] *= m
        end
        apply_diffusion!(dx, x, p, t, s, pwc_density(xs))
    end
end

function abstract_velocities_gen(
    x::ArrayPartition{F, T},
    p::AbstractModel,
    t
) where {
    F,
    T <: Tuple{Vararg{AbstractVector{<:Real}}}
}
    dx = similar(x)
    abstract_velocities_gen!(dx, x, p, t)
    dx
end

@generated function abstract_velocities_gen!(
    dx::ArrayPartition{F, T},
    x::ArrayPartition{F, T},
    p::AbstractModel,
    t
) where {
    F,
    T <: Tuple{Vararg{AbstractVector{<:Real}}},
}
quote
    dens = pwc_densities(x.x)
    dens_diff = similar(x)
    for s in eachspecies(p)
        dens_diff.x[s] .= dens[s][s, 2, :] .- dens[s][s, 1, :]
    end
    $((quote
        xs = species(x, $s)
        dxs = species(dx, $s)
        apply_velocity!(dx, x, p, t, $s)
        $((quote
            apply_interaction!(dx, x, p, t, $s, $o, dens_diff.x[$o])
        end for o in eachspecies(p))...)
        d = dens[$s]
        mob = mobility(p, $s)
        for i in eachindex(dxs)
            if dxs[i] < 0
                m = mob($((:(d[$j, 1, i]) for j in eachspecies(p))...))
            else
                m = mob($((:(d[$j, 2, i]) for j in eachspecies(p))...))
            end
            dxs[i] *= m
        end
        apply_diffusion!(dx, x, p, t, $s, pwc_density(xs))
    end for s in eachspecies(p))...)
end
end

export velocities_diff, velocities_diff!

function velocities_diff(
    x::ArrayPartition{F, T},
    p::Union{DiffusiveSampledModel, DiffusiveIntegratedModel},
    t
) where {
    F,
    T <: Tuple{Vararg{AbstractVector{<:Real}}}
}
    dx = similar(x)
    velocities_diff!(dx, x, p, t)
    dx
end

function velocities_diff!(
    dx::ArrayPartition{F, T},
    x::ArrayPartition{F, T},
    p::DiffusiveSampledModel{N, TVs, TWprimes, Tmobilities, Tdiffusions},
    t
) where {
    F,
    T <: Tuple{Vararg{AbstractVector{<:Real}}},
    N, TVs, TWprimes, Tmobilities, Tdiffusions
}
    dens = pwc_densities(x.x)
    for spec in 1:N
        dx.x[spec] .= p.Vs[spec].(t, x.x[spec])
        for other in 1:N, i in eachindex(x.x[spec])
            dx.x[spec][i] += sampled_interaction(t, x.x[spec][i], p.Wprimes[spec][other], x.x[other])
        end
        d = dens[spec]
        for i in eachindex(dx.x[spec])
            if dx.x[spec][i] < 0
                mob = p.mobilities[spec](d[:, 1, i]...)
            else
                mob = p.mobilities[spec](d[:, 2, i]...)
            end
            dx.x[spec][i] *= mob
        end
        diffuse!(dx.x[spec], x.x[spec], pwc_density(x.x[spec]), p.diffusions[spec])
    end
end

function velocities_diff!(
    dx::ArrayPartition{F, T},
    x::ArrayPartition{F, T},
    p::DiffusiveIntegratedModel{N, TVs, TWs, Tmobilities, Tdiffusions},
    t
) where {
    F,
    T <: Tuple{Vararg{AbstractVector{<:Real}}},
    N, TVs, TWs, Tmobilities, Tdiffusions
}
    dens = pwc_densities(x.x)
    dens_diff = similar(x)
    for s in 1:N
        dens_diff.x[s] .= dens[s][s, 2, :] .- dens[s][s, 1, :]
    end
    for spec in 1:N
        dx.x[spec] .= p.Vs[spec].(t, x.x[spec])
        for other in 1:N, i in eachindex(x.x[spec])
            dx.x[spec][i] += integrated_interaction(t, x.x[spec][i], p.Ws[spec][other], x.x[other], dens_diff.x[other])
        end
        d = dens[spec]
        for i in eachindex(dx.x[spec])
            if dx.x[spec][i] < 0
                mob = p.mobilities[spec](d[:, 1, i]...)
            else
                mob = p.mobilities[spec](d[:, 2, i]...)
            end
            dx.x[spec][i] *= mob
        end
        diffuse!(dx.x[spec], x.x[spec], pwc_density(x.x[spec]), p.diffusions[spec])
    end
end


text:
module testCI

"""
	lol()

Demo-function to test documentation.

### Examples
```julia
julia> lol(5)
5
```
"""
function lol(x)
	return x
end

"""
	lol2()

Says hello as our new toolbox.

### Examples
```julia
julia> lol2()=="hallo"
true
```
"""
function lol2()
	return "hallo"
end

end

ext: jl

max_stars_repo_name: renerichter/testCI

max_issues_repo_name: renerichter/testCI

max_forks_repo_name: renerichter/testCI


text:
@testset "SDP and Exp Atoms: $solver" for solver in solvers
    if can_solve_sdp(solver) && can_solve_exp(solver)
        tol = 1e-2
        @testset "log det atom" begin
            x = Variable(2, 2)
            p = maximize(logdet(x), [x[1, 1] == 1, x[2, 2] == 1])
            @test vexity(p) == ConvexVexity()
            solve!(p, solver)
            @test p.optval ≈ 0 atol=tol
            @test evaluate(logdet(x)) ≈ 0 atol=tol
        end
    end
end


text:
module TwodWrapper

using JavaCall

const AbstractRegion = @jimport org.hipparchus.geometry.partitioning.AbstractRegion
const AbstractSubHyperplane = @jimport org.hipparchus.geometry.partitioning.AbstractSubHyperplane
const BSPTree = @jimport org.hipparchus.geometry.partitioning.BSPTree
const BoundaryProjection = @jimport org.hipparchus.geometry.partitioning.BoundaryProjection
const Class = @jimport java.lang.Class
const Collection = @jimport java.util.Collection
const DiskGenerator = @jimport org.hipparchus.geometry.euclidean.twod.DiskGenerator
const EnclosingBall = @jimport org.hipparchus.geometry.enclosing.EnclosingBall
const Euclidean1D = @jimport org.hipparchus.geometry.euclidean.oned.Euclidean1D
const Euclidean2D = @jimport org.hipparchus.geometry.euclidean.twod.Euclidean2D
const Field = @jimport org.hipparchus.Field
const FieldPosition = @jimport java.text.FieldPosition
const FieldVector2D = @jimport org.hipparchus.geometry.euclidean.twod.FieldVector2D
const Hyperplane = @jimport org.hipparchus.geometry.partitioning.Hyperplane
const JString = @jimport java.lang.String
const JVector = @jimport org.hipparchus.geometry.Vector
const Line = @jimport org.hipparchus.geometry.euclidean.twod.Line
const List = @jimport java.util.List
const Locale = @jimport java.util.Locale
const NumberFormat = @jimport java.text.NumberFormat
const Object = @jimport java.lang.Object
const ParsePosition = @jimport java.text.ParsePosition
const Point = @jimport org.hipparchus.geometry.Point
const PolygonsSet = @jimport org.hipparchus.geometry.euclidean.twod.PolygonsSet
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const Region = @jimport org.hipparchus.geometry.partitioning.Region
const Region_Location = @jimport org.hipparchus.geometry.partitioning.Region$Location
const Segment = @jimport org.hipparchus.geometry.euclidean.twod.Segment
const Space = @jimport org.hipparchus.geometry.Space
const StringBuffer = @jimport java.lang.StringBuffer
const SubHyperplane = @jimport org.hipparchus.geometry.partitioning.SubHyperplane
const SubHyperplane_SplitSubHyperplane = @jimport org.hipparchus.geometry.partitioning.SubHyperplane$SplitSubHyperplane
const SubLine = @jimport org.hipparchus.geometry.euclidean.twod.SubLine
const Transform = @jimport org.hipparchus.geometry.partitioning.Transform
const Vector1D = @jimport org.hipparchus.geometry.euclidean.oned.Vector1D
const Vector2D = @jimport org.hipparchus.geometry.euclidean.twod.Vector2D
const Vector2DFormat = @jimport org.hipparchus.geometry.euclidean.twod.Vector2DFormat
const VectorFormat = @jimport org.hipparchus.geometry.VectorFormat

export DiskGenerator
export Euclidean2D
export FieldVector2D
export Line
export PolygonsSet
export Segment
export SubLine
export Vector2D
export Vector2DFormat
export add
export angle
export apply_transform
export ball_on_support
export build_new
export check_point
export contains
export copy_self
export cross_product
export distance
export distance1
export distance_inf
export distance_sq
export dot_product
export empty_hyperplane
export equals
export format
export get_angle
export get_available_locales
export get_barycenter
export get_boundary_size
export get_class
export get_dimension
export get_end
export get_format
export get_hyperplane
export get_instance
export get_length
export get_line
export get_minus_i
export get_minus_j
export get_na_n
export get_negative_infinity
export get_norm
export get_norm1
export get_norm_inf
export get_norm_sq
export get_offset
export get_origin_offset
export get_plus_i
export get_plus_j
export get_point_at
export get_positive_infinity
export get_prefix
export get_remaining_region
export get_reverse
export get_segments
export get_separator
export get_size
export get_space
export get_start
export get_sub_space
export get_suffix
export get_tolerance
export get_transform
export get_tree
export get_vector2d_format
export get_vertices
export get_x
export get_y
export get_zero
export hash_code
export intersection
export is_empty
export is_full
export is_infinite
export is_na_n
export is_parallel_to
export negate
export normalize
export notify
export notify_all
export orientation
export parse
export project
export project_to_boundary
export reset
export reunite
export revert_self
export same_orientation_as
export scalar_multiply
export set_angle
export set_origin_offset
export split
export subtract
export to_array
export to_space
export to_string
export to_sub_space
export to_vector2d
export translate_to_point
export wait
export whole_hyperplane
export whole_space

include("HullWrapper/HullWrapper.jl")
include("disk_generator.jl")
include("euclidean2_d.jl")
include("field_vector2_d.jl")
include("line.jl")
include("polygons_set.jl")
include("segment.jl")
include("sub_line.jl")
include("vector2_d.jl")
include("vector2_d_format.jl")

end


text:
using PEG
using Test

include("./calculator.jl")
include("./meta.jl")
include("./misc.jl")


text:
function make_complete!(graph, edges)
    for edge in combinations(edges, 2)
        add_edge!(graph, edge...)
    end
end

function getsaturate(g, v)::Vector{Vector{Int}}
    return collect(combinations(neighbors(g, v), 2))
end

function mindeg(g)
    return argmin(degree(g))
end

function add_edges!(g, edges::Vector{Vector{T}}) where T
    added = Vector{Bool}(undef, 0)
    for e in edges
        push!(added, add_edge!(g, e[1], e[2]))
    end
    return added
end

function vertices_meta(g::MetaGraph)
    return Iterators.Stateful([get_prop(g, v, :v) for v in vertices(g)])
end

function edges_meta(g::MetaGraph)
    return Iterators.Stateful([edges_m2s(e) for e in edges(g)])
end

function vertex_m2s(g::MetaGraph, v)
    return get_prop(g, v, :v)
end

function vertex_s2m(g::MetaGraph, v)
    # Here the index are integers so g[v, :v] wouldn't return the node index
    # There seems to be no function to retrieve the index of an integer indexed prop
    # so I access the field `metaindex` directly
    return g.metaindex[:v][v]
end

function edge_m2s(g::MetaGraph, src, dst)
    return [get_prop(g, src, :v), get_prop(g, dst, :v)]
end

function edge_s2m(g::MetaGraph, src, dst)
    return [(g[src, :v], g[dst, :v])]
end

function min_local_clustering(g::AbstractGraph)
    return argmin(local_clustering_coefficient(g))
end

function max_local_clustering(g::AbstractGraph)
    return argmax(local_clustering_coefficient(g))
end

function min_triangles(g::AbstractGraph)
    return argmin(triangles(g))
end

function max_triangles(g::AbstractGraph)
    return argmax(triangles(g))
end

function edge_distance(g1::AbstractGraph, g2::AbstractGraph)
    vertices(g1) != vertices(g2) && error("`edge_distance` only works on graphs with equal vertices set")
    diff_edges = union(setdiff(edges(g1), edges(g2)), setdiff(edges(g2), edges(g1)))
    return length(diff_edges)
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
#  Internal errors to rdkafka:
#  Begin internal error codes
const RD_KAFKA_RESP_ERR__BEGIN = -200
#  Received message is incorrect
const RD_KAFKA_RESP_ERR__BAD_MSG = -199
#  Bad/unknown compression
const RD_KAFKA_RESP_ERR__BAD_COMPRESSION = -198
#  Broker is going away
const RD_KAFKA_RESP_ERR__DESTROY = -197
#  Generic failure
const RD_KAFKA_RESP_ERR__FAIL = -196
#  Broker transport failure
const RD_KAFKA_RESP_ERR__TRANSPORT = -195
#  Critical system resource
const RD_KAFKA_RESP_ERR__CRIT_SYS_RESOURCE = -194
#  Failed to resolve broker
const RD_KAFKA_RESP_ERR__RESOLVE = -193
#  Produced message timed out*/
const RD_KAFKA_RESP_ERR__MSG_TIMED_OUT = -192
#  Reached the end of the topic+partition queue on
#  the broker. Not really an error.
const RD_KAFKA_RESP_ERR__PARTITION_EOF = -191
#  Permanent: Partition does not exist in cluster.
const RD_KAFKA_RESP_ERR__UNKNOWN_PARTITION = -190
#  File or filesystem error
const RD_KAFKA_RESP_ERR__FS = -189
#  Permanent: Topic does not exist in cluster.
const RD_KAFKA_RESP_ERR__UNKNOWN_TOPIC = -188
#  All broker connections are down.
const RD_KAFKA_RESP_ERR__ALL_BROKERS_DOWN = -187
#  Invalid argument, or invalid configuration
const RD_KAFKA_RESP_ERR__INVALID_ARG = -186
#  Operation timed out
const RD_KAFKA_RESP_ERR__TIMED_OUT = -185
#  Queue is full
const RD_KAFKA_RESP_ERR__QUEUE_FULL = -184
#  ISR count < required.acks
const RD_KAFKA_RESP_ERR__ISR_INSUFF = -183
#  Broker node update
const RD_KAFKA_RESP_ERR__NODE_UPDATE = -182
#  SSL error
const RD_KAFKA_RESP_ERR__SSL = -181
#  Waiting for coordinator to become available.
const RD_KAFKA_RESP_ERR__WAIT_COORD = -180
#  Unknown client group
const RD_KAFKA_RESP_ERR__UNKNOWN_GROUP = -179
#  Operation in progress
const RD_KAFKA_RESP_ERR__IN_PROGRESS = -178
#  Previous operation in progress, wait for it to finish.
const RD_KAFKA_RESP_ERR__PREV_IN_PROGRESS = -177
#  This operation would interfere with an existing subscription
const RD_KAFKA_RESP_ERR__EXISTING_SUBSCRIPTION = -176
#  Assigned partitions (rebalance_cb)
const RD_KAFKA_RESP_ERR__ASSIGN_PARTITIONS = -175
#  Revoked partitions (rebalance_cb)
const RD_KAFKA_RESP_ERR__REVOKE_PARTITIONS = -174
#  Conflicting use
const RD_KAFKA_RESP_ERR__CONFLICT = -173
#  Wrong state
const RD_KAFKA_RESP_ERR__STATE = -172
#  Unknown protocol
const RD_KAFKA_RESP_ERR__UNKNOWN_PROTOCOL = -171
#  Not implemented
const RD_KAFKA_RESP_ERR__NOT_IMPLEMENTED = -170
#  Authentication failure*/
const RD_KAFKA_RESP_ERR__AUTHENTICATION = -169
#  No stored offset
const RD_KAFKA_RESP_ERR__NO_OFFSET = -168
#  Outdated
const RD_KAFKA_RESP_ERR__OUTDATED = -167
#  Timed out in queue
const RD_KAFKA_RESP_ERR__TIMED_OUT_QUEUE = -166
#  Feature not supported by broker
const RD_KAFKA_RESP_ERR__UNSUPPORTED_FEATURE = -165
#  Awaiting cache update
const RD_KAFKA_RESP_ERR__WAIT_CACHE = -164
#  Operation interrupted (e.g., due to yield))
const RD_KAFKA_RESP_ERR__INTR = -163
#  Key serialization error
const RD_KAFKA_RESP_ERR__KEY_SERIALIZATION = -162
#  Value serialization error
const RD_KAFKA_RESP_ERR__VALUE_SERIALIZATION = -161
#  Key deserialization error
const RD_KAFKA_RESP_ERR__KEY_DESERIALIZATION = -160
#  Value deserialization error
const RD_KAFKA_RESP_ERR__VALUE_DESERIALIZATION = -159
#  Partial response
const RD_KAFKA_RESP_ERR__PARTIAL = -158
#  Modification attempted on read-only object
const RD_KAFKA_RESP_ERR__READ_ONLY = -157
#  No such entry / item not found
const RD_KAFKA_RESP_ERR__NOENT = -156
#  Read underflow
const RD_KAFKA_RESP_ERR__UNDERFLOW = -155

#  End internal error codes
const RD_KAFKA_RESP_ERR__END = -100

##  Kafka broker errors:
#  Unknown broker error
const RD_KAFKA_RESP_ERR_UNKNOWN = -1
#  Success
const RD_KAFKA_RESP_ERR_NO_ERROR = 0
#  Offset out of range
const RD_KAFKA_RESP_ERR_OFFSET_OUT_OF_RANGE = 1
#  Invalid message
const RD_KAFKA_RESP_ERR_INVALID_MSG = 2
#  Unknown topic or partition
const RD_KAFKA_RESP_ERR_UNKNOWN_TOPIC_OR_PART = 3
#  Invalid message size
const RD_KAFKA_RESP_ERR_INVALID_MSG_SIZE = 4
#  Leader not available
const RD_KAFKA_RESP_ERR_LEADER_NOT_AVAILABLE = 5
#  Not leader for partition
const RD_KAFKA_RESP_ERR_NOT_LEADER_FOR_PARTITION = 6
#  Request timed out
const RD_KAFKA_RESP_ERR_REQUEST_TIMED_OUT = 7
#  Broker not available
const RD_KAFKA_RESP_ERR_BROKER_NOT_AVAILABLE = 8
#  Replica not available
const RD_KAFKA_RESP_ERR_REPLICA_NOT_AVAILABLE = 9
#  Message size too large
const RD_KAFKA_RESP_ERR_MSG_SIZE_TOO_LARGE = 10
#  StaleControllerEpochCode
const RD_KAFKA_RESP_ERR_STALE_CTRL_EPOCH = 11
#  Offset metadata string too large
const RD_KAFKA_RESP_ERR_OFFSET_METADATA_TOO_LARGE = 12
#  Broker disconnected before response received
const RD_KAFKA_RESP_ERR_NETWORK_EXCEPTION = 13
#  Group coordinator load in progress
const RD_KAFKA_RESP_ERR_GROUP_LOAD_IN_PROGRESS = 14
#  Group coordinator not available
const RD_KAFKA_RESP_ERR_GROUP_COORDINATOR_NOT_AVAILABLE = 15
#  Not coordinator for group
const RD_KAFKA_RESP_ERR_NOT_COORDINATOR_FOR_GROUP = 16
#  Invalid topic
const RD_KAFKA_RESP_ERR_TOPIC_EXCEPTION = 17
#  Message batch larger than configured server segment size
const RD_KAFKA_RESP_ERR_RECORD_LIST_TOO_LARGE = 18
#  Not enough in-sync replicas
const RD_KAFKA_RESP_ERR_NOT_ENOUGH_REPLICAS = 19
#  Message(s) written to insufficient number of in-sync replicas
const RD_KAFKA_RESP_ERR_NOT_ENOUGH_REPLICAS_AFTER_APPEND = 20
#  Invalid required acks value
const RD_KAFKA_RESP_ERR_INVALID_REQUIRED_ACKS = 21
#  Specified group generation id is not valid
const RD_KAFKA_RESP_ERR_ILLEGAL_GENERATION = 22
#  Inconsistent group protocol
const RD_KAFKA_RESP_ERR_INCONSISTENT_GROUP_PROTOCOL = 23
#  Invalid group.id
const RD_KAFKA_RESP_ERR_INVALID_GROUP_ID = 24
#  Unknown member
const RD_KAFKA_RESP_ERR_UNKNOWN_MEMBER_ID = 25
#  Invalid session timeout
const RD_KAFKA_RESP_ERR_INVALID_SESSION_TIMEOUT = 26
#  Group rebalance in progress
const RD_KAFKA_RESP_ERR_REBALANCE_IN_PROGRESS = 27
#  Commit offset data size is not valid
const RD_KAFKA_RESP_ERR_INVALID_COMMIT_OFFSET_SIZE = 28
#  Topic authorization failed
const RD_KAFKA_RESP_ERR_TOPIC_AUTHORIZATION_FAILED = 29
#  Group authorization failed
const RD_KAFKA_RESP_ERR_GROUP_AUTHORIZATION_FAILED = 30
#  Cluster authorization failed
const RD_KAFKA_RESP_ERR_CLUSTER_AUTHORIZATION_FAILED = 31
#  Invalid timestamp
const RD_KAFKA_RESP_ERR_INVALID_TIMESTAMP = 32
#  Unsupported SASL mechanism
const RD_KAFKA_RESP_ERR_UNSUPPORTED_SASL_MECHANISM = 33
#  Illegal SASL state
const RD_KAFKA_RESP_ERR_ILLEGAL_SASL_STATE = 34
#  Unuspported version
const RD_KAFKA_RESP_ERR_UNSUPPORTED_VERSION = 35
#  Topic already exists
const RD_KAFKA_RESP_ERR_TOPIC_ALREADY_EXISTS = 36
#  Invalid number of partitions
const RD_KAFKA_RESP_ERR_INVALID_PARTITIONS = 37
#  Invalid replication factor
const RD_KAFKA_RESP_ERR_INVALID_REPLICATION_FACTOR = 38
#  Invalid replica assignment
const RD_KAFKA_RESP_ERR_INVALID_REPLICA_ASSIGNMENT = 39
#  Invalid config
const RD_KAFKA_RESP_ERR_INVALID_CONFIG = 40
#  Not controller for cluster
const RD_KAFKA_RESP_ERR_NOT_CONTROLLER = 41
#  Invalid request
const RD_KAFKA_RESP_ERR_INVALID_REQUEST = 42
#  Message format on broker does not support request
const RD_KAFKA_RESP_ERR_UNSUPPORTED_FOR_MESSAGE_FORMAT = 43
#  Isolation policy volation
const RD_KAFKA_RESP_ERR_POLICY_VIOLATION = 44
#  Broker received an out of order sequence number
const RD_KAFKA_RESP_ERR_OUT_OF_ORDER_SEQUENCE_NUMBER = 45
#  Broker received a duplicate sequence number
const RD_KAFKA_RESP_ERR_DUPLICATE_SEQUENCE_NUMBER = 46
#  Producer attempted an operation with an old epoch
const RD_KAFKA_RESP_ERR_INVALID_PRODUCER_EPOCH = 47
#  Producer attempted a transactional operation in an invalid state
const RD_KAFKA_RESP_ERR_INVALID_TXN_STATE = 48
#  Producer attempted to use a producer id which is not
#  currently assigned to its transactional id
const RD_KAFKA_RESP_ERR_INVALID_PRODUCER_ID_MAPPING = 49
#  Transaction timeout is larger than the maximum
#  value allowed by the broker's max.transaction.timeout.ms
const RD_KAFKA_RESP_ERR_INVALID_TRANSACTION_TIMEOUT = 50
#  Producer attempted to update a transaction while another
#  concurrent operation on the same transaction was ongoing
const RD_KAFKA_RESP_ERR_CONCURRENT_TRANSACTIONS = 51
#  Indicates that the transaction coordinator sending a
#  WriteTxnMarker is no longer the current coordinator for a
#  given producer
const RD_KAFKA_RESP_ERR_TRANSACTION_COORDINATOR_FENCED = 52
#  Transactional Id authorization failed
const RD_KAFKA_RESP_ERR_TRANSACTIONAL_ID_AUTHORIZATION_FAILED = 53
#  Security features are disabled
const RD_KAFKA_RESP_ERR_SECURITY_DISABLED = 54
#  Operation not attempted
const RD_KAFKA_RESP_ERR_OPERATION_NOT_ATTEMPTED = 55


text:
# indicator of nonnegative orthant

"""
  IndNonnegative()

Returns the indicator function the nonnegative orthant, that is

  `g(x) = 0 if x ⩾ 0, +∞ otherwise`
"""

immutable IndNonnegative <: IndicatorConvexCone end

function (f::IndNonnegative){R <: Real}(x::AbstractArray{R})
  for k in eachindex(x)
    if x[k] < 0
      return +Inf
    end
  end
  return 0.0
end

function prox!{R <: Real}(f::IndNonnegative, x::AbstractArray{R}, y::AbstractArray{R}, gamma::Real=1.0)
  for k in eachindex(x)
    if x[k] < 0
      y[k] = zero(R)
    else
      y[k] = x[k]
    end
  end
  return 0.0
end

fun_name(f::IndNonnegative) = "indicator of the Nonnegative cone"
fun_dom(f::IndNonnegative) = "AbstractArray{Real}"
fun_expr(f::IndNonnegative) = "x ↦ 0 if all(0 ⩽ x), +∞ otherwise"
fun_params(f::IndNonnegative) = "none"

function prox_naive{R <: Real}(f::IndNonnegative, x::AbstractArray{R}, gamma::Real=1.0)
  y = max.(zero(R), x)
  return y, 0.0
end


text:
function remove_char(s::String, char::String)
    return replace(s, Regex("[$char]") => "")
end

function prepare_workspace!(folder::String)

    !isdir(folder) && error("Specified folder is not valid")

    cd(folder)
    global_path = joinpath(folder)
    isdir(global_path) && mkpath(global_path)
    _sim_path =  remove_char("$(round(Dates.now(),Dates.Minute))", ":")
    simulation_path = joinpath(global_path,)
    raw_ouput = joinpath(simulation_path)
    mkpath(raw_ouput)

    return raw_ouput

end


text:
"""

```
wget(PDBid; selection)
```
Retrieves a PDB file from the protein data bank. Selections may be applied.

### Example

```julia-repl
julia> protein = wget("1LBD","chain A")
   Array{Atoms,1} with 1870 atoms with fields:
   index name resname chain   resnum  residue        x        y        z  beta occup model segname index_pdb
       1    N     SER     A      225        1   45.228   84.358   70.638 67.05  1.00     1       -         1
       2   CA     SER     A      225        1   46.080   83.165   70.327 68.73  1.00     1       -         2
       3    C     SER     A      225        1   45.257   81.872   70.236 67.90  1.00     1       -         3
                                                       ⋮
    1868  OG1     THR     A      462      238  -27.462   74.325   48.885 79.98  1.00     1       -      1868
    1869  CG2     THR     A      462      238  -27.063   71.965   49.222 78.62  1.00     1       -      1869
    1870  OXT     THR     A      462      238  -25.379   71.816   51.613 84.35  1.00     1       -      1870

```
"""
function wget(pdb_id::String, selection::String)
    query = parse_query(selection)
    return wget(pdb_id, only = atom -> apply_query(query, atom))
end

function wget(pdb_id::String; only = all)
    file = download("https://files.rcsb.org/download/$(pdb_id).pdb")
    return readPDB(file, only = only)
end

max_stars_repo_name: mcubeg/PDBTools

max_issues_repo_name: mcubeg/PDBTools

max_forks_repo_name: mcubeg/PDBTools


text:
include("stop.jl")
include("eas.jl")
include("simulated-annealing.jl")
export name
export optimize, ea1pλwith2rates, ea1p1, ea1p1SD, ea1pλSASD, RLSSDstar, RLS12, RLSSDstarstar
export fastSD
export SDCounter, SDCounterEstimation, RLSSDCounter, threshold_gen
export FixedBudget, niterations
export FixedCooling, temperature, SimulatedAnnealing, AbstractCooling, FixedCoolingTag

name(algo::AbstractAlgorithm) = algo.name

function Base.show(io::IO, algo::AbstractAlgorithm)
    print(io, algo.name)
    first = true
    for p in deepparameters(algo)
        if first
            print(io, " with ")
            first = false
        else
            print(io, " and ")
        end

        print(io, string(p[1]))
        print(io, ":")
        show(io, p[2])
    end
end

max_stars_repo_name: DTUComputeTONIA/SDRLS

max_issues_repo_name: DTUComputeTONIA/SDRLS

max_forks_repo_name: DTUComputeTONIA/SDRLS


text:
using UncertaintyQuantification, Random
Random.seed!(2547);
N = 200
mu = [0, 0]
sig = [1 -0.55; -0.55 1]
mvnormal = MvNormal(mu, sig)

sam = zeros(200,2)
for i in 1:200 sam[i,:] = rand(mvnormal) end

Random.seed!()

function likelihood(x)
    if x[1] >1 || x[1]<-1
        rv = 0
    else
        rv = 1
        for i in 1:N
            rv = rv/(2*pi*sqrt(1-x[1]^2))*exp(-1/(2*(1-x[1]^2))*(sam[i,1]^2-2*x[1]*
            sam[i,1]*sam[i,2]+sam[i,2]^2))
        end
    end
    return rv
end

function prior(x)
    if x[1] >1 || x[1]<-1
        rv = 0
    else
        rv = 1/(1-x[1]^2)^(3/2)
    end
    return rv
end
priorsample(n) = rand(Uniform(-1,1),(n))

tmcmcsample = tmcmc(likelihood, prior, priorsample, 2000, 0.2)


text:
"""
Actor-Critic learner (core)

Critic implemented by td-lambda
This implementation follows Sutton & Barto,
Chapter 7.7 "Eligibility Traces for Actor-Critic Methods"
The policy update follows the scheme given in chapter 6.6 where
the policy changes inversely proportional to the action selection probability
"""

mutable struct ActorCritic
    nr_states::Int
    nr_actions::Int
    α_critic::Float64
    α_actor::Float64
    α_actor_eff::Float64
    γ::Float64
    λ_critic::Float64
    λ_actor::Float64
    initvalue::Float64
    V::Array{Float64, 1}
    policyparameters::Array{Float64, 2}
    eligTrace_critic::Array{Float64, 1}
    eligTrace_actor::Array{Float64, 2}
    td_error::Array{Float64, 1}
    use_replacing_trace::Bool
    do_debug_print::Bool
    actioncost::Float64
end
function ActorCritic(; nr_states = 10, nr_actions = 4,
    α_critic = 0.1, α_actor = 0.1, γ = .9,
    λ_critic = 0.9, λ_actor = 0.9, initvalue = 0.,
    V = zeros(Float64, nr_states) .+ initvalue,
    policyparameters = zeros(Float64, nr_states, nr_actions),
    eligTrace_critic = zeros(Float64, nr_states),
    eligTrace_actor = zeros(Float64, nr_states, nr_actions),
    use_replacing_trace = true,
    do_debug_print = false,
    actioncost = 0.)

    ActorCritic(nr_states, nr_actions, α_critic, α_actor, α_actor, γ,
        λ_critic, λ_actor, initvalue, V,
        policyparameters, eligTrace_critic, eligTrace_actor, [0.],
        use_replacing_trace, do_debug_print, actioncost)
end
export ActorCritic
function do_debug_print_final(learner::Union{ActorCritic, ActorCriticModulated})
    println("Final:")
    for s in 1:size(learner.V,1)
        @show learner.V[s]
    end
    for s in 1:size(learner.policyparameters,1)
        @show learner.policyparameters[s,:]
    end
    println("--------------")
end

""" Update: main function """
function update!(learner::Union{ActorCritic, ActorCriticModulated},
                state, action, next_state, reward,
                p_action, modulatingfactor)

    if learner.use_replacing_trace # replacing trace
        learner.eligTrace_critic[state] = 1.
        learner.eligTrace_actor[state, action] = (1. - p_action) #see eq. 7.12/7.13 in Sutton & Barto
    else # (accumulating) eligibility trace
        learner.eligTrace[state] += 1.
        learner.eligTrace_actor[state, action] += (1. - p_action) #see eq. 7.12 in Sutton & Barto
    end

    getlearningrate!(learner, modulatingfactor)

    learner.td_error[1] = reward - learner.actioncost + (learner.γ * learner.V[next_state]) - learner.V[state]
    learner.V[:] .+= learner.α_critic  * learner.td_error[1] * learner.eligTrace_critic
    learner.policyparameters[:, :] .+= learner.α_actor_eff * learner.td_error[1] * learner.eligTrace_actor

    if reward>0
        learner.eligTrace_critic = zeros(Float64, learner.nr_states)
        learner.eligTrace_actor = zeros(Float64, learner.nr_states, learner.nr_actions)
    else #non goal state
        learner.eligTrace_actor[:, :] .*= learner.λ_actor * learner.γ
        learner.eligTrace_critic[:] .*= learner.λ_critic * learner.γ
    end
    if learner.do_debug_print
        do_debug_print_final(learner)
    end
end
export update!
function getlearningrate!(learner::ActorCritic, modulatingfactor)
    learner.α_actor_eff = learner.α_actor
end


text:
using FeedbackParticleFilters, Random, Statistics, LinearAlgebra, Distributions

println("Testing FPF.jl:")

@testset "FPF.jl" begin

    print("  inner constructor for FPFState")
    ens   = UnweightedParticleEnsemble(randn(3,10))
    h(x)  = [x[1] * x[2], x[2] * x[3]]
    eq    = PoissonEquation(h, ens)
    state = FPFState(ens, eq)
    print(".")
    @test state.ensemble == ens
    print(".")
    @test state.eq == eq
    println("DONE.")

    print("  inner constructor for FPF")
    f(x)   = -x
    g(x)   = hcat([1.; -2.; 1.])
    init   = [0., 0., 0.]
    st_mod = DiffusionStateModel(f, g, init)
    ob_mod = DiffusionObservationModel{Float64, Float64, typeof(h)}(3, 2, h)
    struct DummyMethod <: GainEstimationMethod end
    method = DummyMethod()
    fpf    = FPF(st_mod, ob_mod, method, 10)
    print(".")
    @test fpf.state_model == st_mod
    print(".")
    @test fpf.obs_model == ob_mod
    print(".")
    @test fpf.gain_method == method
    print(".")
    @test fpf.N == 10

    struct DummyModel <: ObservationModel{Vector{Float64}, Vector{Float64}, ContinuousTime} end
    ob_mod2 = DummyModel()
    print(".")
    @test_throws ErrorException FPF(st_mod, ob_mod2, method, 10) # wrong type of observation model
    println("DONE.")

    print("  method initial_condition")
    print(".")
    @test initial_condition(fpf) ≈ [0., 0., 0.]
    println("DONE.")

    print("  method no_of_particles")
    print(".")
    @test no_of_particles(fpf) == 10
    print(".")
    @test no_of_particles(state) == 10
    println("DONE.")

    print("  method state_model")
    print(".")
    @test state_model(fpf) == st_mod
    println("DONE.")

    print("  method obs_model")
    print(".")
    @test obs_model(fpf) == ob_mod
    println("DONE.")

    print("  method gain_estimation_method")
    print(".")
    @test gain_estimation_method(fpf) == method
    println("DONE.")

    print("  method for Statistics.mean")
    print(".")
    @test Statistics.mean(state) ≈ mean(ens)
    println("DONE.")

    print("  method for Statistics.cov")
    print(".")
    @test Statistics.cov(state) ≈ cov(ens)
    println("DONE.")

    print("  method for Statistics.var")
    print(".")
    @test Statistics.var(state) ≈ var(ens)
    println("DONE.")

    print("  outer constructors for FPFState")
    state2 = FPFState(fpf)
    print(".")
    @test no_of_particles(state2) == 10
    print(".")
    @test state2.ensemble.positions == zeros(3,10)
    filt_prob = FilteringProblem(st_mod, ob_mod)
    state3 = FPFState(filt_prob, 20)
    print(".")
    @test no_of_particles(state3) == 20
    println("DONE.")

    print("  outer constructors for FPF")
    fpf2 = FPF(filt_prob, method, 100)
    print(".")
    @test no_of_particles(fpf2) == 100
    println("DONE.")

    print("  method initialize")
    state4 = initialize(fpf)
    print(".")
    @test state4 isa FPFState
    print(".")
    @test state4.ensemble.positions ≈ zeros(3,10)
    println("DONE.")

    print("  method assimilate!")
    function FeedbackParticleFilters.solve!(eq::PoissonEquation, method::DummyMethod)
        eq.gain .= Statistics.mean(eq.positions) * ones(size(eq.gain))
    end

    init   = MvNormal(3, 1.)
    st_mod = DiffusionStateModel(f, g, init)
    ob_mod = DiffusionObservationModel{Float64, Float64, typeof(h)}(3, 2, h)
    method = DummyMethod()
    fpf    = FPF(st_mod, ob_mod, method, 10)
    state  = initialize(fpf)

    oldstate = deepcopy(state)
    solve!(oldstate.eq, method)
    print(".")
    @test oldstate.eq.gain == Statistics.mean(oldstate.eq.positions) * ones(3,10,2) # test DummyMethod

    assimilate!([0.01, -0.02], state, fpf, 0.01)
    error = [0.01, -0.02] .- 0.01*oldstate.eq.H/2 .- 0.01*oldstate.eq.mean_H/2
    heun!(oldstate.eq, oldstate.ensemble, error, method)
    print(".")
    @test state.eq.gain ≈ oldstate.eq.gain

    print(".")
    @test state.eq.positions ≈ state.ensemble.positions
    println("DONE.")

    print("  method propagate!")
    init   = [0., 0., 0.]
    st_mod = DiffusionStateModel(f, g, init)
    fpf   = FPF(st_mod, ob_mod, method, 2)
    state = initialize(fpf)
    Random.seed!(0)
    propagate!(state, fpf, 0.01)
    print(".")
    @test state.ensemble.positions ≈ [0.06791074260357777 -0.013485387193052173; -0.1656826965800072 -0.11732341492662196; -0.03530074003005963 0.029733585084941616]
    println("DONE.")

    print("  method update!")
    init   = [0., 0., 0.]
    st_mod = DiffusionStateModel(f, g, init)
    fpf    = FPF(st_mod, ob_mod, method, 2)
    state  = initialize(fpf)
    state2 = initialize(fpf)
    Random.seed!(0)
    update!(state, fpf, [0.01, -0.02], 0.01)
    Random.seed!(0)
    propagate!(state2, fpf, 0.01)
    assimilate!([0.01, -0.02], state2, fpf, 0.01)
    print(".")
    @test state.ensemble.positions ≈ state2.ensemble.positions
    print(".")
    @test state.eq.gain ≈ state2.eq.gain
    println("DONE.")

end; #FPF.jl


text:
module SynapseClient

using PyCall

const synapseclient = PyNULL()
const pyclient = PyNULL()
const pyannotations = PyNULL()
const pydict_object = PyNULL()
const pycache = PyNULL()

function __init__()
	try
		copy!(synapseclient, pyimport("synapseclient"))
		copy!(pyclient,      synapseclient.client)
		copy!(pyannotations, synapseclient.annotations)
		copy!(pydict_object, synapseclient.dict_object)
		copy!(pycache,       synapseclient.cache)
		copy!(pyhasattr, pybuiltin(:hasattr))
	catch
		@warn "Please install the python package \"synapseclient\"."
	end
end

import Base: get

export
	Synapse,

	# SynapseClient
	chunkedquery,
	delete,
	deleteprovenance,
	downloadtablecolumns,
	downloadtablefile,
	#get, # NB: it makes more sense to extend Base.get since it is a dictionary lookup
	getannotations,
	getchildren,
	getcolumn,
	getcolumns,
	getconfigfile,
	getevaluation,
	getevaluationbycontentsource,
	getevaluationbyname,
	getpermissions,
	getprovenance,
	getsubmission,
	getsubmissionbundles,
	getsubmissionstatus,
	getsubmissions,
	gettablecolumns,
	getteam,
	getteammembers,
	getuserprofile,
	getwiki,
	getwikiheaders,
	hasattr,
	invalidateapikey,
	login,
	logout,
	md5query,
	onweb,
	printentity,
	restdelete,
	restget,
	restpost,
	restput,
	sendmessage,
	setannotations,
	setendpoints,
	setpermissions,
	setprovenance,
	store,
	submit,
	tablequery,
	updateactivity,
	uploadfile,

	# Entity
	create,

	# Activity
	used,
	executed

include("wrapper.jl")
include("types.jl")
include("synapse.jl")
include("utils.jl")
include("entity.jl")
include("annotations.jl")
include("cache.jl")

end


text:
using ADI
using AstroTime
using CSV
using DataFrames
using LinearAlgebra
using Measurements
using PyCall
using PyPlot
using SiriusB

fits = pyimport("astropy.io.fits")
driver = pyimport("orbitize.driver")
pro = pyimport("proplot")

# rcParams
pro.rc["image.origin"] = "lower"
pro.rc["image.cmap"] = "inferno"
pro.rc["grid"] = false

res_cubes = [
    fits.getdata(datadir("epoch_2020feb04", "residuals", "2020feb04_sirius-b_residual_median.fits")),
    fits.getdata(datadir("epoch_2020nov21", "residuals", "2020nov21_sirius-b_residual_median.fits")),
    fits.getdata(datadir("epoch_2020nov28", "residuals", "2020nov28_sirius-b_residual_annular_pca-2.fits")),
]
angles = [
    fits.getdata(datadir("epoch_2020feb04", "processed", "2020feb04_sirius-b_pa.fits")),
    fits.getdata(datadir("epoch_2020nov21", "processed", "2020nov21_sirius-b_pa.fits")),
    fits.getdata(datadir("epoch_2020nov28", "processed", "2020nov28_sirius-b_pa.fits")),
]
fwhms = [
    7.985399074538929,
    7.644946278252633,
    8.216162087416189
]
# stims = [stimmap(res_cube, angle) for (res_cube, angle) in zip(res_cubes, angles)]
flat_res = [collapse(res_cube, angle) for (res_cube, angle) in zip(res_cubes, angles)]

parallax = SiriusB.parallax * 1e3 # mas
pxscale = SiriusB.pxscale * 1e3 # mas/px

errs = fwhms ./ 2

# guesses from GreeDS outputs of three epochs
pos = [
    110.21 ± errs[1]  93 ± errs[1];
    89.5 ± errs[2]  105 ± errs[2];
    92.5 ± errs[3]  111 ± errs[3];
] # px coordinates

ctr = hcat(center(flat_res[1])...)

delta_px = pos .- ctr
delta = delta_px .* pxscale

###
# plot the epochs
# first, zoom in
crop_res = @. crop(flat_res, 50)

tick_locs = range(3, 46, length=5)
tick_labs = @. string(round(pxscale / parallax * (tick_locs - 24.5), digits=1))

pycenter = hcat(center(crop_res[1])...) .- 1
pos_py = Measurements.value.(delta_px .+ pycenter)

py"""
import proplot as pro
import matplotlib.pyplot as plt
fig, axs = pro.subplots(ncols=3, wspace="1em", figwidth="7.5in")

axs[0].imshow($(crop_res[1]))
axs[0].add_artist(plt.Circle( $(pos_py[1, :]), $(errs[1]), lw=2, color="w", fill=False))
axs[0].format(title="Epoch 2020-02-04")

axs[1].imshow($(crop_res[2]))
axs[1].add_artist(plt.Circle( $(pos_py[2, :]), $(errs[2]), lw=2, color="w", fill=False))
axs[1].format(title="Epoch 2020-11-21")

axs[2].imshow($(crop_res[3]))
axs[2].add_artist(plt.Circle( $(pos_py[3, :]), $(errs[3]), lw=2, color="w", fill=False))
axs[2].format(title="Epoch 2020-11-28")
axs.format(
    xticks=$tick_locs,
    xticklabels=$tick_labs,
    yticks=$tick_locs,
    yticklabels=$tick_labs,
    xlabel="x [AU]",
    ylabel="y [AU]",
)

fig.save($(figuredir("orbit_frames.pdf")))
pro.close(fig)
"""

###

# translate into offset and angle

dist = mapslices(norm, delta, dims=2)
angs = mapslices(v -> atand(reverse(v)...) - 90, delta, dims=2)

epochs = [
    "2020feb04",
    "2020nov21",
    "2020nov28"
]
epochs_utc = [
    TTEpoch(2020, 2, 4),
    TTEpoch(2020, 11, 21),
    TTEpoch(2020, 11, 28)
]

epochs_mjd = @. AstroTime.value(modified_julian(epochs_utc))

data = DataFrame(
    epoch = epochs_mjd,
    object = fill(1, length(epochs)),
    sep = Measurements.value.(dist) |> vec,
    sep_err = Measurements.uncertainty.(dist) |> vec,
    pa = Measurements.value.(angs) |> vec,
    pa_err = Measurements.uncertainty.(angs) |> vec
)
@info data

fname = datadir("sirius-b_astrometry_data.csv")
data |> CSV.write(fname)

# now start setting up orbitize
mydriver = driver.Driver(
    fname,
    "OFTI",
    1, # number of bodies
    3.081, # total system mass M_sun https://iopscience.iop.org/article/10.3847/1538-4357/aa6af8/pdf
    parallax, # parallax in mas (EDR3)
    mass_err = 0.034,
    plx_err = 0.2313347
)

s = mydriver.sampler
orbits = s.run_sampler(Int(1e4))

# orbit plot
s.results.plot_orbits(start_mjd=epochs_mjd[1], num_orbits_to_plot=1000)
savefig(figuredir("orbit.pdf"))

# corner plot
s.results.plot_corner(param_list=["sma1", "ecc1", "inc1", "aop1", "pan1","tau1"])
savefig(figuredir("orbit_corner.pdf"))

max_stars_repo_name: mileslucas/sirius-b

max_issues_repo_name: mileslucas/sirius-b

max_forks_repo_name: mileslucas/sirius-b


text:
using Documenter
using Enigma

makedocs(
    # See https://github.com/JuliaDocs/Documenter.jl/issues/868
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    strict = true,
    sitename = "Enigma",
    pages = [
        "Home" => "index.md",
        "Tutorial" => "tutorial.md",
        "How-To" => "how_to.md",
        "P5js visualization" => "p5js.md",
        "Explanation" => "explanation.md",
        "Reference" => "reference.md",
    ]
)

deploydocs(
    repo = "github.com/Wikunia/Enigma.jl.git",
)


text:
function frule(
    (_, _, ΔA),
    ::typeof(select),
    op::Union{Function, SelectUnion},
    A::GBArray
)
    Ω = select(op, A)
    ∂Ω = mask(unthunk(ΔA), Ω, structural = true)
    return Ω, ∂Ω
end

function frule(
    (_, _, ΔA, _),
    ::typeof(select),
    op::Union{Function, SelectUnion},
    A::GBArray,
    thunk::Union{GBScalar, Nothing, valid_union}
)
    Ω = select(op, A, thunk)
    ∂Ω = mask(unthunk(ΔA), Ω, structural = true)
    return Ω, ∂Ω
end

function rrule(
    ::typeof(select),
    op::Union{Function, SelectUnion},
    A::GBArray
)
    out = select(op, A)
    function selectback(ΔΩ)
        ∂A = mask(unthunk(ΔΩ), out, structural = true)
        return NoTangent(), NoTangent(), ∂A
    end
    return out, selectback
end

function rrule(
    ::typeof(select),
    op::Union{Function, SelectUnion},
    A::GBArray,
    thunk::Union{GBScalar, Nothing, valid_union}
)
    out = select(op, A, thunk)
    function selectback(ΔΩ)
        ∂A = mask(unthunk(ΔΩ), out, structural = true)
        return NoTangent(), NoTangent(), ∂A, NoTangent()
    end
    return out, selectback
end


text:
### Constructors for setting up a Monte Carlo job

function MCJob{M<:MCModel, S<:MCSampler, T<:MCTuner}(m::M, s::S, r::SerialMC, t::T, j::Symbol)
  mcjob::MCJob
  if j == :plain
    mcjob = PlainMCJob(m, s, r, t)
  elseif j == :task
    mcjob = TaskMCJob(m, s, r, t)
  else
    error("Only :plain and :task jobs are available.")
  end
  mcjob
end

MCJob{M<:MCModel, S<:MCSampler, T<:MCTuner}(m::M, s::S, r::SerialMC; tuner::T=VanillaMCTuner(), job::Symbol=:task) =
  MCJob(m, s, r, tuner, job)

### Constructors for setting up a vector of Monte Carlo jobs

function MCJob{M<:MCModel, S<:MCSampler, T<:MCTuner}(m::Vector{M}, s::Vector{S}, r::Vector{SerialMC}, t::Vector{T},
  j::Vector{Symbol})
  @assert length(m) == length(s) == length(r) == length(t) == length(j)
    "Number of models, samplers, runners, tuners and job types not equal."
  map(MCJob, m, s, r, t, j)
end

MCJob{M<:MCModel, S<:MCSampler, T<:MCTuner}(m::Vector{M}, s::Vector{S}, r::Vector{SerialMC};
  tuners::Vector{T}=fill(VanillaMCTuner(), length(m)), jobs::Vector{Symbol}=fill(:task, length(m))) =
  map(MCJob, m, s, r, tuners, jobs)

### Functions for running jobs

function run(j::MCJob)
  if j.dim == 1
    if isa(j.runner[1], SerialMC)
      run(j.model[1], j.sampler[1], j.runner[1], j.tuner[1], j.jobtype)
    end
  end
end

run{J<:MCJob}(j::Vector{J}) = map(run, j)

### Functions for resuming jobs

function resume!(j::MCJob, c::MCChain; nsteps::Int=100)
  if j.dim == 1
    if isa(j.runner[1], SerialMC)
      resume!(j.model[1], j.sampler[1], j.runner[1], c, j.tuner[1], j.jobtype; nsteps=nsteps)
    end
  end
end

function resume(j::MCJob, c::MCChain; nsteps::Int=100)
  if j.dim == 1
    if isa(j.runner[1], SerialMC)
      resume!(deepcopy(j.model)[1], j.sampler[1], j.runner[1], c, j.tuner[1], j.jobtype; nsteps=nsteps)
    end
  end
end


text:
g(x) = x + 123456
function f(x)
    []
    Base.invokelatest(g, 0)
    Base.invokelatest(g, x)
    []
end
f(1.23)

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
module CGS

export C, GRAV, BOLTZMANN, STEFANBOLTZMANN, MHYDROGEN, YEAR, AU, PC, MSUN, PLANCK
export RSUN, MEARTH, REARTH, MJUPITER, RJUPITER, AJUPITER, KILOGRAM, METER, MMOON, RMOON

PLANCK = 6.62607004e-27
GRAV = 6.67384e-8
BOLTZMANN = 1.3806488e-16
STEFANBOLTZMANN = 5.670373e-5
MHYDROGEN = 1.00794*1.66053886e-24
YEAR = 3.15569e7
C = 29979245800e0

#astrophysical constants
AU = 1.49598e13
PC = 3.08568025e18

MSUN = 1.98892e33
RSUN = 6.955e10

MEARTH = 5.9742e27
REARTH = 6.3781e8
MMOON = 7.342e25
RMOON = 1737.1e5

MJUPITER = 1.8986e30
RJUPITER = 6.9911e9
AJUPITER = 4.84143144246472090 # AU

KILOGRAM = 1e3
METER = 1e2

end


text:
using RealisticLESbrary

import Interpolations
using Interpolations: Gridded, Linear, interpolate, extrapolate

function interpolate_surface_forcings(sose_surface_forcings, times; array_type=Array{Float64})
    Δt = times[2] - times[1]
    t_max = last(times)
    return (
         τx = InterpolatedProfile(array_type(sose_surface_forcings.τx),  times, Δt, t_max),
         τy = InterpolatedProfile(array_type(sose_surface_forcings.τy),  times, Δt, t_max),
         Qθ = InterpolatedProfile(array_type(sose_surface_forcings.Qθ),  times, Δt, t_max),
         Qs = InterpolatedProfile(array_type(sose_surface_forcings.Qs),  times, Δt, t_max),
        mld = InterpolatedProfile(array_type(sose_surface_forcings.mld), times, Δt, t_max)
    )
end

function interpolate_profile(profile, sose_grid, grid, times, array_type)
    Δt = times[2] - times[1]
    t_max = last(times)
    z_max = grid.zC[grid.Nz]

    ## First we interpolate from the SOSE stretched grid onto the Oceananigans.jl regular grid.

    # Coordinates needs to be in increasing order for Interpolations.jl.
    sose_zC = reverse(sose_grid.zC)
    profile = reverse(profile, dims=2)

    ℑprofile = interpolate((times, sose_zC), profile, Gridded(Linear()))

    # We might go out of bounds if the Oceananigans is fine. In that case, grab the nearest value.
    ℑprofile = extrapolate(ℑprofile, Interpolations.Flat())

    interior_zC = grid.zC[1:grid.Nz]
    interpolated_data = ℑprofile.(times', interior_zC)

    ## Then we construct and return an InterpolatedProfileTimeSeries for fast linear interpolation in kernels.

    return InterpolatedProfileTimeSeries(array_type(interpolated_data), interior_zC, times, grid.Δz, Δt, z_max, t_max)

    return ℑprofile
end

function interpolate_profiles(sose_profiles, sose_grid, grid, times; array_type=Array{Float64})
    return (
        U    = interpolate_profile(sose_profiles.U,    sose_grid, grid, times, array_type),
        V    = interpolate_profile(sose_profiles.V,    sose_grid, grid, times, array_type),
        Θ    = interpolate_profile(sose_profiles.Θ,    sose_grid, grid, times, array_type),
        S    = interpolate_profile(sose_profiles.S,    sose_grid, grid, times, array_type),
        Ugeo = interpolate_profile(sose_profiles.Ugeo, sose_grid, grid, times, array_type),
        Vgeo = interpolate_profile(sose_profiles.Vgeo, sose_grid, grid, times, array_type)
    )
end


text:
using Nextion
using Test

@testset "picture" begin

    nexSerial = NexSerial("/dev/ttyUSB0")

    nexPicture = NexPicture(nexSerial, Name("p0"), cid=ComponentID(1))

    send(nexSerial, "page pg_pic")

    for i in 1:4
        sleep(1)
        #nexPicture.picture = i
        nexPicture.picture = Picture(i)
    end

    close(nexSerial)
    @test true

end


text:
import Base

include("production.jl")

"""
A grammar, represented as a tuple ``G=(N,T,P,S)``
"""
struct Grammar
    "The nonterminal symbols Set"
    N::Set
    "The terminal symbols Set"
    T::Set
    "The productions Array"
    P::Array
    "The starting symbol"
    S::AbstractString
    iscontextfree::Bool
end

### Constructors ###

"""
    Grammar(N::Set, T::Set, P::Array, S::AbstractString)::Grammar
Builds a [`Grammar`](@ref) obtained from the given string of productions.
"""
function Grammar(N, T, P, S)
    cf = all(x -> isa(x.left, AbstractString), P)
    if (N ∩ T) ≠ Set()
        throw(ArgumentError("The set of terminals and nonterminals are not disjoint, but have " * string(collect(N ∩ T)) * " in common"))
    end
    if S ∉ N
        throw(ArgumentError("The start symbol is not a nonterminal."))
    end
    if cf
        badprods = [p for p ∈ P if p.left ∉ N]
        if !isempty(badprods)
            throw(ArgumentError("The following productions have a left-hand side that is not a nonterminal: " * string(badprods)))
        end
    end
    badprods = [p for p ∈ P if (Set(astype0(p).left) ∪ Set(p.right)) ⊈ (N ∪ T ∪ Set(["ε"]))]
    if ~isempty(badprods)
        throw(ArgumentError("The following productions contain symbols that are neither terminals or nonterminals: " * string(badprods)))
    end
    Grammar(N, T, P, S, cf)
end

"""
    Grammar(prods::AbstractString, iscontextfree = true)::Grammar
Builds a [`Grammar`](@ref) obtained from the given string of productions.
"""
function Grammar(prods::AbstractString, iscontextfree = true)::Grammar
    P = parseproduction(prods, iscontextfree)
    S = nothing
    N = nothing
    T = nothing
    if iscontextfree
        S = P[1].left
        N = Set(map(x -> x.left, P))
        T = Set(vcat(map(x -> x.right, P)...)) - N - "ε"
    else
        S = P[1].left[1]
        symbols = Set(vcat(map(x -> x.left, P)...)) ∪ Set(vcat(map(x -> x.right, P)...))
        N = Set(filter(x -> isuppercase(x[1]), symbols))
        T = symbols - N - "ε"
    end
    G = Grammar(N, T, P, S)
    if iscontextfree
        if ~G.iscontextfree
            throw(ArgumentError("The resulting grammar is not context-free, even if so requested."))
        end
    end
    return G
end

### Functions ###

"""
    alternatives(g::Grammar, N::Array)::Array
Returns all the right-hand sides alternatives matching the given nonterminal.
"""
alternatives(g::Grammar, N::Union{AbstractString, Iterable})::Array = [P.right for P ∈ g.P if P.left == (isa(N,Iterable) ? collect(N) : N)]

"""
    restrict(g::Grammar, symbols::Set)
Returns a [`Grammar`](@ref) using only the given symbols.
"""
restrict(g::Grammar, symbols::Set) = if g.S ∉ symbols throw(ArgumentError("The start symbol must be present among the symbols to keep.")) else Grammar(g.N ∩ symbols, g.T ∩ symbols, [P for P ∈ g.P if (Set([P.left]) ∪ Set(P.right)) ≤ symbols], g.S) end

### Operators ###

Base.:(==)(x::Grammar, y::Grammar) = (x.N, x.T, sort(x.P), x.S) == (y.N, y.T, sort(y.P), y.S)

Base.hash(g::Grammar) = Base.hash((g.N, g.T, sort(g.P), g.S))

Base.show(io::IO, g::Grammar) = Base.show(io, string("Grammar(N=", g.N, ", T=", g.T, ", P=", g.P, "S=", g.S, ")"))


text:
using Plots; plotlyjs()
using FFTW

# Construct a signal
f = 1
fs = 5
ts = 1 / fs
l = 2^10 + 1
t = (-(l - 1) ÷ 2 : (l - 1) ÷ 2) * ts
x = sinc.(f * t)

# Construct the frequency spectrum
xf = fft(x)
ff = (0 : 1 / (l - 1) : 1) * fs

# Signal energy
Δt = ts
Δf = fs / (l - 1)
Et = sum(abs.(x).^2) * Δt
Ef = sum(abs.(xf).^2) * Δf / (fs^2)
@show Et, Ef

# Plots
plt = plot(layout=2)
plot!(t, x, subplot=1)
plot!(ff, abs.(xf), subplot=2)
display(plt)


text:
# Spatial calculations

export
    vector1D,
    vector,
    wrapcoords,
    wrapcoordsvec

"""
    vector1D(c1, c2, side_length)

Displacement between two 1D coordinate values from c1 to c2, accounting for
the bounding box.
The minimum image convention is used, so the displacement is to the closest
version of the coordinate accounting for the periodic boundaries.
"""
function vector1D(c1, c2, side_length)
    if c1 < c2
        return (c2 - c1) < (c1 - c2 + side_length) ? (c2 - c1) : (c2 - c1 - side_length)
    else
        return (c1 - c2) < (c2 - c1 + side_length) ? (c2 - c1) : (c2 - c1 + side_length)
    end
end

"""
    vector(c1, c2, box_size)

Displacement between two coordinate values, accounting for the bounding box.
The minimum image convention is used, so the displacement is to the closest
version of the coordinates accounting for the periodic boundaries.
"""
vector(c1, c2, box_size) = vector1D.(c1, c2, box_size)

@generated function vector(c1::SVector{N}, c2::SVector{N}, box_size) where N
    quote
        Base.Cartesian.@ncall $N SVector{$N} i->vector1D(c1[i], c2[i], box_size[i])
    end
end

sqdistance(i, j, coords, box_size) = sum(abs2, vector(coords[i], coords[j], box_size))

"""
    wrapcoords(c, side_length)

Ensure a 1D coordinate is within the simulation box and return the coordinate.
"""
wrapcoords(c, side_length) = c - floor(c / side_length) * side_length

"""
    wrapcoordsvec(c, box_size)

Ensure a coordinate is within the simulation box and return the coordinate.
"""
wrapcoordsvec(v, box_size) = wrapcoords.(v, box_size)


text:
@userplot SurfPlot

RecipesBase.@recipe function f(p::SurfPlot)
    if length(p.args) != 1 || !(typeof(p.args[1]) <: PSAStruct)
        error("surfplot must be given a PSAStruct. Got $(typeof(p.args))")
    end
    ps_data = p.args[1]
    zoom = ps_data.zoom_list[ps_data.zoom_pos]
    if !zoom.computed
        return
    end

    @series begin
        seriestype := :surface
        zoom.x, zoom.y, -log10.(zoom.Z)
    end
end


text:
include("FTypesModule.jl")
include("BoxModule.jl")
include("MatrixUtilityModule.jl")
include("PhysicalUnitModule.jl")
include("RotationUtilModule.jl")
include("CSysModule.jl")
include("FESetModule.jl")
include("FENodeSetModule.jl")
include("FENodeToFEMapModule.jl")
include("FieldModule.jl")
include("GeneralFieldModule.jl")
include("NodalFieldModule.jl")
include("ElementalFieldModule.jl")
include("MeshUtilModule.jl")
include("MeshSelectionModule.jl")
include("MeshExportModule.jl")
include("MeshModificationModule.jl")
include("MeshImportModule.jl")
include("MeshQuadrilateralModule.jl")
include("MeshLineModule.jl")
include("MeshTriangleModule.jl")
include("MeshHexahedronModule.jl")
include("MeshTetrahedronModule.jl")
include("VectorCacheModule.jl")
include("ForceIntensityModule.jl")
include("SurfaceNormalModule.jl")
include("MatHeatDiffModule.jl")
include("MatAcoustFluidModule.jl")
include("DeforModelRedModule.jl")
include("MatDeforModule.jl")
include("MatDeforElastIsoModule.jl")
include("MatDeforElastOrthoModule.jl")
include("AssemblyModule.jl")
include("IntegRuleModule.jl")
include("IntegDataModule.jl")
include("FEMMBaseModule.jl")
include("FEMMHeatDiffModule.jl")
include("FEMMHeatDiffSurfModule.jl")
include("FEMMAcoustModule.jl")
include("FEMMAcoustSurfModule.jl")
include("FEMMDeforLinearBaseModule.jl")
include("FEMMDeforLinearModule.jl")
include("FEMMDeforWinklerModule.jl")
include("FEMMDeforSurfaceDampingModule.jl")
include("FEMMDeforLinearMSModule.jl")
include("VoxelBoxModule.jl")
include("TetRemeshingModule.jl")
include("VoxelTetMeshingModule.jl")
include("AlgoBaseModule.jl")
include("AlgoAcoustModule.jl")
include("AlgoHeatDiffModule.jl")
include("AlgoDeforLinearModule.jl")


text:
# AUTO GENERATED FILE - DO NOT EDIT

export sidemenu

"""
    sidemenu(;kwargs...)

A SideMenu component.
ExampleComponent is an example component.
It takes a property, `label`, and
displays it.
It renders an input with the property `value`
which is editable by the user.
Keyword arguments:
- `id` (String; optional): The ID used to identify this component in Dash callbacks.
- `items` (Array; required): An array of id, label
- `className` (String; optional): Navbar class name
"""
function sidemenu(; kwargs...)
        available_props = Symbol[:id, :items, :className]
        wild_props = Symbol[]
        return Component("sidemenu", "SideMenu", "sm", available_props, wild_props; kwargs...)
end

max_stars_repo_name: lrahal/dash-component-boilerplate

max_issues_repo_name: lrahal/dash-component-boilerplate

max_forks_repo_name: lrahal/dash-component-boilerplate


text:
function useplots()
  res1 = ""
  try
  	res = Pkg.installed("Plots")
  catch
    println("Package Plots not installed, please run Pkg.add(\"Plots\")")
    return
  end

  res2 = ""
  try
  	res = Pkg.installed("GR")
  catch
    println("Package GR.jl not installed, please run Pkg.add(\"GR\")")
    return
  end

  if typeof(res1) == VersionNumber && typeof(res2) == VersionNumber
    eval(quote
      using Plots
      using GR
    end)
  end
  typeof(res1) == VersionNumber && typeof(res2) == VersionNumber
end


text:
using TopOpt, Nonconvex, Zygote, FiniteDifferences, LinearAlgebra, Test, Random, SparseArrays
const FDM = FiniteDifferences

Random.seed!(1)

@testset "Compliance" begin
    nels = (10, 10)
    problem = HalfMBB(Val{:Linear}, nels, (1.0, 1.0), 1.0, 0.3, 1.0)
    for p in (1.0, 2.0, 3.0)
        solver = FEASolver(Displacement, Direct, problem, xmin = 0.01, penalty = TopOpt.PowerPenalty(p))
        comp = Compliance(problem, solver)
        for i in 1:3
            x = clamp.(rand(prod(nels)), 0.1, 1.0)
            val1, grad1 = Nonconvex.value_gradient(comp, x)
            val2, grad2 = comp(x), Zygote.gradient(comp, x)[1]
            grad3 = FDM.grad(central_fdm(5, 1), comp, x)[1]
            @test val1 == val2
            @test norm(grad1 - grad2) == 0
            @test norm(grad2 - grad3) <= 1e-5
        end
    end
end

@testset "Volume" begin
    nels = (10, 10)
    problem = HalfMBB(Val{:Linear}, nels, (1.0, 1.0), 1.0, 0.3, 1.0)
    for p in (1.0, 2.0, 3.0)
        solver = FEASolver(Displacement, Direct, problem, xmin = 0.01, penalty = TopOpt.PowerPenalty(p))
        vol = Volume(problem, solver)
        constr = x -> vol(x) - 0.3
        for i in 1:3
            x = rand(prod(nels))
            val1, grad1 = Nonconvex.value_gradient(constr, x)
            val2, grad2 = constr(x), Zygote.gradient(constr, x)[1]
            grad3 = FDM.grad(central_fdm(5, 1), constr, x)[1]
            @test val1 == val2
            @test norm(grad1 - grad2) == 0
            @test norm(grad2 - grad3) <= 1e-5
        end
    end
end

@testset "DensityFilter" begin
    nels = (10, 10)
    problem = HalfMBB(Val{:Linear}, nels, (1.0, 1.0), 1.0, 0.3, 1.0)
    for p in (1.0, 2.0, 3.0)
        solver = FEASolver(Displacement, Direct, problem, xmin = 0.01, penalty = TopOpt.PowerPenalty(p))
        filter = TopOpt.DensityFilter(solver, rmin = 4.0)
        for i in 1:3
            x = rand(prod(nels))
            v = rand(prod(nels))
            f = FunctionWrapper(x -> dot(filter(x), v), 1)
            val1, grad1 = Nonconvex.value_gradient(f, x)
            val2, grad2 = f(x), Zygote.gradient(f, x)[1]
            grad3 = FDM.grad(central_fdm(5, 1), f, x)[1]
            @test val1 == val2
            @test norm(grad1 - grad2) == 0
            @test norm(grad2 - grad3) <= 1e-5
        end
    end
end

@testset "SensFilter" begin
    nels = (10, 10)
    problem = PointLoadCantilever(Val{:Linear}, nels, (1.0, 1.0))
    solver = FEASolver(
        Displacement, Direct, problem, xmin = 1e-3, penalty = PowerPenalty(3.0),
    )
    sensfilter = SensFilter(solver, rmin = 4.0)
    x = rand(length(solver.vars))
    y = rand(length(x))
    @test sensfilter(x) == x
    grad = Zygote.gradient(x -> dot(sensfilter(x), y), x)[1]
    @test grad != y
end

@testset "Block compliance" begin
    nels = (10, 10)
    nloads = 10
    base_problem = PointLoadCantilever(Val{:Linear}, nels, (1.0, 1.0), 1.0, 0.3, 1.0)
    dense_load_inds = vec(TopOpt.TopOptProblems.get_surface_dofs(base_problem))
    dense_rank = 3
    F = spzeros(TopOpt.Ferrite.ndofs(base_problem.ch.dh), nloads)
    Fsize = size(F)
    for i in 1:dense_rank
        F += sparsevec(dense_load_inds, randn(length(dense_load_inds)) / dense_rank, Fsize[1]) * randn(Fsize[2])'
    end
    problem = MultiLoad(base_problem, F)
    for p in (1.0, 2.0, 3.0)
        solver = FEASolver(Displacement, Direct, problem, xmin = 0.01, penalty = TopOpt.PowerPenalty(p))
        exact_svd_block = BlockCompliance(problem, solver, method = :exact)
        constr = Nonconvex.FunctionWrapper(x -> exact_svd_block(x) .- 1000.0, length(exact_svd_block(solver.vars)))
        for i in 1:3
            x = clamp.(rand(prod(nels)), 0.1, 1.0)
            v = rand(nloads)
            f = FunctionWrapper(x -> dot(constr(x), v), 1)
            val1, grad1 = Nonconvex.value_gradient(f, x)
            val2, grad2 = f(x), Zygote.gradient(f, x)[1]
            grad3 = FDM.grad(central_fdm(5, 1), f, x)[1]
            @test val1 == val2
            @test norm(grad1 - grad2) == 0
            @test norm(grad2 - grad3) <= 1e-4
        end
    end
end

@testset "Local stress" begin
    nels = (10, 10)
    problem = PointLoadCantilever(Val{:Quadratic}, nels, (1.0, 1.0), 1.0, 0.3, 1.0)
    for F in (MacroVonMisesStress, MicroVonMisesStress)
        for p in (1.0, 2.0, 3.0)
            solver = FEASolver(Displacement, Direct, problem, xmin = 0.001, penalty = TopOpt.PowerPenalty(p))
            stress = F(solver)
            for i in 1:3
                x = clamp.(rand(prod(nels)), 0.1, 1.0)
                v = rand(prod(nels))
                f = FunctionWrapper(x -> dot(stress(x), v), 1)
                val1, grad1 = Nonconvex.value_gradient(f, x)
                val2, grad2 = f(x), Zygote.gradient(f, x)[1]
                grad3 = FDM.grad(central_fdm(5, 1), f, x)[1]
                @test val1 == val2
                @test norm(grad1 - grad2) == 0
                @test norm(grad1 - grad3) <= 1e-5
            end
        end
    end

    nels = (10, 10, 2)
    problem = PointLoadCantilever(Val{:Quadratic}, nels, (1.0, 1.0, 1.0), 1.0, 0.3, 1.0)
    for F in (MacroVonMisesStress, MicroVonMisesStress)
        for p in (1.0, 2.0, 3.0)
            solver = FEASolver(Displacement, Direct, problem, xmin = 0.001, penalty = TopOpt.PowerPenalty(p))
            stress = F(solver)
            for i in 1:3
                x = clamp.(rand(prod(nels)), 0.1, 1.0)
                v = rand(prod(nels))
                f = FunctionWrapper(x -> dot(stress(x), v), 1)
                val1, grad1 = Nonconvex.value_gradient(f, x)
                val2, grad2 = f(x), Zygote.gradient(f, x)[1]
                grad3 = FDM.grad(central_fdm(5, 1), f, x)[1]
                @test val1 == val2
                @test norm(grad1 - grad2) == 0
                @test norm(grad1 - grad3) <= 1e-5
            end
        end
    end
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

#Various sizes of hexagonal clusters that are useful in lenslet design. To see what they look like use the functions in HexTilings to make lattices with color properties and then draw them with Vis.draw
#example: Vis.draw(hex3RGB())

function hex3()
    clusterelements = SVector((0, 0), (-1, 0), (-1, 1))
    eltlattice = Repeat.HexBasis1()
    clusterbasis = Repeat.LatticeBasis((-1, 2), (2, -1))
    return Repeat.LatticeCluster(clusterbasis, eltlattice, clusterelements)

end
export hex3

function hex4()
    clusterelements = SVector((-1, 0), (-1, 1), (0, 0), (0, 1))
    eltlattice = Repeat.HexBasis1()
    clusterbasis = Repeat.LatticeBasis((2, 0), (0, 2))
    return Repeat.LatticeCluster(clusterbasis, eltlattice, clusterelements)
end
export hex4

function hex7()
    clusterelements = SVector((1, -1), (0, -1), (-1, 0), (0, 0), (-1, 1), (0, 1), (1, 0), )
    eltlattice = Repeat.HexBasis1()
    clusterbasis = Repeat.LatticeBasis((3, 0), (0, 3))
    clusterbasis = Repeat.LatticeBasis((2, 2), (-3, 2))
    return Repeat.LatticeCluster(clusterbasis, eltlattice, clusterelements)
end
export hex7

function hex9()
    clusterelements = SVector((-1, -1), (-1, 0), (-2, 1), (0, -1), (0, 0), (-1, 1), (1, -1), (1, 0), (0, 1))
    eltlattice = Repeat.HexBasis1()
    clusterbasis = Repeat.LatticeBasis((3, 0), (0, 3))
    return Repeat.LatticeCluster(clusterbasis, eltlattice, clusterelements)
end
export hex9

function hex12()
    clusterelements = SVector(
        (-1, 1),(0, 1),
        (-1, 0),(0, 0),(1, 0),
        (-1, -1),(0, -1),(1, -1),(2, -1),
        (0, -2),(1, -2),(2, -2)
    )

    eltlattice = Repeat.HexBasis3()
    clusterbasis = Repeat.LatticeBasis((2, 2), (-3, 2))
    return Repeat.LatticeCluster(clusterbasis, eltlattice, clusterelements)
end
export hex12

function hex19()
    clusterelements = SVector(
        (0, 0), # ring 0
        (-1, 0),(0, -1),(1, -1),(1, 0),(0, 1),(-1, 1), # ring1
        (-2, 0),(-1, -1),(0, -2),(1, -2),(2, -2),(2, -1),(2, 0),(1, 1),(0, 2),(-1, 2),(-2, 2),(-2, 1) # ring2
        )

    eltlattice = Repeat.HexBasis3()
    clusterbasis = Repeat.LatticeBasis((5, 0), (-2, 4))
    return Repeat.LatticeCluster(clusterbasis, eltlattice, clusterelements)
end
export hex19


text:
"""
    $(SIGNATURES)

Convenience function to load yaml config files, commonly used to prepare the required beamforming computational objects.
"""
function loadConfigFile(file::A) where {A <: AbstractString}
  YAML.load(open(file))
end

"""
    $(SIGNATURES)

Prepare a sas2d factor to use in the factor graph where `totalPhones` is the size of the SAS array.  Pass a known `cfgd::Dict{String,} for faster load times.`
"""
function prepareSAS2DFactor(totalPhones::Int, csvWaveData::Array{<:Real};
              cfgd::Dict=loadConfigFile(joinpath(dirname(pathof(Caesar)),"config","SAS2D.yaml")), rangemodel=:Rayleigh, chirpFile=joinpath(dirname(pathof(Caesar)),"config","chirp250.txt"))
  # # @assert size(csvWaveData,2) == totalPhones
  #Constant Acoustic Params
  fFloor = cfgd["freq_lower_cutoff"];
  fCeil = cfgd["freq_upper_cutoff"];
  fSampling = cfgd["freq_sampling"];
  soundSpeed = cfgd["sound_speed"];
  nFFT_full = cfgd["nfft_full"];
  nFFT_czt = cfgd["nfft_czt"];
  azimuthDivs = cfgd["azimuth_divs"];
  rngSNRfloor = cfgd["range_snr_floor"];
  aziSNRfloor = cfgd["azi_snr_floor"];

  #Params subject to change/ Heuristics
  nPhones = totalPhones - 1;
  azimuths = LinRange(0,360,azimuthDivs)*pi/180;
  #Data preprocessing MF + CZT
  w = exp(-2im*pi*(fCeil-fFloor)/(nFFT_czt*fSampling))
  a = exp(2im*pi*fFloor/fSampling)

  chirpIn = readdlm(chirpFile,',',Float64,'\n')

  #Matched Filter on Data In
  nFFT_full = nextpow(2,size(csvWaveData,1))  # MF
  mfData = zeros(Complex{Float64}, nFFT_full, totalPhones)
  mf = prepMF(chirpIn,nFFT_full,totalPhones) # MF
  mf(csvWaveData,mfData) # MF

  cztData = zeros(Complex{Float64}, nFFT_czt,totalPhones)
  filterCZT = prepCZTFilter(nFFT_full,totalPhones,w,nFFT_czt,a)
  filterCZT(mfData,cztData)

  #CBF Filter shared by

  FFTfreqs = collect(LinRange(fFloor,fCeil,nFFT_czt))
  cfgCBF_init = CBFFilterConfig(fFloor,fCeil,nFFT_czt,totalPhones,azimuths,soundSpeed, FFTfreqs)
  cfgCBF_init_LIE = CBFFilterConfig(fFloor,fCeil,nFFT_czt,nPhones,Float64[0.0;],soundSpeed, FFTfreqs)

  sas2d = SASBearing2D()
  if rangemodel == :Rayleigh
    sas2d.rangemodel = Rayleigh(100.0)
  elseif rangemodel == :Correlator
      range_mfData = zeros(Complex{Float64}, nFFT_full, totalPhones)
      range_mf = prepMF(chirpIn,nFFT_full,totalPhones) # MF
      range_mf(csvWaveData,range_mfData) # matched filter
      ranget = [0:8000-500;]*soundSpeed/fSampling
      sas2d.rangemodel = Vector{IIF.AliasingScalarSampler}(undef, totalPhones)
      for i in 1:totalPhones
        sas2d.rangemodel[i] = IIF.AliasingScalarSampler(ranget, exp.(norm.(range_mfData[1:7501,i])), SNRfloor=rngSNRfloor)
      end
  else
      error("No can do: what is a $rangemodel rangemodel?")
  end

  #thread safe
  sas2d.cfgTotal = cfgCBF_init
  sas2d.cfgLIE = cfgCBF_init_LIE
  sas2d.cfg = deepcopy(cfgd)
  sas2d.waveformsRaw = deepcopy(csvWaveData)
  sas2d.waveformsIn = cztData
  sas2d.debugging = false

  #create memory
  prepareSASMemory(sas2d)

  return sas2d
end

function prepareSASMemory(sas2d::SASBearing2D)::Nothing
  #Constant Acoustic Params
  fFloor = sas2d.cfg["freq_lower_cutoff"];
  fCeil = sas2d.cfg["freq_upper_cutoff"];
  fSampling = sas2d.cfg["freq_sampling"];
  soundSpeed = sas2d.cfg["sound_speed"];
  nFFT_full = sas2d.cfg["nfft_full"];
  nFFT_czt = sas2d.cfg["nfft_czt"];
  azimuthDivs = sas2d.cfg["azimuth_divs"];
  rngSNRfloor = sas2d.cfg["range_snr_floor"];
  aziSNRfloor = sas2d.cfg["azi_snr_floor"];
  totalPhones = size(sas2d.waveformsRaw,2)
  azimuths = LinRange(0,360,azimuthDivs)*pi/180;

  #create memory
  sas2d.threadreuse = Vector{SASREUSE}(undef, Threads.nthreads())
  for thritr in 1:Threads.nthreads()
      sas2d.threadreuse[thritr] = SASREUSE()
      sas2d.threadreuse[thritr].CBFFull = zeros(Complex{Float64}, getCBFFilter2Dsize(sas2d.cfgTotal));
      sas2d.threadreuse[thritr].CBFLIE = zeros(Complex{Float64}, getCBFFilter2Dsize(sas2d.cfgLIE));
      sas2d.threadreuse[thritr].BFOutFull = zeros(length(azimuths))
      sas2d.threadreuse[thritr].BFOutLIE = zeros(Complex{Float64},nFFT_czt)
      sas2d.threadreuse[thritr].BFOutLIETemp = zeros(Complex{Float64},nFFT_czt)
      sas2d.threadreuse[thritr].phaseshiftLOO = zeros(Complex{Float64},nFFT_czt)
      sas2d.threadreuse[thritr].arrayPos = zeros(Float64,totalPhones,2)
      sas2d.threadreuse[thritr].arrayPosLIE = zeros(Float64,totalPhones-1,2)
      sas2d.threadreuse[thritr].waveformsLIE = zeros(Complex{Float64},nFFT_czt,totalPhones-1)
      sas2d.threadreuse[thritr].waveformsLOOc = zeros(Complex{Float64},nFFT_czt)

      sas2d.threadreuse[thritr].BFtemp = zeros(Complex{Float64}, totalPhones,nFFT_czt)
      sas2d.threadreuse[thritr].LIEtemp = zeros(Complex{Float64}, totalPhones-1,nFFT_czt)

      # caching hack
      sas2d.threadreuse[thritr].hackazi = (Int[0;], Float64[0.0;])
      sas2d.threadreuse[thritr].oncebackidx = Int[0;]

      # a few more memory reuse positions
      sas2d.threadreuse[thritr].sourceXY = zeros(2)

      # beamformer update
      sas2d.threadreuse[thritr].temp1 = zeros(Complex{Float64},nFFT_czt);
      sas2d.threadreuse[thritr].temp2 = zeros(Complex{Float64},size(sas2d.waveformsIn));
      sas2d.threadreuse[thritr].temp = zeros(Complex{Float64},nFFT_czt,totalPhones-1);

      sas2d.threadreuse[thritr].dbg = SASDebug()
      reset!(sas2d.threadreuse[thritr].dbg)
  end
end

# plot(x=azimuths, y=acceptRatio, Geom.path())
# careful -- this great function depends on how you discretize the azimuths in CBF
function sampleDiscrete!(smpls::Vector{T}, BFresult::Vector{T}) where {T <: Union{Float64, Float32}}
  N = length(smpls)
  y = 1/2pi
  acceptRatio = 2pi*BFresult

  sample = zeros(1)
  sampleInt = zeros(Int, 1)

  d2r = 2*pi/180.0
  i = 0
  while true
     sample[1] = 179*rand()
     sampleInt[1] = ceil.(Int, sample[1])
     # acceptRatio[sampleInt[1]]
      if acceptRatio[sampleInt[1]] > rand()
         i += 1
         smpls[i] = sample[1] * d2r
        if i >= N
          break;
        end
      end
  end
  nothing
end

# this is a direct copy from TransformUtils
function wrapRad(th::Float64)
  if th >= pi
    th -= 2.0*pi
  end
  if th < -pi
    th += 2.0*pi
  end
  return th
end