text:
struct LoadedDataset
    name::String
    parent::Group
    header_offset#::RelOffset
    header_size::Int
    data_offset#::RelOffset
    data_length::Int
    dataspace
    datatype_class
    datatype_offset
    filter_id
    attrs
    layout_class # Compact / Contiguous / Chunked
    #storage_message
    juliatype
    odr

end

read_dataset(f::JLDFile, name::AbstractString) =
    read_dataset(f.root_group, name)

function read_dataset(g::Group, name::AbstractString)
    f = g.f
    f.n_times_opened == 0 && throw(ArgumentError("file is closed"))

    (g, name) = pathize(g, name, false)

    roffset = lookup_offset(g, name)
    if roffset == UNDEFINED_ADDRESS
        haskey(g.unwritten_child_groups, name) && return g.unwritten_child_groups[name]
        throw(KeyError(name))
    end

    if isgroup(f, roffset)
        let loaded_groups = f.loaded_groups
            get!(()->load_group(f, roffset), loaded_groups, roffset)
        end
    else
        read_dataset(g, roffset, name)
    end
end

function read_dataset(group::Group, offset::RelOffset, name)
    f = group.f
    io = f.io
    seek(io, fileoffset(f, offset))
    cio = begin_checksum_read(io)
    sz = read_obj_start(cio)
    pmax = position(cio) + sz

    # Messages
    dataspace = ReadDataspace()
    attrs = EMPTY_READ_ATTRIBUTES
    datatype_class::UInt8 = 0
    datatype_offset::Int64 = 0
    data_offset::Int64 = 0
    data_length::Int = -1
    chunked_storage::Bool = false
    filter_id::UInt16 = 0

    storage_type = 0xFF
    while position(cio) <= pmax-4
        msg = jlread(cio, HeaderMessage)
        endpos = position(cio) + msg.size
        if msg.msg_type == HM_DATASPACE
            dataspace = read_dataspace_message(cio)
        elseif msg.msg_type == HM_DATATYPE
            datatype_class, datatype_offset = read_datatype_message(cio, f, (msg.flags & 2) == 2)
        elseif msg.msg_type == HM_FILL_VALUE
            (jlread(cio, UInt8) == 3 && jlread(cio, UInt8) == 0x09) || throw(UnsupportedFeatureException())
        elseif msg.msg_type == HM_DATA_LAYOUT
            jlread(cio, UInt8) == 4 || throw(UnsupportedVersionException())
            storage_type = jlread(cio, UInt8)
            if storage_type == LC_COMPACT_STORAGE
                data_length = jlread(cio, UInt16)
                data_offset = position(cio)
            elseif storage_type == LC_CONTIGUOUS_STORAGE
                data_offset = fileoffset(f, jlread(cio, RelOffset))
                data_length = jlread(cio, Length)
            elseif storage_type == LC_CHUNKED_STORAGE
                # TODO: validate this
                flags = jlread(cio, UInt8)
                dimensionality = jlread(cio, UInt8)
                dimensionality_size = jlread(cio, UInt8)
                skip(cio, Int(dimensionality)*Int(dimensionality_size))

                chunk_indexing_type = jlread(cio, UInt8)
                chunk_indexing_type == 1 || throw(UnsupportedFeatureException("Unknown chunk indexing type"))
                data_length = jlread(cio, Length)
                jlread(cio, UInt32)
                data_offset = fileoffset(f, jlread(cio, RelOffset))
                chunked_storage = true
            else
                throw(UnsupportedFeatureException("Unknown data layout"))
            end
        elseif msg.msg_type == HM_FILTER_PIPELINE
            version = jlread(cio, UInt8)
            version == 2 || throw(UnsupportedVersionException("Filter Pipeline Message version $version is not implemented"))
            nfilters = jlread(cio, UInt8)
            nfilters == 1 || throw(UnsupportedFeatureException())
            filter_id = jlread(cio, UInt16)
            issupported_filter(filter_id) || throw(UnsupportedFeatureException("Unknown Compression Filter $filter_id"))
        elseif msg.msg_type == HM_ATTRIBUTE
            if attrs === EMPTY_READ_ATTRIBUTES
                attrs = ReadAttribute[read_attribute(cio, f)]
            else
                push!(attrs, read_attribute(cio, f))
            end
        elseif (msg.flags & 2^3) != 0
            throw(UnsupportedFeatureException())
        end
        seek(cio, endpos)
    end
    seek(cio, pmax)

    filter_id != 0 && !chunked_storage && throw(InvalidDataException("Compressed data must be chunked"))

    # Checksum
    end_checksum(cio) == jlread(io, UInt32) || throw(InvalidDataException("Invalid Checksum"))

    headersize = jlsizeof(ObjectStart) + size_size(sz) + sz + 4
    LoadedDataset(
        name, group, offset, headersize, data_offset, data_length,
        dataspace, datatype_class, datatype_offset, filter_id, attrs,
         storage_type,nothing, nothing)

end

function load_ds(ds::LoadedDataset)
    ds.parent.f.n_times_opened == 0 && throw(ArgumentError("file is closed"))

    read_data(ds.parent.f,
        ds.dataspace,
        ds.datatype_class,
        ds.datatype_offset,
        ds.data_offset,
        ds.data_length,
        ds.filter_id,
        ds.header_offset,
        ds.attrs)
end

function edit!(ds::LoadedDataset, data, wsession=JLDWriteSession())
    @assert ds.parent.f.writable "Cannot edit in read-only mode"

    odr = objodr(data)
    f = ds.parent.f
    io = f.io
    datatype_class = ds.datatype_class
    datatype_offset = ds.datatype_offset
    data_offset = ds.data_offset
    if datatype_class == typemax(UInt8) # Committed datatype
        rr = jltype(f, f.datatype_locations[h5offset(f, datatype_offset)])
    else
        seek(io, datatype_offset)
        @read_datatype io datatype_class dt begin
            rr = jltype(f, dt)
        end
    end

    wodr = typeof(rr).parameters[2]
    @assert wodr == odr "Types do not match $wodr (written) $odr (new)"
    @assert ds.dataspace.dataspace_type == DS_SCALAR
    seek(io, data_offset)
    write_data(io, f, data, odr, datamode(odr), wsession)

    # Update checksum for compact storage
    if ds.layout_class == LC_COMPACT_STORAGE
        offset = fileoffset(f, ds.header_offset)
        update_checksum(io, offset, offset+ds.header_size)
    end
    return data
end

function JLD2.update_checksum(io, start_offset::Int, checksum_offset)
    seek(io, start_offset)
    cio = begin_checksum_read(io)
    seek(cio, checksum_offset-4)
    seek(io, checksum_offset-4)
    jlwrite(io, end_checksum(cio))
end


text:
using UnicodeData
using Test

@testset "Testing Proper Translation" begin
    @test udata('A') == "Latin Capital Letter A"
    @test udata('A') == udata(0x0041)

    @test udata('\uFF01') == "Fullwidth Exclamation Mark"
    @test udata('\uFF00') == nothing

    # Test the higher planes as well.
    @test udata(0x13201) == "Egyptian Hieroglyph N018a"
    @test udata(0x2FA00) == "CJK Compatibility Ideograph-2Fa00"
    @test udata(0xE005F) == "Tag Low Line"
end

ext: jl


text:
import GLAbstraction: gluniform

struct UniformCalculator <: System end

function Overseer.requested_components(::UniformCalculator)
    return (Spatial, Shape, ModelMat, Dynamic, Camera3D, UpdatedComponents, Rotation)
end

function Overseer.update(::UniformCalculator, m::AbstractLedger)
    uc        = singleton(m, UpdatedComponents)
    m_updated = false
    modelmat  = m[ModelMat]
    dyn       = m[Dynamic]
    camera    = m[Camera3D]
    spatial   = m[Spatial]
    shape     = m[Shape]
    rotation  = m[Rotation]

    for e in @entities_in(spatial)
        if !in(e, modelmat) ||
           in(e, dyn) ||
           in(e, camera) ||
           in(Spatial, uc) ||
           in(Rotation, uc)
            m_updated = true

            tmat = translmat(spatial[e].position)
            if in(e, rotation)
                tmat = tmat * Mat4f0(rotation[e].q)
            end

            if in(e, shape)
                tmat = tmat * scalemat(shape[e].scale)
            end
            modelmat[e] = ModelMat(tmat)
        end
    end
    if m_updated
        push!(uc, ModelMat)
    end
end

# function Overseer.update(::UniformCalculator, m::AbstractLedger)
#     uc        = m[UpdatedComponents][1]
#     m_updated = false
#     modelmat  = m[ModelMat]
#     dyn       = m[Dynamic]
#     camera    = m[Camera3D]
#     spatial   = m[Spatial]
#     shape     = m[Shape]
#     g1 = group(m, ModelMat, Spatial, Shape)
#     g2 = group(m, ModelMat, Spatial)
#     if length(g1) == 0
#         for e in @entities_in(spatial && shape)
#             g1[e] = ModelMat(translmat(spatial[e].position) * scalemat(Vec3f0(shape[e].scale)))
#         end
#         for e in @entities_in(spatial && !shape)
#             g2[e] = ModelMat(translmat(spatial[e].position))
#         end
#     else
#         @inbounds for i = 1:length(g2)
#             modelmat.data[i] = ModelMat(translmat(spatial[i].position))
#         end
#         @inbounds for i = 1:length(g1)
#             modelmat.data[i] = ModelMat(modelmat[i].modelmat * scalemat(Vec3f0(shape[i].scale)))
#         end
#     end
#     push!(uc, ModelMat)
# end

function set_entity_uniforms_func(render_program::Union{DefaultProgram,PeelingProgram,
                                                        LineProgram}, system::System)
    prog     = render_program.program
    material = component(system, Material)
    modelmat = component(system, ModelMat)
    ucolor   = component(system, UniformColor)
    return e -> begin
        gluniform(prog, :material, Vec2(material[e].specpow, material[e].specint))
        gluniform(prog, :modelmat, modelmat[e].modelmat)
        if has_entity(ucolor, e)
            gluniform(prog, :uniform_color, ucolor[e].color)
            gluniform(prog, :is_uniform, true)
        else
            gluniform(prog, :is_uniform, false)
            gluniform(prog, :specpow, material[e].specpow)
        end
    end
end

function gluniform(program::GLA.Program, spatial, camera::Camera3D)
    gluniform(program, :projview, camera.projview)
    return gluniform(program, :campos, spatial.position)
end

function gluniform(program::GLA.Program, pointlight::PointLight, color::UniformColor,
                   spatial::Spatial)
    gluniform(program, Symbol("plight.color"), RGB(color.color))
    gluniform(program, Symbol("plight.position"), spatial.position)
    gluniform(program, Symbol("plight.amb_intensity"), pointlight.ambient)
    gluniform(program, Symbol("plight.specular_intensity"), pointlight.specular)
    return gluniform(program, Symbol("plight.diff_intensity"), pointlight.diffuse)
end


text:
# a tutorial adapted from the fipy diffusion 1D example
# see: http://www.ctcms.nist.gov/fipy/examples/diffusion/index.html
# to run this file, first load JFVM and GR packages
# using JFVM, GR

# define the domain
L = 5  # domain length
Nx = 100 # number of cells
meshstruct = createMesh1D(Nx, L)
BC = createBC(meshstruct) # all Neumann boundary condition structure
BC.left.a[:] = 0.0
BC.left.b[:]=1.0
BC.left.c[:]=1.0 # left boundary
BC.right.a[:] = 0.0
BC.right.b[:]=1.0
BC.right.c[:]=0.0 # right boundary
x = meshstruct.cellcenters.x
## define the transfer coeffs
D_val = 1.0
alfa = createCellVariable(meshstruct, 1)
Dave = createFaceVariable(meshstruct, D_val)
## define initial values
c_old = createCellVariable(meshstruct, 0, BC) # initial values
c = c_old
## loop
dt = 0.001 # time step
final_t = 1.5
for t=dt:dt:final_t
    # step 1: calculate divergence term
    RHS = divergenceTerm(Dave.*gradientTerm(c_old))
    c = solveExplicitPDE(c_old, dt, RHS, BC, alfa)
    # step 2: calculate the new value for internal cells
#     c_old_int=internalCells(c_old)
#     c_val_int = dt*excludeGhostRHS(meshstruct, RHS+constantSourceTerm(c_old))+c_old_int[:]
    c_analytical = 1-erf(x/(2*sqrt(D_val*t)))
    GR.plot(x, [c.value[2:end-1] c_analytical])
    c_old = c
    print(t)
end


text:
# Inverse Wishart distribution
#
#   following Wikipedia parametrization
#

immutable InverseWishart{T<:Real, ST<:AbstractPDMat} <: ContinuousMatrixDistribution
    df::T     # degree of freedom
    Ψ::ST           # scale matrix
    c0::T     # log of normalizing constant
end

#### Constructors

function InverseWishart{T<:Real}(df::T, Ψ::AbstractPDMat{T})
    p = dim(Ψ)
    df > p - 1 || error("df should be greater than dim - 1.")
    c0 = _invwishart_c0(df, Ψ)
    R = Base.promote_eltype(T, c0)
    prom_Ψ = convert(AbstractArray{R}, Ψ)
    InverseWishart{R, typeof(prom_Ψ)}(R(df), prom_Ψ, R(c0))
end

function InverseWishart(df::Real, Ψ::AbstractPDMat)
    T = Base.promote_eltype(df, Ψ)
    InverseWishart(T(df), convert(AbstractArray{T}, Ψ))
end

InverseWishart(df::Real, Ψ::Matrix) = InverseWishart(df, PDMat(Ψ))

InverseWishart(df::Real, Ψ::Cholesky) = InverseWishart(df, PDMat(Ψ))

function _invwishart_c0(df::Real, Ψ::AbstractPDMat)
    h_df = df / 2
    p = dim(Ψ)
    h_df * (p * typeof(df)(logtwo) - logdet(Ψ)) + logmvgamma(p, h_df)
end

#### Properties

insupport(::Type{InverseWishart}, X::Matrix) = isposdef(X)
insupport(d::InverseWishart, X::Matrix) = size(X) == size(d) && isposdef(X)

dim(d::InverseWishart) = dim(d.Ψ)
size(d::InverseWishart) = (p = dim(d); (p, p))
params(d::InverseWishart) = (d.df, d.Ψ, d.c0)
@inline partype{T<:Real}(d::InverseWishart{T}) = T

### Conversion
function convert{T<:Real}(::Type{InverseWishart{T}}, d::InverseWishart)
    P = Wishart{T}(d.Ψ)
    InverseWishart{T, typeof(P)}(T(d.df), P, T(d.c0))
end
function convert{T<:Real}(::Type{InverseWishart{T}}, df, Ψ::AbstractPDMat, c0)
    P = Wishart{T}(Ψ)
    InverseWishart{T, typeof(P)}(T(df), P, T(c0))
end

#### Show

show(io::IO, d::InverseWishart) = show_multline(io, d, [(:df, d.df), (:Ψ, full(d.Ψ))])

#### Statistics

function mean(d::InverseWishart)
    df = d.df
    p = dim(d)
    r = df - (p + 1)
    if r > 0.0
        return full(d.Ψ) * (1.0 / r)
    else
        error("mean only defined for df > p + 1")
    end
end

mode(d::InverseWishart) = d.Ψ * inv(d.df + dim(d) + 1.0)

#### Evaluation

function _logpdf(d::InverseWishart, X::AbstractMatrix)
    p = dim(d)
    df = d.df
    Xcf = cholfact(X)
    # we use the fact: trace(Ψ * inv(X)) = trace(inv(X) * Ψ) = trace(X \ Ψ)
    Ψ = full(d.Ψ)
    -0.5 * ((df + p + 1) * logdet(Xcf) + trace(Xcf \ Ψ)) - d.c0
end

#### Sampling

rand(d::InverseWishart) = inv(cholfact!(rand(Wishart(d.df, inv(d.Ψ)))))

function _rand!{M<:Matrix}(d::InverseWishart, X::AbstractArray{M})
    wd = Wishart(d.df, inv(d.Ψ))
    for i in 1:length(X)
        X[i] = inv(cholfact!(rand(wd)))
    end
    return X
end


text:
using ManifoldsBase
using LinearAlgebra
using Test

struct CustomValidationManifoldRetraction <: ManifoldsBase.AbstractRetractionMethod end

function ManifoldsBase.injectivity_radius(
    ::ManifoldsBase.DefaultManifold,
    ::CustomValidationManifoldRetraction,
)
    return 10.0
end
function ManifoldsBase.injectivity_radius(
    ::ManifoldsBase.DefaultManifold,
    p,
    ::CustomValidationManifoldRetraction,
)
    return 11.0
end

@testset "Validation manifold" begin
    M = ManifoldsBase.DefaultManifold(3)
    A = ValidationManifold(M)
    x = [1.0, 0.0, 0.0]
    y = 1 / sqrt(2) * [1.0, 1.0, 0.0]
    z = [0.0, 1.0, 0.0]
    v = log(M, x, y)
    x2 = ValidationMPoint(x)
    y2 = ValidationMPoint(y)
    v2 = log(A, x, y) # auto convert
    y2 = exp(A, x, v2)
    w = log(M, x, z)
    w2 = log(A, x, z; atol = 10^(-15))
    @testset "Checks forward" begin
        @test ManifoldsBase.check_size(A, x2, v2) === ManifoldsBase.check_size(M, x, v)
        @test ManifoldsBase.check_size(A, x2) === ManifoldsBase.check_size(M, x)
    end
    @testset "Types and Conversion" begin
        @test convert(typeof(M), A) == M
        @test convert(typeof(A), M) == A
        @test base_manifold(A) == M
        @test base_manifold(base_manifold(A)) == base_manifold(A)
        @test ManifoldsBase.representation_size(A) == ManifoldsBase.representation_size(M)
        @test ManifoldsBase.representation_size(A) == (3,)
        @test manifold_dimension(A) == manifold_dimension(M)
        @test manifold_dimension(A) == 3
        for T in [ValidationMPoint, ValidationTVector, ValidationCoTVector]
            p = T(x)
            @test convert(typeof(x), p) == x
            @test convert(typeof(p), y) == T(y)
            @test number_eltype(typeof(p)) == eltype(x)
            @test number_eltype(p) == eltype(x)
            @test typeof(allocate(p)) == typeof(p)
            @test typeof(allocate(p, eltype(x))) == typeof(p)
            @test allocate(p) isa T
            @test allocate(p, Float32) isa T
            @test number_eltype(allocate(p, Float32)) == Float32
            @test similar(p) isa T
            @test similar(p, Float32) isa T
            @test number_eltype(similar(p, Float32)) == Float32
            q = allocate(p)
            if T == ValidationMPoint
                copyto!(A, q, p)
            else
                copyto!(A, q, ValidationMPoint(x), p) # generate base point “on the fly”
            end
            @test isapprox(A, q, p)
            @test ManifoldsBase.array_value(p) == x
            @test ManifoldsBase.array_value(x) == x
        end
    end
    @testset "Vector functions" begin
        for T in [ValidationTVector, ValidationCoTVector]
            a = T(v)
            b = T(w)
            @test isapprox(A, a + b, T(v + w))
            @test isapprox(A, (a - b), T(v - w))
            @test isapprox(A, -b, T(-w))
            @test isapprox(A, 2 * a, T(2 .* v))
            @test isapprox(A, 2 .* a .+ b, T(2 .* v .+ w))
            c = similar(a)
            c .= a .+ b
            @test isapprox(A, c, a .+ b)
        end
    end
    @testset "AbstractManifold functions" begin
        @test manifold_dimension(A) == manifold_dimension(M)
        @test isapprox(y2.value, y)
        @test distance(A, x, y) == distance(M, x, y)
        @test norm(A, x, v) == norm(M, x, v)
        @test inner(A, x, v2, w2; atol = 10^(-15)) == inner(M, x, v, w)
        @test isapprox(A, x2, y2) == isapprox(M, x, y)
        @test isapprox(A, x, y) == isapprox(A, x2, y2)
        @test isapprox(A, x, v2, v2) == isapprox(M, x, v, v)
        v2s = similar(v2)
        project!(A, v2s, x2, v2)
        @test isapprox(A, v2, v2s)
        y2s = similar(y2)
        exp!(A, y2s, x2, v2)
        @test isapprox(A, y2s, y2)
        log!(A, v2s, x, y)
        @test mid_point(A, x, y) == 0.5 * (x .+ y)
        mp = similar(x)
        mid_point!(A, mp, x, y)
        @test mp == 0.5 * (x .+ y)
        @test isapprox(A, x, v2s, v2)
        @test isapprox(A, exp(A, x, v), y2)
        @test isapprox(A, zero_vector(A, x), zero_vector(M, x))
        vector_transport_to!(A, v2s, x2, v2, y2)
        @test isapprox(A, x2, v2, v2s)
        @test isapprox(A, x2, v2, vector_transport_to(A, x2, v2, y2))
        slt = ManifoldsBase.SchildsLadderTransport()
        vector_transport_to!(A, v2s, x2, v2, y2, slt)
        @test isapprox(A, x2, v2, v2s)
        @test isapprox(A, x2, v2, vector_transport_to(A, x2, v2, y2, slt))
        plt = ManifoldsBase.PoleLadderTransport()
        vector_transport_to!(A, v2s, x2, v2, y, plt)
        @test isapprox(A, x2, v2, vector_transport_to(A, x2, v2, y2, plt))
        @test isapprox(A, x2, v2, v2s)
        pt = ManifoldsBase.ProjectionTransport()
        vector_transport_to!(A, v2s, x2, v2, y2, pt)
        @test isapprox(A, x2, v2, v2s)
        @test isapprox(A, x2, v2, vector_transport_to(A, x2, v2, y2, pt))
        zero_vector!(A, v2s, x)
        @test isapprox(A, x, v2s, zero_vector(M, x))
        c2 = [x2]
        v3 = similar(v2)
        @test isapprox(A, x2, v2, vector_transport_along!(A, v3, x2, v2, c2))
        @test isapprox(A, x2, v2, vector_transport_along(A, x2, v2, c2))
        @test injectivity_radius(A) == Inf
        @test injectivity_radius(A, x) == Inf
        @test injectivity_radius(A, ManifoldsBase.ExponentialRetraction()) == Inf
        @test injectivity_radius(A, x, ManifoldsBase.ExponentialRetraction()) == Inf
        @test injectivity_radius(A, CustomValidationManifoldRetraction()) == 10
        @test injectivity_radius(A, x, CustomValidationManifoldRetraction()) == 11
    end

    @testset "ValidationManifold basis" begin
        b = [Matrix(I, 3, 3)[:, i] for i in 1:3]
        for BT in (DefaultBasis, DefaultOrthonormalBasis, DefaultOrthogonalBasis)
            @testset "Basis $(BT)" begin
                cb = BT()
                @test b == get_vectors(M, x, get_basis(A, x, cb))
                v = similar(x)
                v2 = similar(x)
                @test_throws ErrorException get_vector(A, x, [1.0], cb)
                @test_throws DomainError get_vector(A, [1.0], [1.0, 0.0, 0.0], cb)
                @test_throws ErrorException get_vector!(A, v, x, [], cb)
                @test_throws DomainError get_vector!(A, v, [1.0], [1.0, 0.0, 0.0], cb)
                @test_throws DomainError get_coordinates(A, x, [1.0], cb)
                @test_throws DomainError get_coordinates!(A, v, x, [], cb)
                @test_throws DomainError get_coordinates!(A, v, [1.0], [1.0, 0.0, 0.0], cb)
                @test get_vector(A, x, [1, 2, 3], cb) ≈ get_vector(M, x, [1, 2, 3], cb)
                @test get_vector!(A, v2, x, [1, 2, 3], cb) ≈
                      get_vector!(M, v, x, [1, 2, 3], cb)
                @test get_coordinates(A, x, [1, 2, 3], cb) ≈
                      get_coordinates(M, x, [1, 2, 3], cb)
                @test get_coordinates!(A, v2, x, [1, 2, 3], cb) ≈
                      get_coordinates!(M, v, x, [1, 2, 3], cb)

                @test_throws ErrorException get_basis(A, x, CachedBasis(cb, [x]))
                @test_throws ErrorException get_basis(A, x, CachedBasis(cb, [x, x, x]))
                @test_throws ErrorException get_basis(A, x, CachedBasis(cb, [2 * x, x, x]))
                if BT <: ManifoldsBase.AbstractOrthogonalBasis
                    @test_throws ArgumentError get_basis(
                        A,
                        x,
                        CachedBasis(
                            cb,
                            [[1.0, 0.0, 0.0], [1.0, 1.0, 0.0], [0.0, 0.0, 1.0]],
                        ),
                    )
                elseif BT <: ManifoldsBase.AbstractOrthonormalBasis
                    @test_throws ArgumentError get_basis(
                        A,
                        x,
                        CachedBasis(
                            cb,
                            [[2.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]],
                        ),
                    )
                end
            end
        end
    end
end


text:
# Inspired from https://github.com/JuliaDatabases/MySQL.jl/blob/36eaf2bfbbdd9a27c408d0b2a734fff0d81b63ad/deps/build.jl
module Anon1 end
module Anon2 end
module Anon3 end

@static if VERSION < v"1.3.0"

using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))

products = [
    LibraryProduct(prefix, ["libcsdp"], :libcsdp)
]

Anon1.include("build_CompilerSupportLibraries.v0.3.3.jl")
Anon2.include("build_OpenBLAS32.v0.3.9.jl")
Anon3.include("build_CSDP.v6.2.0.jl")

# Finally, write out a deps.jl file
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=true)

end # VERSION


text:
function loadDiabetes(path)
    cd(path)
    mainData = CSV.read("Diabetes/diabetes.csv", delim = ";", nullable=false)
    return mainData
end

function loadWineQuality(path)
    cd(path)
    mainData = CSV.read("WineQualityRed/winequality-red.csv", delim = ";", nullable=false, types=fill(Float64,12))
    return mainData
end

function loadHousingData(path)
    cd(path)
    mainData = CSV.read("AmesHousingModClean.csv", delim = ';', nullable=false)
    return mainData
end

function loadConcrete(path)
    cd(path)
    mainData = CSV.read("Concrete_Data.csv", delim = ',', nullable=false, types=Dict(3=>Float64))
    return mainData
end

function loadIndexDataNoDur(path)
    cd(path)
    mainData = CSV.read("monthlyNoDurReturn.csv", header=["NoDur","Durbl",
        "Manuf","Enrgy","HiTec","Telcm","Shops", "Hlth","Utils","Other","Index",
        "D12","E12","b.m","tbl","AAA","BAA","lty","ntis","Rfree","infl","ltr",
        "corpr","svar","csp","CRSP_SPvw","CRSP_SPvwx","Ycol"],
        delim = ',', nullable=false, types=Dict(25=>Float64))
    return mainData
end

function loadIndexDataNoDurLOGReturn(path)
    cd(path)
    mainData = CSV.read("monthlyNoDurLOGReturn.csv", header=["NoDur","Durbl",
        "Manuf","Enrgy","HiTec","Telcm","Shops", "Hlth","Utils","Other","Index",
        "D12","E12","b.m","tbl","AAA","BAA","lty","ntis","Rfree","infl","ltr",
        "corpr","svar","csp","CRSP_SPvw","CRSP_SPvwx","VIX","Ycol","Date","Resession"],
        delim = ',', nullable=false, types=Dict(25=>Float64, 28=>Float64))
    return mainData
end

#=  For index data where index is absolute value and not diff
function loadIndexDataLOGReturn(industry, path)
    cd(path)
    mainData = CSV.read("monthly"*industry*"LOGReturn.csv", header=["NoDur","Durbl",
        "Manuf","Enrgy","HiTec","Telcm","Shops", "Hlth","Utils","Other","Index",
        "D12","E12","b.m","tbl","AAA","BAA","lty","ntis","Rfree","infl","ltr",
        "corpr","svar","csp","CRSP_SPvw","CRSP_SPvwx","VIX","Ycol","Date","Resession"],
        delim = ',', nullable=false, types=Dict(25=>Float64, 28=>Float64))
    return mainData
end
=#
function loadIndexDataLOGReturn(industry, path)
    cd(path)
    mainData = CSV.read("monthly"*industry*"LOGReturn2.csv", header=["NoDur","Durbl",
        "Manuf","Enrgy","HiTec","Telcm","Shops", "Hlth","Utils","Other","Index",
        "D12","E12","b.m","tbl","AAA","BAA","lty","ntis","Rfree","infl","ltr",
        "corpr","svar","csp","CRSP_SPvw","CRSP_SPvwx","VIX","Ycol","Date","Resession"],
        delim = ',', nullable=false, types=Dict(25=>Float64, 28=>Float64))
    return mainData
end

function loadRiskFreeRate(industry, path)
    cd(path)
    mainData = CSV.read("riskFreeRateLOG.csv",header =["RfreeRate"],delim = ',', nullable = false)
    return mainData
end

function loadIndexDataNoDurVIX(path)
    cd(path)
    mainData = CSV.read("monthlyNoDurReturnVIX.csv", header=["NoDur","Durbl",
        "Manuf","Enrgy","HiTec","Telcm","Shops", "Hlth","Utils","Other","Index",
        "D12","E12","b.m","tbl","AAA","BAA","lty","ntis","Rfree","infl","ltr",
        "corpr","svar","csp","CRSP_SPvw","CRSP_SPvwx","VIX","Ycol"],
        delim = ',', nullable=false, types=Dict(25=>Float64))
    return mainData
end

function loadIndexDataOther(path)
    cd(path)
    mainData = CSV.read("monthlyOtherReturn.csv", header=["NoDur","Durbl",
        "Manuf","Enrgy","HiTec","Telcm","Shops", "Hlth","Utils","Other","Index",
        "D12","E12","b.m","tbl","AAA","BAA","lty","ntis","Rfree","infl","ltr",
        "corpr","svar","csp","CRSP_SPvw","CRSP_SPvwx","Ycol"],
        delim = ',', nullable=false, types=Dict(25=>Float64))
    return mainData
end

function loadIndexDataDailyNoDur(path)
    cd(path)
    mainData = CSV.read("dailyNoDurReturn.csv", header=["NoDur","Durbl",
        "Manuf","Enrgy","HiTec","Telcm","Shops", "Hlth","Utils","Other","Index"],
        delim = ',', nullable=false, types=Dict(11=>Float64))
    return mainData
end

function loadIndexDataDailyDurbl(path)
    cd(path)
    mainData = CSV.read("dailyDurblReturn.csv", header=["NoDur","Durbl",
        "Manuf","Enrgy","HiTec","Telcm","Shops", "Hlth","Utils","Other","Index"],
        delim = ',', nullable=false, types=Dict(11=>Float64))
    return mainData
end

function loadIndexDataDailyManuf(path)
    cd(path)
    mainData = CSV.read("dailyManufReturn.csv", header=["NoDur","Durbl",
        "Manuf","Enrgy","HiTec","Telcm","Shops", "Hlth","Utils","Other","Index"],
        delim = ',', nullable=false, types=Dict(11=>Float64))
    return mainData
end

function loadIndexDataDailyEnrgy(path)
    cd(path)
    mainData = CSV.read("dailyEnrgyReturn.csv", header=["NoDur","Durbl",
        "Manuf","Enrgy","HiTec","Telcm","Shops", "Hlth","Utils","Other","Index"],
        delim = ',', nullable=false, types=Dict(11=>Float64))
    return mainData
end

function loadIndexDataDailyHiTec(path)
    cd(path)
    mainData = CSV.read("dailyHiTecReturn.csv", header=["NoDur","Durbl",
        "Manuf","Enrgy","HiTec","Telcm","Shops", "Hlth","Utils","Other","Index"],
        delim = ',', nullable=false, types=Dict(11=>Float64))
    return mainData
end

function loadIndexDataDailyTelcm(path)
    cd(path)
    mainData = CSV.read("dailyTelcmReturn.csv", header=["NoDur","Durbl",
        "Manuf","Enrgy","HiTec","Telcm","Shops", "Hlth","Utils","Other","Index"],
        delim = ',', nullable=false, types=Dict(11=>Float64))
    return mainData
end

function loadIndexDataDailyShops(path)
    cd(path)
    mainData = CSV.read("dailyShopsReturn.csv", header=["NoDur","Durbl",
        "Manuf","Enrgy","HiTec","Telcm","Shops", "Hlth","Utils","Other","Index"],
        delim = ',', nullable=false, types=Dict(11=>Float64))
    return mainData
end

function loadIndexDataDailyHlth(path)
    cd(path)
    mainData = CSV.read("dailyHlthReturn.csv", header=["NoDur","Durbl",
        "Manuf","Enrgy","HiTec","Telcm","Shops", "Hlth","Utils","Other","Index"],
        delim = ',', nullable=false, types=Dict(11=>Float64))
    return mainData
end

function loadIndexDataDailyUtils(path)
    cd(path)
    mainData = CSV.read("dailyUtilsReturn.csv", header=["NoDur","Durbl",
        "Manuf","Enrgy","HiTec","Telcm","Shops", "Hlth","Utils","Other","Index"],
        delim = ',', nullable=false, types=Dict(11=>Float64))
    return mainData
end

function loadIndexDataDailyOther(path)
    cd(path)
    mainData = CSV.read("dailyOtherReturn.csv", header=["NoDur","Durbl",
        "Manuf","Enrgy","HiTec","Telcm","Shops", "Hlth","Utils","Other","Index"],
        delim = ',', nullable=false, types=Dict(11=>Float64))
    return mainData
end

function loadCPUData(path)
    cd(path)
    mainData = CSV.read("machine.csv", header=["MYCT",
    	"MMIN","MMAX","CACH","CHMIN","CHMAX","PRP"], datarow=1, nullable=false)
    return mainData
end

function loadElevatorData(path)
    cd(path)
    mainData = CSV.read("Elevators/elevators.data", header=["climbRate", "Sgz", "p", "q", "curRoll", "absRoll", "diffClb",
    	"diffRollRate", "diffDiffClb", "SaTime1", "SaTime2", "SaTime3", "SaTime4", "diffSaTime1", "diffSaTime2",
    	"diffSaTime3", "diffSaTime4", "Sa", "Goal"], datarow=1, nullable=false)
    testData = CSV.read("Elevators/elevators.test", header=["climbRate", "Sgz", "p", "q", "curRoll", "absRoll", "diffClb",
    	"diffRollRate", "diffDiffClb", "SaTime1", "SaTime2", "SaTime3", "SaTime4", "diffSaTime1", "diffSaTime2",
    	"diffSaTime3", "diffSaTime4", "Sa", "Goal"], datarow=1, nullable=false)
    return mainData, testData
end

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_name: skipperafrosenborg/Thesis


text:
# simple program to test affinity propagation

using Base.Test
using Distances
using Clustering

srand(34568)

d = 10
n = 500
x = rand(d, n)
S = -pairwise(Euclidean(), x, x)

# set diagonal value to median value
S = S - diagm(diag(S)) + median(S)*eye(size(S,1))

R = affinityprop(S)

@test isa(R, AffinityPropResult)
k = length(R.exemplars)
@test 0 < k < n
@test length(R.assignments) == n
@test all(R.assignments .>= 1) && all(R.assignments .<= k)
@test all(R.assignments[R.exemplars] .== collect(1:k))

@test length(R.counts) == k
@test sum(R.counts) == n
for i = 1:k
    @test R.counts[i] == countnz(R.assignments .== i)
end


text:
pointstr, foldstr = split(read("input.txt", String), "\n\n")
points = map(x -> parse.(Int, split(x, ',')), split.(pointstr, '\n', keepempty=false))
folds = map(x -> split(x, "="), split.(foldstr, '\n', keepempty=false))

function fold!(points, directions)
    dim = directions[1] == "fold along x" ? 1 : 2
    mag = parse(Int, directions[2])
    for p in points
        Δ = p[dim] - mag
        Δ > 0 && (p[dim] -= 2*Δ)
    end
    return points
end

function display_code(points)
    points |> sort! |> unique!
    X, Y = maximum(reduce(hcat, points), dims=2)
    for y in 0:Y
        println()
        for x in 0:X
            print([x, y] in points ? repeat(Char(0x2588), 2) : "  ")
        end
    end
end

for d in folds
    fold!(points, d)
    d == first(folds) && @info "Part 1: $(length(unique(points)))"
end
display_code(points)


text:
using OrdinaryDiffEq

export ode_solver

"""
    ode_solver(f!, k, ts::Vector{Float64}, x0::Vector{Float64}, method; savedt=true)

    Solves differential equations!

    Input
    f!		set of differential equations
    k		reactions rates for the network
    ts		time grid vector
    x0		initial condition
    method	integration algorithm. If not specified: default is Tsit5()

	You can find a list of available solvers at https://diffeq.sciml.ai/stable/solvers/ode_solve/

	savedt	If true, the result will be saved after each time step.
			If false, an adaptive time step is used.
			For method==Euler, this has to be true and dt will be set to the
			discretization in the time grid vector `ts`

    Returns the result in a struct with time and mean copy
    numbers of molecules.

"""
function ode_solver(f!, k, ts::Vector{Float64}, x0::Vector{Float64}, method; savedt=true)

	tspan = (ts[1], ts[end])
	prob = ODEProblem(f!, x0, tspan, k)

	if method==Euler
		println("Euler algoritm needs choice of dt, setting savedt=true")
		delta_t = (ts[2]-ts[1])
		sol = solve(prob, Euler(), dt=delta_t, saveat=delta_t)
	elseif savedt
		sol = solve(prob, method(), saveat=(ts[2]-ts[1]))
	else
		sol = solve(prob, method())
	end

	col = []
	for i in 1:length(x0)
	    push!(col, map(A->A[i], sol.u))
	end

    xar = reshape(cat(col..., dims=1), length(sol.t), length(x0))
    return Result(sol.t, xar')
end

function ode_solver(f!, k, tspan::Tuple, x0)
	println("ode_solver: This is a old version! Please use a vector for tspan
	and a specified method for the integration algorithm (default: Tsit5()). ")
	ts = [tspan[1], tspan[end]]
	ode_solver(f!, k, ts, x0, Tsit5; savedt=false)
end

function ode_solver(f!, k, ts::Vector{Float64}, x0::Vector{Float64})
	println("No integration algorithm specified! Using Tsit5() as default!")
	ode_solver(f!, k, ts, x0, Tsit5)
end

max_stars_repo_name: mkurjahn/CRN

max_issues_repo_name: mkurjahn/CRN

max_forks_repo_name: mkurjahn/CRN


text:
## Sensor model that returns position and velocities of other traffic participants
# Gaussian Noise
# FP/FN
# Late detection?
# Position noise depends on distance

"""
Sensor model that returns position and velocities of other traffic participants
Gaussian Noise
FP/FN
Late detection?
Position noise depends on distance
"""
@with_kw struct GaussianSensor{P <: AbstractNoiseModel, V <: AbstractNoiseModel} <: AbstractSensor
    pos_noise::P = LinearNoise()
    vel_noise::V = LinearNoise()
    false_positive_rate::Float64 = 0.1
    false_negative_rate::Float64 = 0.1
    rng::AbstractRNG = MersenneTwister(1)
end

function measure(sensor::GaussianSensor, ego::Vehicle, scene::Scene, roadway::Roadway, obstacles::Vector{ConvexPolygon})
    obs = Vector{Vehicle}()
    sizehint!(obs, 10) #XXX
    for veh in scene
        if veh.id == ego.id
            continue
        end
        zveh = measure(sensor, ego, veh, roadway, obstacles, sensor.rng)
        if zveh != nothing
            push!(obs, Vehicle(zveh, veh.def, veh.id))
        end
    end
    if isempty(obs)
        false_car = car_false_positive(sensor.false_positive_rate, roadway, scene, obstacles, sensor.rng)
        false_car == nothing ? nothing : push!(obs, false_car)
        false_ped = ped_false_positive(sensor.false_positive_rate, roadway, scene, obstacles, sensor.rng)
        false_ped == nothing ? nothing : push!(obs, false_ped)
    end
    return obs
end

function measure(sensor::GaussianSensor, ego::Vehicle, veh::Vehicle, roadway::Roadway, obstacles::Vector{ConvexPolygon}, rng::AbstractRNG)
    occluded = occlusion_checker(ego.state, veh.state, obstacles)
    FN = rand(rng) < sensor.false_negative_rate
    if occluded || FN
        return nothing
    else
        dist = sqrt((ego.state.posG.x - veh.state.posG.x)^2 + (ego.state.posG.y - veh.state.posG.y)^2)
        pos_std = noise(sensor.pos_noise, dist)
        vel_std = noise(sensor.vel_noise, dist)
        z_pos = VecSE2(veh.state.posG.x + pos_std*randn(rng),
                       veh.state.posG.y + pos_std*randn(rng),
                       veh.state.posG.θ)
        z_v = veh.state.v + vel_std*randn(rng)
        return VehicleState(z_pos, roadway, z_v)
    end
end

function obs_weight(sensor::GaussianSensor, ego::VehicleState, obs::Union{Nothing, VehicleState}, veh::Union{Nothing, VehicleState}, obstacles::Vector{ConvexPolygon})
    weight = 1.0
    visible = veh != nothing && !occlusion_checker(ego, veh, obstacles)
    if !visible
        if obs == nothing
            return weight*(1. - sensor.false_positive_rate)
        else
            return weight*sensor.false_positive_rate # should handle false positive
        end
    else
        if obs == nothing
            return sensor.false_positive_rate
        else
            dist = sqrt((ego.posG.x - veh.posG.x)^2 + (ego.posG.y - veh.posG.y)^2)
            pos_std = noise(sensor.pos_noise, dist)
            vel_std = noise(sensor.vel_noise, dist)
            weight *= (1 - sensor.false_positive_rate)
            weight *= 1/(sqrt(2*pi*pos_std*vel_std))*exp(-1/pos_std*(veh.posG.x - obs.posG.x)^2 -1/pos_std*(veh.posG.y - obs.posG.y)^2 - 1/vel_std*(veh.v - obs.v)^2)
            return weight
        end
    end
end

@with_kw struct GaussianSensorOverlay <: SceneOverlay
    sensor::GaussianSensor = GaussianSensor()
    o::Vector{Vehicle} = Vector{Vehicle}()
    color::Colorant = RGBA(0.976, 0.592, 0.122, 0.7) # orange
end

function AutoViz.render!(rendermodel::RenderModel, overlay::GaussianSensorOverlay, scene::Scene, roadway::R) where R
    for veh in overlay.o
        p = veh.state.posG
        add_instruction!(rendermodel, render_vehicle, (p.x, p.y, p.θ, veh.def.length, veh.def.width, overlay.color))
    end
    return rendermodel
end


text:
module QuickCheck2

using Base.Test
import Base.Test: Result, Pass, Fail, Error, record, get_testset

export property, condproperty, quantproperty
export @test_formany, @test_forall, @test_exists, @not

for (dir, filename) in [

		(".", "test-formany.jl"),
		(".", "test-exists.jl"),
		(".", "test-forall.jl"),
		(".", "parse-argument-data.jl"),
		(".", "custom-generator.jl"),
		(".", "quickcheck-original.jl"),
	]

include(joinpath(dir, filename))
end

end #of module


text:
# this part is just defined for debug convenience
using DiffRules, SpecialFunctions, NaNMath

for (M, f, arity) in DiffRules.diffrules()
    op = :($M.$f)

    if arity == 1
        @eval $op(x::Expression) = Term($op, x)
    elseif arity == 2
        @eval $op(x::Expression, y::Expression) = Term($op, x, y)
        @eval $op(x::Number, y::Expression) = Term($op, Numeric(x), y)
        @eval $op(x::Expression, y::Number) = Term($op, x, Numeric(y))

        @eval $op(x::Irrational, y::Expression) = Term($op, Constant(x), y)
        @eval $op(x::Expression, y::Irrational) = Term($op, x, Constant(y))
    end
end


text:
### Soil heat model

export SoilHeatModel, PrescribedTemperatureModel

abstract type AbstractHeatModel <: AbstractSoilComponentModel end

"""
    PrescribedTemperatureModel{F1} <: AbstractHeatModel

Model structure for a prescribed temperature model.
"""
struct PrescribedTemperatureModel{F1} <: AbstractHeatModel
    "Temperature"
    T::F1
end

"""
    PrescribedTemperatureModel(
        T::Function = (aux,t) -> eltype(aux)(0.0)
    )
Outer constructor for the PrescribedTemperatureModel defining default values.
The functions supplied by the user are point-wise evaluated and are
evaluated in the Balance Law functions compute_gradient_argument,
 nodal_update, etc. whenever the prescribed temperature content variables are
needed by the water model.
"""
function PrescribedTemperatureModel(T::Function = (aux, t) -> eltype(aux)(0.0))
    return PrescribedTemperatureModel{typeof(T)}(T)
end

"""
    SoilHeatModel{FT, FiT, BCT} <: AbstractHeatModel
The necessary components for the Heat Equation in a soil water matrix.
# Fields
$(DocStringExtensions.FIELDS)
"""
struct SoilHeatModel{FT, FiT, BCT} <: AbstractHeatModel
    "Initial conditions for temperature"
    initialT::FiT
    "Boundary Condition Type"
    boundaries::BCT
end

"""
    SoilHeatModel(
        ::Type{FT};
        initialT::FT = FT(NaN),
        boundaries::AbstractBoundaryConditions,
    ) where {FT}

Constructor for the SoilHeatModel.
"""
function SoilHeatModel(
    ::Type{FT};
    initialT = (aux) -> FT(NaN),
    boundaries::AbstractBoundaryConditions,
) where {FT}
    args = (initialT, boundaries)
    return SoilHeatModel{FT, typeof.(args)...}(args...)
end

"""
    function get_temperature(
        heat::SoilHeatModel
        aux::Vars,
        t::Real
    )
Returns the temperature when the heat model chosen is the dynamical SoilHeatModel.
This is necessary for fully coupling heat and water.
"""
function get_temperature(heat::SoilHeatModel, aux::Vars, t::Real)
    T = aux.soil.heat.T
    return T
end

"""
    function get_temperature(
        heat::PrescribedTemperatureModel,
        aux::Vars,
        t::Real
    )
Returns the temperature when the heat model chosen is a user prescribed one.
This is useful for driving Richard's equation without a back reaction on temperature.
"""
function get_temperature(heat::PrescribedTemperatureModel, aux::Vars, t::Real)
    T = heat.T(aux, t)
    return T
end

"""
    function get_initial_temperature(
        m::SoilHeatModel
        aux::Vars,
        t::Real
    )
Returns the temperature from the SoilHeatModel.
Needed for soil_init_aux! of SoilWaterModel.
"""
function get_initial_temperature(m::SoilHeatModel, aux::Vars, t::Real)
    return m.initialT(aux)
end

"""
    function get_initial_temperature(
        m::PrescribedTemperatureModel,
        aux::Vars,
        t::Real
    )
Returns the temperature from the prescribed model.
Needed for soil_init_aux! of SoilWaterModel.
"""
function get_initial_temperature(
    m::PrescribedTemperatureModel,
    aux::Vars,
    t::Real,
)
    return m.T(aux, t)
end

vars_state(heat::SoilHeatModel, st::Prognostic, FT) = @vars(ρe_int::FT)
vars_state(heat::SoilHeatModel, st::Auxiliary, FT) = @vars(T::FT)
vars_state(heat::SoilHeatModel, st::Gradient, FT) = @vars(T::FT)
vars_state(heat::SoilHeatModel, st::GradientFlux, FT) =
    @vars(κ∇T::SVector{3, FT})

function soil_init_aux!(
    land::LandModel,
    soil::SoilModel,
    heat::SoilHeatModel,
    aux::Vars,
    geom::LocalGeometry,
)
    aux.soil.heat.T = heat.initialT(aux)
end

function land_nodal_update_auxiliary_state!(
    land::LandModel,
    soil::SoilModel,
    heat::SoilHeatModel,
    state::Vars,
    aux::Vars,
    t::Real,
)

    ϑ_l, θ_i = get_water_content(land.soil.water, aux, state, t)
    eff_porosity = soil.param_functions.porosity - θ_i
    θ_l = volumetric_liquid_fraction(ϑ_l, eff_porosity)
    ρc_ds = soil.param_functions.ρc_ds
    ρc_s = volumetric_heat_capacity(θ_l, θ_i, ρc_ds, land.param_set)
    aux.soil.heat.T = temperature_from_ρe_int(
        state.soil.heat.ρe_int,
        θ_i,
        ρc_s,
        land.param_set,
    )
end

function compute_gradient_argument!(
    land::LandModel,
    soil::SoilModel,
    heat::SoilHeatModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)

    ϑ_l, θ_i = get_water_content(land.soil.water, aux, state, t)
    eff_porosity = soil.param_functions.porosity - θ_i
    θ_l = volumetric_liquid_fraction(ϑ_l, eff_porosity)
    ρc_ds = soil.param_functions.ρc_ds
    ρc_s = volumetric_heat_capacity(θ_l, θ_i, ρc_ds, land.param_set)
    transform.soil.heat.T = temperature_from_ρe_int(
        state.soil.heat.ρe_int,
        θ_i,
        ρc_s,
        land.param_set,
    )
end

function compute_gradient_flux!(
    land::LandModel,
    soil::SoilModel,
    heat::SoilHeatModel,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ϑ_l, θ_i = get_water_content(land.soil.water, aux, state, t)
    eff_porosity = soil.param_functions.porosity - θ_i
    θ_l = volumetric_liquid_fraction(ϑ_l, eff_porosity)
    κ_dry = k_dry(land.param_set, soil.param_functions)
    S_r = relative_saturation(θ_l, θ_i, soil.param_functions.porosity)
    kersten = kersten_number(θ_i, S_r, soil.param_functions)
    κ_sat = saturated_thermal_conductivity(
        θ_l,
        θ_i,
        soil.param_functions.κ_sat_unfrozen,
        soil.param_functions.κ_sat_frozen,
    )
    diffusive.soil.heat.κ∇T =
        thermal_conductivity(κ_dry, kersten, κ_sat) * ∇transform.soil.heat.T
end

function flux_second_order!(
    land::LandModel,
    soil::SoilModel,
    heat::SoilHeatModel,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    hyperdiffusive::Vars,
    aux::Vars,
    t::Real,
)
    ρe_int_l = volumetric_internal_energy_liq(aux.soil.heat.T, land.param_set)
    diffusive_water_flux =
        -ρe_int_l .* get_diffusive_water_term(soil.water, diffusive)
    diffusive_heat_flux = -diffusive.soil.heat.κ∇T
    flux.soil.heat.ρe_int += diffusive_heat_flux + diffusive_water_flux

end

max_forks_repo_name: jm-c/CLIMA


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

__precompile__(true)
module Qux
this = "Qux"
end

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


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

using Markdown
using InteractiveUtils

# ╔═╡ 63ba08cc-59a8-11eb-0a0f-27efac60d779
using Pkg, DrWatson

# ╔═╡ 6db218c6-59a8-11eb-2a8b-7107354cf590
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample
	using StatisticalRethinking
end

# ╔═╡ 51fc19b8-59a8-11eb-2214-15aca59b807b
md" ## Clip-08-07-12s.jl"

# ╔═╡ f346a9a1-de24-40ef-bf2f-eb046ed800d5
#@quickactivate "StatisticalRethinkingStan"

# ╔═╡ 8aaa4bcc-59a8-11eb-2003-f1213b116565
begin
	df = CSV.read(sr_datadir("rugged.csv"), DataFrame)
	dropmissing!(df, :rgdppc_2000)
	dropmissing!(df, :rugged)
	df.log_gdp = log.(df[:, :rgdppc_2000])
	df.log_gdp_s = df.log_gdp / mean(df.log_gdp)
	df.rugged_s = df.rugged / maximum(df.rugged)
	df.cid = [df.cont_africa[i] == 1 ? 1 : 2 for i in 1:size(df, 1)]
	data = (N = size(df, 1), K = length(unique(df.cid)),
		G = df.log_gdp_s, R = df.rugged_s, cid=df.cid)
	PRECIS(df[:, [:log_gdp, :log_gdp_s, :rugged, :rugged_s, :cid]])
end

# ╔═╡ 4c3738f0-ba84-456a-a53f-c2690be9c077
df[1:8, [:log_gdp, :rugged, :cont_africa, :log_gdp_s, :rugged_s]]

# ╔═╡ 60053fa4-f1f9-4754-825e-170a095a284d
md"
```
8×5 DataFrame
 Row │ log_gdp  rugged   cont_africa  log_gdp_std  rugged_std
     │ Float64  Float64  Int64        Float64      Float64
─────┼────────────────────────────────────────────────────────
   1 │ 7.49261    0.858            1     1.00276    0.138342
   2 │ 6.43238    1.78             1     0.860869   0.287004
   3 │ 6.86612    0.141            1     0.918918   0.0227346
   4 │ 6.90617    0.236            1     0.924278   0.0380522
   5 │ 8.9493     0.181            1     1.19772    0.0291841
   6 │ 7.04582    0.197            1     0.942968   0.0317639
   7 │ 7.36241    0.224            1     0.985338   0.0361174
   8 │ 7.54046    0.515            1     1.00917    0.0830377
```"

# ╔═╡ e890c59c-6a68-11eb-0d33-4d52a21d0ccf
stan8_1 = "
data {
	int N;
	vector[N] G;
	vector[N] R;
}

parameters {
	real a;
	real b;
	real<lower=0> sigma;
}

transformed parameters {
	vector[N] mu;
	mu = a + b * (R - 0.125);
}

model {
	a ~ normal(1, 0.1);
	b ~ normal(0, 0.3);
	sigma ~ exponential(1);
	G ~ normal(mu, sigma);
}
generated quantities {
	vector[N] log_lik;
	for (i in 1:N)
		log_lik[i] = normal_lpdf(G[i] | mu[i], sigma);
}
";

# ╔═╡ 8e815d5c-6a66-11eb-21dc-99734c31f4e1
begin
	m8_1s = SampleModel("m8.1s", stan8_1)
	rc8_1s = stan_sample(m8_1s; data)
	if success(rc8_1s)
		post8_1s_df = read_samples(m8_1s, :dataframe)
		PRECIS(post8_1s_df[:, [:a, :b, :sigma]])
	end
end

# ╔═╡ 04a39568-6a6c-11eb-19ea-f55be3914a51
begin
	x = range(0, stop=1, length=50)
	plot(;leg=false)
	for i in 1:10:4000
		y = post8_1s_df.a[i] .+ post8_1s_df.b[i] .* x
		plot!(x, y, color=:lightgrey)
	end
	plot!(x, mean(post8_1s_df.a) .+ mean(post8_1s_df.b) .* x, color=:darkblue)
end

# ╔═╡ 7e9e3674-6aea-11eb-3206-c1f775fe2d78
stan8_2 = "
data {
	int N;
	int K;
	vector[N] G;
	vector[N] R;
	int cid[N];
}

parameters {
	vector[K] a;
	real b;
	real<lower=0> sigma;
}

transformed parameters {
	vector[N] mu;
	for (i in 1:N)
		mu[i] = a[cid[i]] + b * (R[i] - 0.215);
}

model {
	a ~ normal(1, 0.1);
	b ~ normal(0, 0.3);
	sigma ~ exponential(1);
	G ~ normal(mu, sigma);
}
generated quantities {
	vector[N] log_lik;
for (i in 1:N)
	log_lik[i] = normal_lpdf(G[i] | mu[i], sigma);
}
";

# ╔═╡ 12db71d0-6aeb-11eb-341c-e3243ba74862
begin
	m8_2s = SampleModel("m8.2s", stan8_2)
	rc8_2s = stan_sample(m8_2s; data)
	if success(rc8_2s)
		post8_2s_df = read_samples(m8_2s, :dataframe)
		PRECIS(post8_2s_df[:, [Symbol("a.1"), Symbol("a.2"), :b, :sigma]])
	end
end

# ╔═╡ 8b889740-7142-11eb-3d74-b1b9dc6f3282
plot_models([m8_1s, m8_2s], [:a, Symbol("a.1"), Symbol("a.2"), :b, :sigma])

# ╔═╡ 4a5070f8-6af4-11eb-31ca-a7a3ae504821
compare([m8_1s, m8_2s], :waic)

# ╔═╡ 4fa4b38c-7069-11eb-1202-7552c5dcd671
plot_models([m8_1s, m8_2s], :waic)

# ╔═╡ f6c60a90-6e29-11eb-1133-99ad7b37e0af
compare([m8_1s, m8_2s], :psis)

# ╔═╡ 7244381c-6e2b-11eb-2373-1d6c5d27de80
md"
```

> PSIS(m8.1)
       PSIS    lppd  penalty  std_err
1 -188.6114 94.3057 2.755778 13.35391

> PSIS(m8.2)
       PSIS     lppd  penalty  std_err
1 -252.1136 126.0568 4.272533 15.29773

> compare( m8.1 , m8.2 )
       WAIC   SE dWAIC   dSE pWAIC weight
m8.2 -252.3 15.3   0.0    NA   4.3      1
m8.1 -188.7 13.3  63.5 15.15   2.7      0

> compare( m8.1 , m8.2, func=PSIS )
Some Pareto k values are high (>0.5). Set pointwise=TRUE to inspect individual points.
       PSIS    SE dPSIS   dSE pPSIS weight
m8.2 -252.1 15.28   0.0    NA   4.3      1
m8.1 -188.7 13.40  63.4 15.16   2.7      0
```
"

# ╔═╡ dad9e286-6e2c-11eb-1e89-9f3cb4e2d06b
l, ls, pk = psisloo(m8_1s)

# ╔═╡ 133600ce-6e2d-11eb-3030-0bad514887eb
sum(ls)

# ╔═╡ 943f0ba2-6af4-11eb-2698-b3e2be715e88
begin
	nt8_2s = read_samples(m8_2s)
	hpdi(nt8_2s.a[1,:] .- nt8_2s.a[2, :])
end

# ╔═╡ 45767e2e-6a63-11eb-3e12-354a7e32a374
md" ## End of clip-08-07-12s.jl"

# ╔═╡ Cell order:
# ╟─51fc19b8-59a8-11eb-2214-15aca59b807b
# ╠═63ba08cc-59a8-11eb-0a0f-27efac60d779
# ╠═f346a9a1-de24-40ef-bf2f-eb046ed800d5
# ╠═6db218c6-59a8-11eb-2a8b-7107354cf590
# ╠═8aaa4bcc-59a8-11eb-2003-f1213b116565
# ╠═4c3738f0-ba84-456a-a53f-c2690be9c077
# ╟─60053fa4-f1f9-4754-825e-170a095a284d
# ╠═e890c59c-6a68-11eb-0d33-4d52a21d0ccf
# ╠═8e815d5c-6a66-11eb-21dc-99734c31f4e1
# ╠═04a39568-6a6c-11eb-19ea-f55be3914a51
# ╠═7e9e3674-6aea-11eb-3206-c1f775fe2d78
# ╠═12db71d0-6aeb-11eb-341c-e3243ba74862
# ╠═8b889740-7142-11eb-3d74-b1b9dc6f3282
# ╠═4a5070f8-6af4-11eb-31ca-a7a3ae504821
# ╠═4fa4b38c-7069-11eb-1202-7552c5dcd671
# ╠═f6c60a90-6e29-11eb-1133-99ad7b37e0af
# ╟─7244381c-6e2b-11eb-2373-1d6c5d27de80
# ╠═dad9e286-6e2c-11eb-1e89-9f3cb4e2d06b
# ╠═133600ce-6e2d-11eb-3030-0bad514887eb
# ╠═943f0ba2-6af4-11eb-2698-b3e2be715e88
# ╟─45767e2e-6a63-11eb-3e12-354a7e32a374


text:
module Models

import Base: +, -, *, isapprox
using LinearAlgebra
using FlightMechanics

export
# Attitude
Attitude, get_euler_angles, get_quaternions,
# Position
Position, PositionLLH, PositionEarth, PositionECEF,
get_position, get_llh, get_xyz_earth, get_xyz_ecef, get_height,
# State
State,
get_attitude,
get_body_velocity, get_horizon_velocity, get_flight_path_angle,
get_body_ang_velocity, get_turn_rate,
get_body_accel, get_horizon_accel,
get_body_ang_accel,
get_sixdof_euler_fixed_mass_state,
# Mass
RigidSolid, PointMass, get_mass, get_cg, get_inertia,
# Pont Forces and moments
PointForcesMoments, translate_forces_moments, rotate,
# Atmosphere
AtmosphereISA,AtmosphereF16,
get_temperature, get_pressure, get_density, get_sound_velocity,
calculate_atmosphere,
# Wind
Wind, ConstantWind,
get_wind_NED, get_wind_body, get_wind_direction, get_wind_intensity,
get_vertical, calculate_wind,
# Gravity
EarthConstantGravity, get_gravity_horizon, get_gravity_accel, get_gravity_body,
calculate_gravity,
# Environment
Environment, calculate_environment, get_atmos, get_wind, get_gravity,
# Aerostate
AeroState, get_alpha, get_beta, get_aero_angles, get_alpha_dot, get_tas,
get_eas, get_cas, get_ias, get_aero_speeds, get_qinf, get_mach,
generate_state_aerostate, check_state_aerostate_env_coherence,
# Controls
Control, RangeControl, DiscreteControl,
get_value, get_value_range, get_value_choices, set_value,
# FCS
FCS,
# Engine
Engine, get_engine_position, get_engine_orientation, get_engine_gyro_effects,
# Propulsion
Propulsion, get_propulsion_position, get_fuel_mass_props, calculate_propulsion,
get_engines, get_gyro_effects,
# Aerodynamics
Aerodynamics, aerodynamics_from_wind_total, aerodynamics_from_wind_coeff,
aerodynamics_from_body_total, aerodynamics_from_body_coeff, get_pfm,
get_wind_pfm, get_wind_adim_pfm, get_body_pfm, get_body_adim_pfm,
calculate_aerodynamics,
# Aircraft
Aircraft, get_mass_props, get_pfm, get_aerodynamics, get_propulsion, get_name,
get_wing_area, get_wing_span, get_chord, get_arp, get_empty_mass_props,
get_payload_mass_props, calculate_aircraft

include("attitude.jl")
include("position.jl")
include("state.jl")
include("mass.jl")
include("point_forces_moments.jl")
include("atmosphere.jl")
include("wind.jl")
include("gravity.jl")
include("environment.jl")
include("aero_state.jl")
include("controls.jl")
include("fcs.jl")
include("propulsion.jl")
include("aerodynamics.jl")
include("aircraft.jl")
end


text:
# Helper funtions
get_R(P, G, Δ) = sqrt.(Δ) * G * P * G' * sqrt.(Δ)
get_Q(F, R) = F' * R * F + 1.0
get_K(F, R, Q) = R * F / Q
get_m(m, G, K, e) = G * m + K * e'
get_P(R, K, Q) = R - K * K' * Q
get_S(S, k, e, Q) = S/k + e*e'/Q

function diagnostics(ssm::StateSpaceModel, output)
    μ = output.predicted.μ
    Σ = output.predicted.Σ
    e = output.predicted.e
    u = output.predicted.u
    ME = mean(e; dims = 1)
    MAE = mean(abs.(e); dims = 1)
    MSSE = mean(u.^2; dims = 1)
    LL = mean([logpdf(MvTDist(ssm.ν, μ[t, :], Σ[t, :, :]), y[t, :]) for t in 1:size(ssm.y, 1)])
    return (ME = ME, MAE = MAE, MSSE = MSSE, LogLik = LL)
end


text:
# # ======================
# Example usage
# using CSV
# using Plots
#
# TODO data for test could be potentially proided with artifacts
# file_name = "sts_for_VP_test.csv"
# csv_matrix = CSV.read(file_name)[:,2:end]
# almost_not_csv_matrix = Matrix(csv_matrix)
#
# file_name2 = "spikes.csv"
# spikes = load_csv_file_to_array(file_name2)
#
# file_name3 = "th_chunks.csv"
# th_chunks = load_csv_file_to_array(file_name3)
#
# sts = generate_spike_matrix(spikes; th_chunks=th_chunks)
# VPd = [get_selfdist(s; n_chan=32, cost=60., dt=0.01) for s in sts]
#
#
# get_selfdist(st_inp; n_chan=32, cost=60., dt=0.01)
# spkd(s1, s2, cost)
#

#=
An example from other library that computes VP distance

(source:
https://elephant.readthedocs.io/en/latest/reference/_spike_train_processing.html#module-elephant.spike_train_dissimilarity
)
=#

cost = 1
st_a = [10, 20, 30] ./ 10
st_b = [12, 24, 30] ./ 10

result = spkd(st_a, st_b, cost)

@test Int(ceil(result*1000)) == 600


text:
module AppliAR

import AppliSales: Order
import AppliGeneralLedger: JournalEntry
using Dates: Date, DateTime
using DataFrames
using CSV
using Serialization

export process, retrieve_unpaid_invoices, retrieve_paid_invoices, read_bank_statements, report

export UnpaidInvoice, PaidInvoice, meta, header, body, id
export PaidInvoice, stm
export BankStatement, date, descr, iban, amount
export MetaInvoice, order_id, training_id, date, currency, currency_ratio
export Header, invoice_nbr, name, address, postal_code, city, country, order_ref, name_contact, email_contact
export OpentrainingItem, name_training, date, price_per_student, students, vat_perc

# first, link to the model
include("./domain/Domain.jl"); using .Domain
include("./api/API.jl"); using .API
include("./infrastructure/Infrastructure.jl"); using .Infrastructure

# next, submodule Reporting
include("Reporting.jl"); using .Reporting

# Aging report
"""
	report

Generate an aging report

# Example
```
julia> using Dates

julia> using AppliSales

julia> using AppliAR

julia> orders = AppliSales.process();

julia> AppliAR.process(orders);

julia> stm1 = BankStatement(Date(2020-01-15), "Duck City Chronicals Invoice A1002", "NL93INGB", 2420.0);

julia> stms = [stm1];

julia> AppliAR.process(unpaid_invoices, stms);

julia> r = report()
2-element Array{Any,1}:
 AppliAR.Reporting.Aging("A1001", "Scrooge Investment Bank", 2020-07-14, 1210.0, 0 days)
 AppliAR.Reporting.Aging("A1003", "Donalds Hardware Store", 2020-07-14, 1210.0, 0 days)
```
"""
function report(;path_unpaid="./test_invoicing.txt", path_paid="./test_invoicing_paid.txt")
	#@info(path)
	x = Reporting.aging(path_unpaid, path_paid)
	return x
end

export report

end # module


text:
module TestDoctest

import SplittablesBase
using Documenter: doctest
using Test

@testset "doctest" begin
    if lowercase(get(ENV, "JULIA_PKGEVAL", "false")) == "true"
        @info "Skipping doctests on PkgEval."
        return
    end
    doctest(SplittablesBase)
end

end  # module


text:
export constructtriad
export constructtorus
export constructsphere
export constructcylinder
export constructbox
export constructwhirl
export constructframe
export getplane

"""
    constructtriad(q [; length])

Construct a triad with the given configuration `q` and `length`.
"""
function constructtriad(q::Biquaternion; length::Float64 = 1.0)
    map(x -> x + gettranslation(q), rotate(ℝ³.([[0; 0; 0], [length; 0; 0],
                                                [0; 0; 0], [0; length; 0],
                                                [0; 0; 0], [0; 0; length]]),
                                           getrotation(q)))
end

"""
    constructtorus(q, r, R)

Construct a torus of revolution with the given configuration `q`, the smaller radius `r`
and the bigger radius `R`.
"""
function constructtorus(q::Biquaternion,
                        r::Real,
                        R::Real;
                        segments::Int = 36)
    array = Array{ℝ³,2}(undef, segments, segments)
    for i in 1:segments
        for j in 1:segments
            ϕ = i * 2pi / (segments - 1)
            θ = j * 2pi / (segments - 1)
            x₁ = (R + r * cos(ϕ)) * cos(θ)
            x₂ = (R + r * cos(ϕ)) * sin(θ)
            x₃ = r * sin(ϕ)
            array[i, j] = ℝ³(x₁, x₂, x₃)
        end
    end
    map(x -> x + gettranslation(q), rotate(array, getrotation(q)))
end

"""
    constructsphere(q, radius)

Construct a sphere with the given configuration `q` and `radius`.
"""
function constructsphere(q::Biquaternion,
                         radius::Real;
                         segments::Int = 36)
    array = Array{ℝ³,2}(undef, segments, segments)
    lspace = collect(range(float(-pi), stop = float(pi), length = segments))
    for i in 1:segments
        for j in 1:segments
            ϕ = lspace[i]
            θ = lspace[j] / 2
            array[i, j] = ℝ³(Cartesian(Geographic(ϕ, θ))) * float(radius)
        end
    end
    map(x -> x + gettranslation(q), rotate(array, getrotation(q)))
end

"""
    constructcylinder(q, height, radius)

Construct a cylinder with the given configuration `q`, `height` and `radius`.
"""
function constructcylinder(q::Biquaternion,
                           height::Float64,
                           radius::Float64;
                           segments::Int = 36)
    array = Array{ℝ³,2}(undef, segments, segments)
    lspace = collect(range(float(-pi), stop = float(pi), length = segments))
    for i in 1:segments
        for j in 1:segments
            α = lspace[j]
            x₁ = radius * cos(α)
            x₂ = radius * sin(α)
            x₃ = (i - 1) * height / (segments - 1)
            array[i, j] = ℝ³(x₁, x₂, x₃)
        end
    end
    map(x -> x + gettranslation(q), rotate(array, getrotation(q)))
end

"""
    constructbox(p, q, s)

Construct a box with the given translation `p`, rotation `q` and scale `s`.
"""
function constructbox(p::ℝ³, q::Quaternion, s::ℝ³)
    x, y, z = vec(s) ./ 2
    points = [ℝ³(-x, -y, -z) ℝ³(x, -y, -z) ℝ³(-x, -y, -z) ℝ³(x, -y, -z);
              ℝ³(-x, -y, +z) ℝ³(x, -y, +z) ℝ³(-x, -y, +z) ℝ³(x, -y, +z);
              ℝ³(-x, +y, +z) ℝ³(x, +y, +z) ℝ³(-x, +y, +z) ℝ³(x, +y, +z);
              ℝ³(-x, +y, -z) ℝ³(x, +y, -z) ℝ³(-x, +y, -z) ℝ³(x, +y, -z)]
    map(x -> x + p, rotate(points, q))
end

"""
    constructwhirl(points, [top, [bottom, [s3rotation, [config, [segments]]]]])

Construct a whirl with the given `points` in the base space, `top`, `bottom`, S³
rotation `s3rotation`, configuration `config` and the number of `segments`.
"""
function constructwhirl(points::Array{<:S²,1};
                        top::S¹ = U1(-pi),
                        bottom::S¹ = U1(pi),
                        s3rotation::S³ = Quaternion(1, 0, 0, 0),
                        config::Biquaternion = Biquaternion(ℝ³(0, 0, 0)),
                        segments::Int = 36)
    # Map from S² into its upper hemisphere
    f(b::S²) = begin
        p = Geographic(b)
        r = sqrt((1 - sin(p.θ)) / 2)
        Geographic(r * cos(p.ϕ), r * sin(p.ϕ))
    end
    # Map from S² into S³
    function σ(b::S², f::S¹)
        p = Geographic(b)
        arg = angle(exp(-im * p.ϕ) * sqrt(1 + sin(p.θ)))
        z = exp(im * angle(f)) * exp(-im * (arg + p.ϕ)) * sqrt(1 + sin(p.θ)) / sqrt(2)
        w = exp(im * angle(f)) * exp(-im * arg) * sqrt(1 - sin(p.θ)) / sqrt(2)
        Quaternion(ComplexPlane(z, w)) # σ ⊂ [|z| ≤ 1, w > 0]
    end
    # Map from S² into S³
    function π(b::S², f::S¹)
        p = Geographic(b)
        arg = angle(exp(-im * p.ϕ) * sqrt(1 + sin(p.θ)))
        z = exp(im * angle(f)) * exp(-im * arg) * sqrt(1 - sin(p.θ)) / sqrt(2)
        w = exp(im * angle(f)) * exp(-im * (arg + p.ϕ)) * sqrt(1 + sin(p.θ)) / sqrt(2)
        Quaternion(ComplexPlane(z, w)) # π ⊂ [z > 0, |w| ≤ 1]
    end
    array = Array{ℝ³,2}(undef, segments, length(points))
    lspace = range(angle(top), stop = angle(bottom), length = segments)
    for (i, α) in enumerate(lspace)
        for (j, point) in enumerate(points)
            array[i, j] = σmap(rotate(s3rotation, σ(f(point), U1(α))))
        end
    end
    applyconfig(array, config)
end

"""
    constructframe(circle, [s3rotation, [config, [segments]]])

Construct a frame that is like an S³ cross-section under stereographic projection with the
given `circle` that determines the phase angle in the fiber space, S³ rotation `s3rotation`,
configuration `config` and the number of `segments`. The Clifford bundle does not have a
cross section and so we restrict the base space coordinate chart to its upper hemisphere.
"""
function constructframe(circle::S¹;
                        s3rotation::S³ = Quaternion(1, 0, 0, 0),
                        config::Biquaternion = Biquaternion(ℝ³(0, 0, 0)),
                        segments::Int = 36)
    # Map from S² into its upper hemisphere
    f(b::S²) = begin
        p = Geographic(b)
        r = sqrt((1 - sin(p.θ)) / 2)
        Geographic(r * cos(p.ϕ), r * sin(p.ϕ))
    end
    # Map from S² into S³
    function σ(b::S², f::S¹)
        p = Geographic(b)
        arg = angle(exp(-im * p.ϕ) * sqrt(1 + sin(p.θ)))
        z = exp(im * angle(f)) * exp(-im * (arg + p.ϕ)) * sqrt(1 + sin(p.θ)) / sqrt(2)
        w = exp(im * angle(f)) * exp(-im * arg) * sqrt(1 - sin(p.θ)) / sqrt(2)
        Quaternion(ComplexPlane(z, w)) # σ ⊂ [|z| ≤ 1, w > 0]
    end
    # Map from S² into S³
    function π(b::S², f::S¹)
        p = Geographic(b)
        arg = angle(exp(-im * p.ϕ) * sqrt(1 + sin(p.θ)))
        z = exp(im * angle(f)) * exp(-im * arg) * sqrt(1 - sin(p.θ)) / sqrt(2)
        w = exp(im * angle(f)) * exp(-im * (arg + p.ϕ)) * sqrt(1 + sin(p.θ)) / sqrt(2)
        Quaternion(ComplexPlane(z, w)) # π ⊂ [z > 0, |w| ≤ 1]
    end
    array = Array{ℝ³,2}(undef, segments, segments)
    lspaceϕ = collect(range(-pi, stop = pi, length = segments))
    lspaceθ = collect(range(-pi / 2, stop = pi / 2, length = segments))
    for (i, θ) in enumerate(lspaceθ)
        for (j, ϕ) in enumerate(lspaceϕ)
            array[i, j] = σmap(rotate(s3rotation, σ(f(Geographic(ϕ, θ)), circle)))
        end
    end
    applyconfig(array, config)
end

getplane(p::Array{Float64}, q::Quaternion, s::Array{Float64}) = begin
    plane = Array{Float64,3}(undef, 2, 2, 3)
    plane[1, 1, :] = vec(rotate(ℝ³([0.0; -1.0; 1.0]), q)) .* s + p
    plane[1, 2, :] = vec(rotate(ℝ³([0.0; 1.0; 1.0]), q)) .* s + p
    plane[2, 1, :] = vec(rotate(ℝ³([0.0; -1.0; -1.0]), q)) .* s + p
    plane[2, 2, :] = vec(rotate(ℝ³([0.0; 1.0; -1.0]), q)) .* s + p
    plane
end


text:
using Test
using DynamicPolynomials
using SetProg, SetProg.Sets

@testset "Sets" begin
    @testset "ConvexPolySet" begin
        @polyvar x y
        P = SetProg.SumOfSquares.SymMatrix(Float64[1, 2, 3], 2)
        Q = SetProg.SumOfSquares.SymMatrix(BigInt[2, 3, 4], 2)
        basis = SetProg.Sets.MonoBasis(monovec([x, y]))
        q = SetProg.Sets.ConvexPolySet(2, SetProg.SumOfSquares.GramMatrix(P, basis), Q)
        @test q isa SetProg.Sets.ConvexPolySet{BigFloat}
    end
    @testset "zero_eliminate" begin
        @polyvar x y z
        p = SetProg.GramMatrix{Float64}((i, j) -> convert(Float64, i + j),
                                           monovec([x, y, z]))
        set = Sets.ConvexPolySet(2, p, nothing)
        el = Sets.zero_eliminate(set, 1:2)
        @test el.p.Q == 6ones(1, 1)
        el = Sets.zero_eliminate(set, 3:3)
        @test el.p.Q == [2 3; 3 4]
        el = Sets.zero_eliminate(set, 2:2)
        @test el.p.Q == [2 4; 4 6]
    end
end


text:
### special-functions.jl
#
# Copyright (C) 2019 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: uncertainty, error propagation, physics
#
# This file is a part of Measurements.jl.
#
# License is MIT "Expat".
#
### Commentary:
#
# This file contains definition of special functions that support
# Measurement objects.
#
### Code:

### Special functions
using .SpecialFunctions

# Error function: erf, erfinv, erfc, erfcinv, erfcx, erfi, dawson

function SpecialFunctions.erf(a::Measurement{T}) where {T<:AbstractFloat}
    aval = a.val
    return result(erf(aval), 2*exp(-abs2(aval))/sqrt(T(pi)), a)
end

function SpecialFunctions.erfinv(a::Measurement{T}) where {T<:AbstractFloat}
    res = erfinv(a.val)
    # For the derivative, see http://mathworld.wolfram.com/InverseErf.html
    return result(res, sqrt(T(pi)) * exp(abs2(res)) / 2, a)
end

function SpecialFunctions.erfc(a::Measurement{T}) where {T<:AbstractFloat}
    aval = a.val
    return result(erfc(aval), -2*exp(-abs2(aval))/sqrt(T(pi)), a)
end

function SpecialFunctions.erfcinv(a::Measurement{T}) where {T<:AbstractFloat}
    res = erfcinv(a.val)
    # For the derivative, see http://mathworld.wolfram.com/InverseErfc.html
    return result(res, -sqrt(T(pi)) * exp(abs2(res)) / 2, a)
end

function SpecialFunctions.erfcx(a::Measurement{T}) where {T<:AbstractFloat}
    aval = a.val
    res = erfcx(aval)
    return result(res, 2 * (aval * res - inv(sqrt(T(pi)))), a)
end

function SpecialFunctions.erfi(a::Measurement{T}) where {T<:AbstractFloat}
    aval = a.val
    return result(erfi(aval), 2*exp(abs2(aval))/sqrt(T(pi)), a)
end

function SpecialFunctions.dawson(a::Measurement{T}) where {T<:AbstractFloat}
    aval = a.val
    res = dawson(aval)
    return result(res, one(T) - 2 * aval * res, a)
end

# Factorial and gamma

function SpecialFunctions.factorial(a::Measurement)
    aval = a.val
    fact = SpecialFunctions.factorial(aval)
    return result(fact, fact*digamma(aval + one(aval)), a)
end

function SpecialFunctions.gamma(a::Measurement)
    aval = a.val
    Γ = gamma(aval)
    return result(Γ, Γ*digamma(aval), a)
end

function SpecialFunctions.lgamma(a::Measurement)
    aval = a.val
    return result(lgamma(aval), digamma(aval), a)
end

function SpecialFunctions.digamma(a::Measurement)
    aval = a.val
    return result(digamma(aval), trigamma(aval), a)
end

function SpecialFunctions.invdigamma(a::Measurement)
    aval = a.val
    res = invdigamma(aval)
    return result(res, inv(trigamma(res)), a)
end

function SpecialFunctions.trigamma(a::Measurement)
    aval = a.val
    return result(trigamma(aval), polygamma(2, aval), a)
end

function SpecialFunctions.polygamma(n::Integer, a::Measurement)
    aval = a.val
    return result(polygamma(n, aval), polygamma(n + 1, aval), a)
end

# Beta function: beta, lbeta

function SpecialFunctions.beta(a::Measurement, b::Measurement)
    aval = a.val
    bval = b.val
    res = beta(aval, bval)
    return result(res,
                  (res*(digamma(aval) - digamma(aval + bval)),
                   res*(digamma(bval) - digamma(aval + bval))),
                  (a, b))
end

function SpecialFunctions.beta(a::Measurement, b::Real)
    aval = a.val
    res = beta(aval, b)
    return result(res, res*(digamma(aval) - digamma(aval + b)), a)
end

SpecialFunctions.beta(a::Real, b::Measurement) = beta(b, a)

function SpecialFunctions.lbeta(a::Measurement, b::Measurement)
    aval = a.val
    bval = b.val
    return result(lbeta(aval, bval),
                  (digamma(aval) - digamma(aval + bval),
                   digamma(bval) - digamma(aval + bval)),
                  (a, b))
end

function SpecialFunctions.lbeta(a::Measurement, b::Real)
    aval = a.val
    return result(lbeta(aval, b), digamma(aval) - digamma(aval + b), a)
end

SpecialFunctions.lbeta(a::Real, b::Measurement) = lbeta(b, a)

# Airy functions

for f in (:airyai, :airybi)
    der = Symbol(string(f) * "prime")
    @eval begin
        function SpecialFunctions.$f(a::Measurement)
            x = a.val
            return result($f(x), $der(x), a)
        end
        function SpecialFunctions.$der(a::Measurement)
            x = a.val
            # Use Airy equation: y'' - xy = 0 => y'' = xy
            return result($der(x), x * $f(x), a)
        end
    end
end

# Bessel functions

function SpecialFunctions.besselj0(a::Measurement)
    x = a.val
    return result(besselj0(x), -besselj1(x), a)
end

function SpecialFunctions.besselj1(a::Measurement)
    x = a.val
    return result(besselj1(x), (besselj0(x) - besselj(2, x)) / 2, a)
end

# XXX: I don't know a closed form expression for the derivative with respect to
# first argument of J_n.  Arguably, there will be more cases where the
# measurement is the second argument, than the first one.  In any case, you can
# use "@uncertain" macro when both arguments are of Measurement type.
function SpecialFunctions.besselj(nu::Real, a::Measurement)
    x = a.val
    return result(besselj(nu, x), (besselj(nu - 1, x) - besselj(nu + 1, x)) / 2, a)
end

function SpecialFunctions.bessely0(a::Measurement)
    x = a.val
    return result(bessely0(x), -bessely1(x), a)
end

function SpecialFunctions.bessely1(a::Measurement)
    x = a.val
    return result(bessely1(x), (bessely0(x) - bessely(2, x)) / 2, a)
end

# XXX: I don't know a closed form expression for the derivative with respect to
# first argument of y_n, see comments about "besselj".
function SpecialFunctions.bessely(nu::Real, a::Measurement)
    x = a.val
    return result(bessely(nu, x), (bessely(nu - 1, x) - bessely(nu + 1, x)) / 2, a)
end

function SpecialFunctions.besselh(nu::Real, k::Integer, a::Measurement)
    x = a.val
    return result(besselh(nu, k, x),
                  (besselh(nu - 1, k, x) - besselh(nu + 1, k, x)) / 2,
                  a)
end

function SpecialFunctions.besseli(nu::Real, a::Measurement)
    x = a.val
    return result(besseli(nu, x), (besseli(nu - 1, x) + besseli(nu + 1, x)) / 2, a)
end

function SpecialFunctions.besselix(nu::Real, a::Measurement)
    x = a.val
    return result(besselix(nu, x),
                  (besseli(nu - 1, x) + besseli(nu + 1, x)) * exp(-abs(x)) / 2 -
                  besseli(nu, x)*sign(x)*exp(-abs(x)),
                  a)
end

function SpecialFunctions.besselk(nu::Real, a::Measurement)
    x = a.val
    return result(besselk(nu, x), -(besselk(nu - 1, x) + besselk(nu + 1, x)) / 2, a)
end

function SpecialFunctions.besselkx(nu::Real, a::Measurement)
    x = a.val
    return result(besselkx(nu, x),
                  -(besselk(nu - 1, x) + besselk(nu + 1, x)) * exp(x) / 2 +
                  besselk(nu, x)*exp(x),
                  a)
end


text:
using Caesar
using RoMEPlotting, KernelDensityEstimatePlotting
using KernelDensityEstimate, IncrementalInference
using DocStringExtensions, DelimitedFiles
using Gadfly, Cairo, Fontconfig

include(joinpath(@__DIR__,"slamUtils.jl"))

wstart = 415;
wlen = 9;

topdir = joinpath(ENV["HOME"],"data","kayaks")
trialstr = "20_gps_pos";
datawindow = collect(wstart:wstart+wlen-1);
datadir = joinpath(topdir,trialstr);
rangedir = joinpath(topdir, "rangeOnly_"*trialstr)

mfin = zeros(7501,2,wlen);
posData = zeros(wlen,2);
for i in 1:wlen
    dataindex = datawindow[i];
    navfile = datadir*"/nav$(dataindex).csv"
    posData[i,:] = readdlm(navfile,',',Float64,'\n')
    rangeFile = rangedir*"/range$(dataindex).txt";
    mfin[:,:,i] = readdlm(rangeFile,',',Float64,'\n')
end

dposData = deepcopy(posData)
cumulativeDrift!(dposData,[0.0;0],[0.2,0.2])

fg = initfg();

poses = [Symbol("x$i") for i in 1:wlen]
for sym in poses
  addVariable!(fg, sym, Point2)
end

priors = [1,8]
rtkCov = Matrix(Diagonal([0.1;0.1].^2))
#Priors
for i in priors
    pp = PriorPoint2(MvNormal(posData[i,:], rtkCov))
    addFactor!(fg, [Symbol("x$i");], pp, autoinit=true)
end

vps = 1:wlen-1
dpσ = Matrix(Diagonal([0.2;0.2].^2))
for i in vps
    dx = dposData[i+1,1] - dposData[i,1];
    dy = dposData[i+1,2] - dposData[i,2];
    p2p2 = Point2Point2(MvNormal([dx;dy],dpσ))
    addFactor!(fg, [Symbol("x$i");Symbol("x$(i+1)")], p2p2, autoinit=true)
end

beacon = :l1
addVariable!(fg, beacon, Point2 )
# manualinit!(fg,beacon,kde!(rand(MvNormal([0;0],Matrix(Diagonal([7.0;7].^2))),100)))

rangewindow = 1:3:wlen
for i in rangewindow
    sym = Symbol("x$i")
    ppR = Point2Point2Range(AliasingScalarSampler(mfin[:,1,i],exp.(mfin[:,2,i]),SNRfloor=0.8))
    addFactor!(fg, [sym;beacon], ppR, autoinit=false)
end

writeGraphPdf(fg, engine="neato")
wipeBuildNewTree!(fg, drawpdf=true, show=true, imgs=false)
# getSolverParams(fg).drawtree = true
#getSolverParams(fg).showtree = true

## solve the factor graph , recordcliqs=[:x1; :l1; :x8]
tree, smt, hist = solveTree!(fg)

# using FunctionalStateMachine
# using IncrementalInference
# csmAnimate(fg,tree,[:x1; :l1; :x8])

# Debugging here

bss = AliasingScalarSampler(mfin[:,1,4],exp.(mfin[:,2,4]),SNRfloor=0.8)
out = rand(bss,1000)
Gadfly.plot(x=out,Geom.histogram,Coord.cartesian(xmin=0, xmax=100, ymin=0, ymax=400,fixed=true))

pdict = getCliqParentMsgDown(tree,tree.cliques[5])
mykde = pdict[:x4][1];

xx = rand(mykde,1000)
plt = Gadfly.plot(layer(x=xx,Geom.histogram(density=true)),Guide.xlabel("X,Y(m)"));  plt |> PDF("/tmp/plt.pdf");

plt = plotKDEContour(mykde,xlbl="X (m)", ylbl="Y (m)",levels=4); plt |> PDF("/tmp/plt.pdf")

getFrontals(tree.cliques[4])
cdict = getCliqChildMsgsUp(tree,tree.cliques[4],BallTreeDensity)

mykde = cdict[:x4][1][1];
# X1 = rand(mykde,400)
# plk = Gadfly.plot(layer(x=X1[1,:],y=X1[2,:],Geom.histogram2d(xbincount=400, ybincount=400)),Theme(key_position = :none));  plk |> PDF("/tmp/plt.pdf");

# plt = plotKDE(marginal(mykde,[1;]),axis=[0.0 100],N=100); plt |> PDF("/tmp/plt.pdf");

xx = rand(marginal(mykde,[2]),1000)
plt = Gadfly.plot(layer(x=xx,Geom.histogram(density=true)));  plt |> PDF("/tmp/plt.pdf");

plt = plotKDE(mykde,levels=4); plt |> "/tmp/plt.pdf"

plk= [];

for sym in poses #plotting all syms labeled
    X1 = getKDEMean(getVertKDE(fg,sym))
    push!(plk, layer(x=[X1[1];],y=[X1[2];], label=["$(sym)";], Geom.point, Geom.label, Theme(default_color=colorant"blue",point_size = 1.5pt,highlight_width = 0pt)))
    K1 = plotKDEContour(getVertKDE(fg,sym),xlbl="", ylbl="",levels=2,layers=true);
    push!(plk,K1...)
    push!(plk,Gadfly.Theme(key_position = :none));
end
push!(plk, Coord.cartesian(xmin=20, xmax=60, ymin=-60, ymax=-30,fixed=true))
plkplot = Gadfly.plot(plk...); plkplot |> PDF("/tmp/test.pdf");

#
# igt = [17.0499;1.7832];
#
# push!(plk,layer(x=[igt[1];],y=[igt[2];], label=String["Beacon Ground Truth";],Geom.point,Geom.label(hide_overlaps=false), order=2, Theme(default_color=colorant"red",highlight_width = 0pt)));
#
# X1 = approxConv(fg,:x4l1f1,:l1,N=200)
# push!(plk,layer(x=X1[1,:],y=X1[2,:],Geom.histogram2d(xbincount=400, ybincount=400)))
# K1 = plotKDEContour(kde!(X1),xlbl="X (m)", ylbl="Y (m)",levels=4,layers=true);
# push!(plk,K1...)
# push!(plk,Gadfly.Theme(key_position = :none));

# L1 = getVal(getVariable(fg, beacon))
# L1 = rand(getVertKDE(fg,:l1),1000)
# push!(plk,layer(x=L1[1,:],y=L1[2,:],Geom.histogram2d(xbincount=300, ybincount=300)))
# K1 = plotKDEContour(getVertKDE(fg,:l1),xlbl="X (m)", ylbl="Y (m)",levels=6,layers=true);
# push!(plk,K1...)
# push!(plk,Gadfly.Theme(key_position = :none));
# push!(plk, Guide.xlabel("X (m)"),Guide.ylabel("Y (m)"))
# push!(plk, Coord.cartesian(xmin=35, xmax=45, ymin=-50, ymax=-35,fixed=true))

Profile.init(n = 10^7, delay = 0.01)
@profile stuff = sandboxCliqResolveStep(tree,:x3,6)
ProfileView.view()
Profile.print(format=:flat)

## Plot the SAS factor pairs
fsym = :l1x1x2x3x4x5f1
pl = plotSASPair(fg, fsym, show=true, filepath="/tmp/testDP_2pp_3vpvp_init.pdf");

plk=[]
X1 = getKDEMean(getVertKDE(fg,:x1))
push!(plk, layer(x=X1[1,:],y=X1[2,:], Geom.point))
tplt = Gadfly.plot(plk...); tplt |> SVG("/tmp/test.svg") ; tplt |> PDF("/tmp/test.pdf")

plotKDE(fg, ls(fg,r"x"), dims=[1;2])

## More debugging below

stuff = IncrementalInference.localproduct(fg, :x1)
pl = plotKDE(stuff[1], dims=[1;2], levels=3, c=["blue"])

## plotting for dehann
# plk= [];
#
# for sym in poses #plotting all syms labeled
#     X1 = getKDEMean(getVertKDE(fg,sym))
#     push!(plk, layer(x=[X1[1];],y=[X1[2];], label=["$(sym)";], Geom.point, Geom.label, Theme(default_color=colorant"blue",point_size = 1.5pt,highlight_width = 0pt)))
#     K1 = plotKDEContour(getVertKDE(fg,sym),xlbl="", ylbl="",levels=2,layers=true);
#     push!(plk,K1...)
#     push!(plk,Gadfly.Theme(key_position = :none));
# end
# push!(plk, Coord.cartesian(xmin=20, xmax=50, ymin=-60, ymax=-30,fixed=true))
# plkplot = Gadfly.plot(plk...)


text:
"""

Functions that wrap the coordinates

They modify the coordinates of the input vector.

```
function wrap!(x::AbstractVector{T}, sides::T, center::T) where T <: AbstractVector
```

Wrap to a given center of coordinates

"""
function wrap!(x::AbstractVector{T}, sides::T, center::T) where {T<:AbstractVector}
    for i in eachindex(x)
        x[i] = wrapone(x[i], sides, center)
    end
    nothing
end

@inline function wrapone(x::T, sides::T, center::T) where {T<:AbstractVector}
    s = @. (x - center) % sides
    s = @. wrapx(s, sides) + center
    return s
end

@inline function wrapx(x, s)
    if x > s / 2
        x = x - s
    elseif x < -s / 2
        x = x + s
    end
    x
end

"""

```
wrap!(x::AbstractVector{T}, sides::T) where T <: AbstractVector
```

Wrap to origin.

"""
function wrap!(x::AbstractVector{T}, sides::T) where {T<:AbstractVector}
    for i in eachindex(x)
        x[i] = wrapone(x[i], sides)
    end
    nothing
end

@inline function wrapone(x::T, sides::T) where {T<:AbstractVector}
    s = @. x % sides
    s = @. wrapx(s, sides)
    return s
end


text:
using Distributions
using LightXML
using Random

global COLLISION_COUNT = 0
global PORT = "2222" # CARLA TCP port
if !@isdefined(CARLA_PROCESS)
    global CARLA_PROCESS = nothing
end

py_ver = read(`python --version`, String)
if py_ver[1:10] != "Python 3.7"
    # run(`../set_python37.bat`)
    error("Please set Python to 3.7")
end

## 1) Run this first. Wait.
function open_carla()
    global PORT, CARLA_PROCESS
    @info "Opening CARLA on port $PORT."

    use_config_py = false

    kill_carla()
    @async run(`..\\..\\CarlaUE4.exe -carla-rpc-port=$PORT`)
    # CARLA_PROCESS = run(`..\\..\\CarlaUE4.exe -carla-rpc-port=$PORT`, wait=false)

    if use_config_py
        # preload the correct map (Town01 for red light scenarios)
        @async run(`python ..\\util\\config.py --port=$PORT --map Town01`)
    end
end

"""
Kill the CARLA process (on Windows x64).
"""
function kill_carla()
    global CARLA_PROCESS
    # if !isnothing(CARLA_PROCESS)
    #     @info "Killing $CARLA_PROCESS"
    #     kill(CARLA_PROCESS)
    # end
    try
        run(`taskkill /f /t /im CarlaUE4-Win64-Shipping.exe`)
    catch err end
end

function add_disturbance(xdoc)
    x = parse(Float64, attribute(root(xdoc)["ego_vehicle"][1], "x"))
    noise_model = Uniform(0, 20)
    ϵ = rand(noise_model)
    x′ = string(x + ϵ)
    @info x′
    set_attribute(root(xdoc)["ego_vehicle"][1], "x", x′)
end

function disturb_file()
    @info "Adding disturbance to initial positions."
    xdoc = parse_file("red_light_config.xml")
    add_disturbance(xdoc)
    disturb_filename = "red_light_config_2.xml"
    save_file(xdoc, disturb_filename)
    return disturb_filename
end

## 2) Then run this. Wait.
function start_scenario_runner(; reload_world=true, tm_port=8001, xml="red_light_config.xml")
    global PORT
    @info "Starting scenario runner."

    # NOTE, fails to run iteratively in "--sync" mode.
    # arguments = ["--scenario", "OppositeVehicleRunningRedLight001", "--configFile", xml, "--trafficManagerPort", "$tm_port", "--timeout", "120", "--port", PORT, "--sync"]

    # Running full set of scenarios
    # arguments = ["--scenario", "group:OppositeVehicleRunningRedLight", "--configFile", "RunningRedLightExample.xml", "--trafficManagerPort", "$tm_port", "--timeout", "120", "--port", PORT]

    arguments = ["--scenario", "OppositeVehicleRunningRedLight001", "--configFile", xml, "--trafficManagerPort", "$tm_port", "--timeout", "120", "--port", PORT]
    if reload_world
        push!(arguments, "--reloadWorld")
    end

    scenerio_runner_cmd = ["python", "scenario_runner.py", arguments...]
    @info scenerio_runner_cmd
    @async run(Cmd(scenerio_runner_cmd))
end

## 3) Then run this.
# @async run(`python ast_no_rendering_mode.py`)
ENV["PYTHON"] = "C:\\Python37\\python.exe"
using PyCall

# Run AST agent without rendering
function run_no_render(rerun=false)
    global PORT
    @info "Running AST agent in 'no render' mode."
    # @async run(`python ast_no_rendering_mode.py --port $PORT`)

    ## NOTE: Comment out "if __name__" in ast_no_rendering_mode.py when using his py"" call style.
    # Async Python call.
    if rerun
        py"""
        rerun_game_loop()
        """
        # py"main($PORT)"
    else
        # Load Python code into enviroment
        py"""
        exec(open('ast_no_rendering_mode.py').read())
        """

        py"""
        main($PORT)
        """
    end
end

function runcarla(newcarla=false; N=10, reload_world=true, tm_port=8001, rerun=false)
    global CARLA_PROCESS

    # Traffic Manager Port: increase every run? (https://github.com/carla-simulator/carla/issues/2789#issuecomment-760951212)
    if newcarla
        # CARLA_PROCESS = nothing
        @warn "Please run `open_carla()` separately, then `runcarla(false)`"
        open_carla()
        sleep(10)
    end

    for n in 1:N
        @info "Running simulation $n"
        Random.seed!(n)

        # 2. scenario runner
        disturb_filename = disturb_file()
        # start_scenario_runner(tm_port=tm_port, xml=disturb_filename, reload_world=reload_world)
        start_scenario_runner(tm_port=tm_port, xml=disturb_filename, reload_world=reload_world)
        # start_scenario_runner(tm_port=tm_port, xml="red_light_config_no_collision.xml", reload_world=reload_world)
        sleep(10) # wait to run "no render" until this finishes (TODO: exact.)

        # 3. AST agent ("no render" mode)
        run_no_render(rerun || n > 1) # TODO. speed up?

        # println("Going to kill...")
        # sleep(20)
        # println("Killing...")
        # kill(process)
        num_collisions = py"NUM_COLLISIONS"
        @info num_collisions, num_collisions / n
    end
end

# function destroy()

##################################################
# Main: Open CARLA before running anything.
##################################################
open_carla()


text:
############################################################################################
# Types and structs - Reweighting Methods
struct Bootstrap <: BaytesCore.ParameterWeighting end

############################################################################################
# Bootstrap Method
function weight!(
    method::Bootstrap,
    weights::BaytesCore.ParameterWeights,
    dataₜ::D,
    kernel::ParticleKernel,
    val::AbstractMatrix{P},
    iter::Integer,
) where {D,P}
    ##Calculate weights
    ℓevidence!(weights.ℓweights, dataₜ, kernel, val, iter)
    ##Normalize weights
    BaytesCore.normalize!(weights)
    return nothing
end

############################################################################################
#export
export Bootstrap, weight!


text:
function check_manifold(M)
    # check random points, vectors
    for i in 1:10
        x = randomMPoint(M)
        @test is_manifold_point(M, x)

        ξ = randomTVector(M, x)
        @test is_tangent_vector(M, x, ξ, atol = 1e-10)
        @test isapprox(norm(M, x, ξ), 1)

        # check projection
        q = rand(representation_size(M)...)

        η = zero_tangent_vector(M, x)
        project!(M, η, x, q)
        @test is_tangent_vector(M, x, η; atol = 1e-10)

        # check retraction
        y = retract(M, x, ξ)
        @test is_manifold_point(M, y)
    end

    return nothing
end


text:
# ---
# title: 1567. Maximum Length of Subarray With Positive Product
# id: problem1567
# author: zhwang
# date: 2022-01-23
# difficulty: Medium
# categories: Greedy
# link: <https://leetcode.com/problems/maximum-length-of-subarray-with-positive-product/description/>
# hidden: true
# ---
#
# Given an array of integers `nums`, find the maximum length of a subarray where
# the product of all its elements is positive.
#
# A subarray of an array is a consecutive sequence of zero or more values taken
# out of that array.
#
# Return  _the maximum length of a subarray with positive product_.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,-2,-3,4]
#     Output: 4
#     Explanation: The array nums already has a positive product of 24.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [0,1,-2,-3,-4]
#     Output: 3
#     Explanation: The longest subarray with positive product is [1,-2,-3] which has a product of 6.
#     Notice that we cannot include 0 in the subarray since that'll make the product 0 which is not positive.
#
# **Example 3:**
#
#
#
#     Input: nums = [-1,-2,-3,0,1]
#     Output: 2
#     Explanation: The longest subarray with positive product is [-1,-2] or [-2,-3].
#
#
# **Example 4:**
#
#
#
#     Input: nums = [-1,2]
#     Output: 1
#
#
# **Example 5:**
#
#
#
#     Input: nums = [1,2,3,5,-6,4,0,10]
#     Output: 4
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 10^5`
#   * `-10^9 <= nums[i] <= 10^9`
#
#
## @lc code=start
using LeetCode

function get_max_len(nums::Vector{Int})::Int
    pos_pre = neg_pre = max_ans = 0
    for num in nums
        if num == 0
            pos_pre = neg_pre = 0
        elseif num > 0
            pos_pre = pos_pre + 1
            neg_pre = neg_pre != 0 ? neg_pre + 1 : 0
        else
            neg_tmp = neg_pre
            neg_pre = pos_pre != 0 ? pos_pre + 1 : 1
            pos_pre = neg_tmp != 0 ? neg_tmp + 1 : 0
        end
        max_ans = max(max_ans, pos_pre)
    end
    return max_ans
end

## @lc code=end


text:
"""
    Call(args...; kwargs...)

Represents a function call.
"""
struct Call{T<:Tuple, P<:Pairs}
    args::T
    kwargs::P

    Call(args...; kwargs...) = new{typeof(args), typeof(kwargs)}(args, kwargs)
end

"""
    Predicate(f)

A predicate that can be used instead of literal call arguments,
useful for situations where the exact value of arguments is not known.
The function `f` is called with the observed argument, and must return a `Bool`.

## Example

```julia
m = Mock()
m(rand(1:10))
@assert called_with(m, Predicate(in(1:10)))
```
"""
struct Predicate{F}
    f::F
end

Base.show(io::IO, ::MIME"text/plain", p::Predicate) = print(io, "Predicate(", p.f, ")")

(p::Predicate)(x) = p.f(x)

Base.:(==)(a::Call{T1, P1}, b::Call{T2, P2}) where {T1, P1, T2, P2} = false
Base.:(==)(a::Call{T, P}, b::Call{T, P}) where {T, P} =
    a.args == b.args && a.kwargs == b.kwargs

function Base.show(io::IO, ::MIME"text/plain", c::Call)
    print(io, "Call(")
    print(io, join(map(repr, c.args), ", "))
    isempty(c.kwargs) ||
        print(io, "; ", join(map(((k, v),) -> "$k=$(repr(v))", collect(c.kwargs)), ", "))
    print(io, ")")
end

"""
    Mock([effect])

Create a new mocking object that can act as a replacement for a function.

## Effects

Use the `effect` argument to determine what happens upon calling the mock.

- If the value is callable, then it is called with the same arguments and keywords.
- If the value is an `Exception`, then the exception is thrown.
- If the value is a `Vector`, then each call consumes the next element.
  Nested `Vector`s are returned, but callables and exceptions are treated as explained above.
- Any other value is returned without modification.

By default, calling a `Mock` returns a new `Mock`.
"""
struct Mock{T}
    id::Symbol
    calls::Vector{Call}
    effect::T
end

Mock(effect=(_args...; _kwargs...) -> Mock()) = Mock(gensym(), Call[], effect)

Base.:(==)(a::Mock, b::Mock) = a.id === b.id

Base.show(io::IO, ::MIME"text/plain", m::Mock) = print(io, "Mock(id=$(m.id))")

"""
    (::Mock)(args...; kwargs...)

Calling a `Mock` records the call in its history and triggers its `effect`.
"""
function (m::Mock)(args...; kwargs...)
    push!(calls(m), Call(args...; kwargs...))
    return do_effect(m.effect, args...; kwargs...)
end

"""
    calls(::Mock) -> Vector{Call}

Return the call history of the [`Mock`](@ref).
"""
calls(m::Mock) = m.calls

"""
    last_call(::Mock) -> Call

Return the last call of the [`Mock`](@ref).
"""
last_call(m::Mock) = last(calls(m))

"""
    ncalls(::Mock) -> Int

Return the number of times that the [`Mock`](@ref) has been called.
"""
ncalls(m::Mock) = length(calls(m))

"""
    called(::Mock) -> Bool

Return whether or not the [`Mock`](@ref) has been called.
"""
called(m::Mock) = !isempty(calls(m))

"""
    called_once(::Mock) -> Bool

Return whether or not the [`Mock`](@ref) has been called exactly once.
"""
called_once(m::Mock) = length(calls(m)) == 1

"""
    called_with(::Mock, args...; kwargs...) -> Bool

Return whether or not the [`Mock`](@ref) has been called with the given arguments.
"""
called_with(m::Mock, args...; kwargs...) = has_call(m, Call(args...; kwargs...))

"""
    called_last_with(::Mock, args...; kwargs...) -> Bool

Return whether or not the [`Mock`](@ref) was last called with the given arguments.
"""
called_last_with(m::Mock, args...; kwargs...) = last_call(m) == Call(args...; kwargs...)

"""
    called_once_with(::Mock, args...; kwargs...) -> Bool

Return whether or not the [`Mock`](@ref) has been called exactly once with the given arguments.
"""
called_once_with(m::Mock, args...; kwargs...) =
    called_once(m) && called_with(m, args...; kwargs...)

"""
    has_call(::Mock, ::Call) -> Bool

Similiar to [`called_with`](@ref), but using a [`Call`](@ref).
"""
has_call(m::Mock, c::Call) = any(call -> call_matches(c, call), calls(m))

"""
    has_calls(::Mock, ::Call...) -> Bool

Return whether or not the [`Mock`](@ref) has a particular ordered sequence of [`Call`](@ref)s.
"""
function has_calls(m::Mock, cs::Call...)
    isempty(cs) && return true
    existing = calls(m)
    length(cs) > length(existing) && return false
    cs = collect(cs)  # Omitting this causes a segfault?!
    n = length(cs) - 1
    for i in 1:(length(existing) - n)
        all(ab -> call_matches(ab...), zip(cs, existing[i:i+n])) && return true
    end
    return false
end

"""
    reset!(::Mock)

Reset a [`Mock`](@ref)'s call history.
Its `effect` is preserved.
"""
function reset!(m::Mock)
    empty!(m.calls)
    return m
end

# Handle the effect.
do_effect(x, args...; kwargs...) = isempty(methods(x)) ? x : x(args...; kwargs...)
do_effect(x::Callable, args...; kwargs...) = x(args...; kwargs...)
do_effect(ex::Exception, _args...; _kwargs...) = throw(ex)
function do_effect(xs::Vector, args...; kwargs...)
    x = popfirst!(xs)
    return x isa Vector ? x : do_effect(x, args...; kwargs...)
end

# Match call arguments.
call_matches(expected, observed) =
    length(expected.args) == length(observed.args) &&
    issetequal(keys(expected.kwargs), keys(observed.kwargs)) &&
    all(ab -> arg_matches(ab...), zip(expected.args, observed.args)) &&
    all(k -> arg_matches(expected.kwargs[k], observed.kwargs[k]), keys(expected.kwargs))

arg_matches(expected, observed) = expected == observed
arg_matches(pred::Predicate, observed) = pred(observed)


text:
a = 2
function bar() println(a); a = 3 end
bar()

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia