text:
function doSPHmagic(x, ρ, u, v, p, L, H)

    # 1 - DirectFind algorithm to find interacting particles
    (pairs, W, ∂x_W, W_0) = directfind(x)

    # 2 - SPH density approximation
    ρ = computedensity(pairs, W, W_0)

    # 3 - Equation of state for ideal gases
    p = (u .* ρ) * (GAMMA - 1)

    # 4,5 - Compute L,H operands
    (L, H) = computeoperands(ρ, p, v, pairs, ∂x_W)

    # 6 - Add artificial viscosity contribution
    (L, H) = correctwithartviscosity(x, ρ, p, v, pairs, ∂x_W, L, H)

    (x, ρ, u, v, p, L, H)
end

function directfind(x)
    pairs = []
    W = []
    ∂x_W = []

    for i in 1:PARTICLES_COUNT-1
        for j in i+1:PARTICLES_COUNT
            r = x[i] - x[j]
            if abs(r) < 2 * 2HSML
                # Saving the interacting pair
                push!(pairs, (i, j))

                # Computing and saving the corresponding smoothing kernel value
                (W_k, ∂x_W_k) = doublecosine(r, HSML)
                push!(W, W_k)
                push!(∂x_W, ∂x_W_k)
            end
        end
    end

    (W_0, _) = doublecosine(0, HSML)

    (pairs, W, ∂x_W, W_0)
end

function computedensity(pairs, W, W_0)
    # Considering the self effect
    ρ = ones(PARTICLES_COUNT) * PARTICLE_MASS * W_0

    # Computing for the neighbours
    for k in 1:length(pairs)
        (i, j) = pairs[k]
        ρ[i] += W[k] * PARTICLE_MASS
        ρ[j] += W[k] * PARTICLE_MASS
    end

    ρ
end

function computeoperands(ρ, p, v, pairs, ∂x_W)
    # Zero-initialization, since ∂x_W(0) = 0
    L = zeros(PARTICLES_COUNT)
    H = zeros(PARTICLES_COUNT)

    # Computing with the neighbours
    for k in 1:length(pairs)
        # Getting pairs
        (i, j) = pairs[k]

        # Calculating the constant term
        p_ij = p[i]/ρ[i]^2 + p[j]/ρ[j]^2

        # Wrapping for L and both particles, remembering ∂x_W_ij = -∂x_W_ji
        L[i] += ∂x_W[k] * p_ij * PARTICLE_MASS
        L[j] -= ∂x_W[k] * p_ij * PARTICLE_MASS

        # Wrapping for L and both particles, remembering ∂x_W_ij = -∂x_W_ji
        H[i] -= ∂x_W[k] * 0.5 * (v[i] - v[j]) * p_ij * PARTICLE_MASS
        H[j] += ∂x_W[k] * 0.5 * (v[j] - v[i]) * p_ij * PARTICLE_MASS
    end

    (L, H)
end

function correctwithartviscosity(x, ρ, p, v, pairs, ∂x_W, L, H)
    # Constants
    qa = 1
    qb = 1
    etq = 0.1

    # Initialization of the corrections
    art_L = zeros(PARTICLES_COUNT)
    art_H = zeros(PARTICLES_COUNT)

    # Looping over pairs
    for k in 1:length(pairs)
        # Getting pairs
        (i, j) = pairs[k]

        x_ij = x[i] - x[j]
        v_ij = v[i] - v[j]

        xpsv = v_ij * x_ij
        xpsx = x_ij^2

        if (xpsv < 0)
            # Intermediary results
            muv = HSML * xpsv / (xpsx + HSML^2 * etq^2)
            mc = (sqrt(abs(GAMMA * p[i]/ρ[i])) + sqrt(abs(GAMMA * p[j]/ρ[j]))) / 2
            ρ_m = (ρ[i] + ρ[j]) / 2
            piv = (qb * muv - qa * mc) * muv/ρ_m
            hx = -piv * ∂x_W[k]

            # Adding to artificial acceleration contribution
            art_L[i] += PARTICLE_MASS * hx
            art_L[j] -= PARTICLE_MASS * hx

            # Adding to artificial energy contribution
            art_H[i] += v_ij / 2 * PARTICLE_MASS * hx
            art_H[j] -= -v_ij / 2 * PARTICLE_MASS * hx
        end
    end

    # Applying the contribution
    L -= art_L
    H -= art_H

    (L, H)
end


text:
using LinearAlgebra
using ManifoldsBase
using ManifoldsBase: DefaultManifold, ℝ, ℂ
using Test
import Base: +, -, *, copyto!, isapprox
import ManifoldsBase: allocate

struct ProjManifold <: Manifold{ℝ} end

ManifoldsBase.inner(::ProjManifold, x, w, v) = dot(w, v)
ManifoldsBase.project!(::ProjManifold, w, x, v) = (w .= v .- dot(x, v) .* x)
ManifoldsBase.representation_size(::ProjManifold) = (2, 3)
ManifoldsBase.manifold_dimension(::ProjManifold) = 5
ManifoldsBase.get_vector(::ProjManifold, x, v, ::DefaultOrthonormalBasis) = reverse(v)

@testset "Dispatch" begin
    @test ManifoldsBase.decorator_transparent_dispatch(
        get_basis,
        DefaultManifold(3),
        [0.0, 0.0, 0.0],
        DefaultBasis(),
    ) === Val(:parent)
    @test ManifoldsBase.decorator_transparent_dispatch(
        get_coordinates,
        DefaultManifold(3),
        [0.0, 0.0, 0.0],
    ) === Val(:parent)
    @test ManifoldsBase.decorator_transparent_dispatch(
        get_coordinates!,
        DefaultManifold(3),
        [0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0],
    ) === Val(:transparent)
    @test ManifoldsBase.decorator_transparent_dispatch(
        get_vector,
        DefaultManifold(3),
        [0.0, 0.0, 0.0],
    ) === Val(:parent)
    @test ManifoldsBase.decorator_transparent_dispatch(
        get_vector!,
        DefaultManifold(3),
        [0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0],
    ) === Val(:transparent)
end

struct ProjectionTestManifold <: Manifold{ℝ} end

ManifoldsBase.inner(::ProjectionTestManifold, ::Any, X, Y) = dot(X, Y)
function ManifoldsBase.project!(::ProjectionTestManifold, Y, p, X)
    Y .= X .- dot(p, X) .* p
    Y[end] = 0
    return Y
end
ManifoldsBase.manifold_dimension(::ProjectionTestManifold) = 100

@testset "Projected and arbitrary orthonormal basis" begin
    M = ProjManifold()
    x = [
        sqrt(2)/2 0.0 0.0
        0.0 sqrt(2)/2 0.0
    ]

    for pb in (ProjectedOrthonormalBasis(:svd), ProjectedOrthonormalBasis(:gram_schmidt))
        pb = get_basis(M, x, pb)
        @test number_system(pb) == ℝ
        @test get_basis(M, x, pb) == pb
        N = manifold_dimension(M)
        @test isa(pb, CachedBasis)
        @test CachedBasis(pb) === pb
        @test length(get_vectors(M, x, pb)) == N
        # test orthonormality
        for i in 1:N
            @test norm(M, x, get_vectors(M, x, pb)[i]) ≈ 1
            for j in (i + 1):N
                @test inner(M, x, get_vectors(M, x, pb)[i], get_vectors(M, x, pb)[j]) ≈ 0 atol =
                    1e-15
            end
        end
        # check projection idempotency
        for i in 1:N
            @test project(M, x, get_vectors(M, x, pb)[i]) ≈ get_vectors(M, x, pb)[i]
        end
    end
    aonb = get_basis(M, x, DefaultOrthonormalBasis())
    @test size(get_vectors(M, x, aonb)) == (5,)
    @test get_vectors(M, x, aonb)[1] ≈ [0, 0, 0, 0, 1]

    @testset "Gram-Schmidt" begin
        # for a basis
        M = ManifoldsBase.DefaultManifold(3)
        p = zeros(3)
        V = [[2.0, 0.0, 0.0], [1.1, 2.2, 0.0], [0.0, 3.3, 4.4]]
        b1 = ManifoldsBase.gram_schmidt(M, p, V)
        b2 = ManifoldsBase.gram_schmidt(M, zeros(3), CachedBasis(DefaultBasis(), V))
        @test b1 == get_vectors(M, p, b2)
        # projected gram schmidt
        tm = ProjectionTestManifold()
        bt = ProjectedOrthonormalBasis(:gram_schmidt)
        p = [sqrt(2) / 2, 0.0, sqrt(2) / 2, 0.0, 0.0]
        @test_throws ErrorException get_basis(tm, p, bt)
        b = get_basis(
            tm,
            p,
            bt;
            return_incomplete_set = true,
            warn_linearly_dependent = true,
        )
        @test length(get_vectors(tm, p, b)) == 3
    end
end

struct NonManifold <: Manifold{ℝ} end
struct NonBasis <: ManifoldsBase.AbstractBasis{ℝ} end

struct NonBroadcastBasisThing{T}
    v::T
end

+(a::NonBroadcastBasisThing, b::NonBroadcastBasisThing) = NonBroadcastBasisThing(a.v + b.v)
*(α, a::NonBroadcastBasisThing) = NonBroadcastBasisThing(α * a.v)
-(a::NonBroadcastBasisThing, b::NonBroadcastBasisThing) = NonBroadcastBasisThing(a.v - b.v)

isapprox(a::NonBroadcastBasisThing, b::NonBroadcastBasisThing) = isapprox(a.v, b.v)

function ManifoldsBase.number_eltype(a::NonBroadcastBasisThing)
    return typeof(reduce(+, one(number_eltype(eti)) for eti in a.v))
end

allocate(a::NonBroadcastBasisThing) = NonBroadcastBasisThing(allocate(a.v))
function allocate(a::NonBroadcastBasisThing, ::Type{T}) where {T}
    return NonBroadcastBasisThing(allocate(a.v, T))
end
allocate(::NonBroadcastBasisThing, ::Type{T}, s::Integer) where {S,T} = Vector{T}(undef, s)

function copyto!(a::NonBroadcastBasisThing, b::NonBroadcastBasisThing)
    copyto!(a.v, b.v)
    return a
end

function ManifoldsBase.log!(
    ::DefaultManifold,
    v::NonBroadcastBasisThing,
    x::NonBroadcastBasisThing,
    y::NonBroadcastBasisThing,
)
    return copyto!(v, y - x)
end

function ManifoldsBase.exp!(
    ::DefaultManifold,
    y::NonBroadcastBasisThing,
    x::NonBroadcastBasisThing,
    v::NonBroadcastBasisThing,
)
    return copyto!(y, x + v)
end

function ManifoldsBase.get_basis(
    M::DefaultManifold,
    p::NonBroadcastBasisThing,
    B::DefaultOrthonormalBasis,
)
    return CachedBasis(
        B,
        [
            NonBroadcastBasisThing(ManifoldsBase._euclidean_basis_vector(p.v, i)) for
            i in eachindex(p.v)
        ],
    )
end
function ManifoldsBase.get_basis(
    M::DefaultManifold,
    p::NonBroadcastBasisThing,
    B::DefaultOrthogonalBasis,
)
    return CachedBasis(
        B,
        [
            NonBroadcastBasisThing(ManifoldsBase._euclidean_basis_vector(p.v, i)) for
            i in eachindex(p.v)
        ],
    )
end
function ManifoldsBase.get_basis(
    M::DefaultManifold,
    p::NonBroadcastBasisThing,
    B::DefaultBasis,
)
    return CachedBasis(
        B,
        [
            NonBroadcastBasisThing(ManifoldsBase._euclidean_basis_vector(p.v, i)) for
            i in eachindex(p.v)
        ],
    )
end

function ManifoldsBase.get_coordinates!(
    M::DefaultManifold,
    Y,
    p::NonBroadcastBasisThing,
    X::NonBroadcastBasisThing,
    B::DefaultOrthonormalBasis,
)
    copyto!(Y, reshape(X.v, manifold_dimension(M)))
    return Y
end

function ManifoldsBase.get_vector!(
    M::DefaultManifold,
    Y::NonBroadcastBasisThing,
    p::NonBroadcastBasisThing,
    X,
    B::DefaultOrthonormalBasis,
)
    copyto!(Y.v, reshape(X, representation_size(M)))
    return Y
end

function ManifoldsBase.inner(
    ::DefaultManifold,
    x::NonBroadcastBasisThing,
    v::NonBroadcastBasisThing,
    w::NonBroadcastBasisThing,
)
    return dot(v.v, w.v)
end

ManifoldsBase._get_vector_cache_broadcast(::NonBroadcastBasisThing) = Val(false)

DiagonalizingBasisProxy() = DiagonalizingOrthonormalBasis([1.0, 0.0, 0.0])

@testset "ManifoldsBase.jl stuff" begin

    @testset "Errors" begin
        m = NonManifold()
        onb = DefaultOrthonormalBasis()

        @test_throws ErrorException get_basis(m, [0], onb)
        @test_throws ErrorException get_basis(m, [0], NonBasis())
        @test_throws ErrorException get_coordinates(m, [0], [0], onb)
        @test_throws ErrorException get_coordinates!(m, [0], [0], [0], onb)
        @test_throws ErrorException get_vector(m, [0], [0], onb)
        @test_throws ErrorException get_vector!(m, [0], [0], [0], onb)
        @test_throws ErrorException get_vectors(m, [0], NonBasis())
    end

    M = DefaultManifold(3)

    _pts = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
    @testset "basis representation" for BT in (
            DefaultBasis,
            DefaultOrthonormalBasis,
            DefaultOrthogonalBasis,
            DiagonalizingBasisProxy,
        ),
        pts in (_pts, map(NonBroadcastBasisThing, _pts))

        if BT == DiagonalizingBasisProxy && pts !== _pts
            continue
        end
        v1 = log(M, pts[1], pts[2])
        @test ManifoldsBase.number_of_coordinates(M, BT()) == 3

        if BT != DiagonalizingBasisProxy
            vb = get_coordinates(M, pts[1], v1, BT())
            @test isa(vb, AbstractVector)
            vbi = get_vector(M, pts[1], vb, BT())
            @test isapprox(M, pts[1], v1, vbi)
        end

        b = get_basis(M, pts[1], BT())
        if BT != DiagonalizingBasisProxy
            if pts[1] isa Array
                @test isa(b, CachedBasis{ℝ,BT{ℝ},Vector{Vector{Float64}}})
            else
                @test isa(
                    b,
                    CachedBasis{ℝ,BT{ℝ},Vector{NonBroadcastBasisThing{Vector{Float64}}}},
                )
            end
        end
        @test get_basis(M, pts[1], b) === b
        N = manifold_dimension(M)
        @test length(get_vectors(M, pts[1], b)) == N
        # check orthonormality
        if BT isa DefaultOrthonormalBasis && pts[1] isa Vector
            for i in 1:N
                @test norm(M, pts[1], get_vectors(M, pts[1], b)[i]) ≈ 1
                for j in (i + 1):N
                    @test inner(
                        M,
                        pts[1],
                        get_vectors(M, pts[1], b)[i],
                        get_vectors(M, pts[1], b)[j],
                    ) ≈ 0
                end
            end
            # check that the coefficients correspond to the basis
            for i in 1:N
                @test inner(M, pts[1], v1, get_vectors(M, pts[1], b)[i]) ≈ vb[i]
            end
        end

        if BT != DiagonalizingBasisProxy
            @test get_coordinates(M, pts[1], v1, b) ≈ get_coordinates(M, pts[1], v1, BT())
            @test get_vector(M, pts[1], vb, b) ≈ get_vector(M, pts[1], vb, BT())
        end

        v1c = Vector{Float64}(undef, 3)
        get_coordinates!(M, v1c, pts[1], v1, b)
        @test v1c ≈ get_coordinates(M, pts[1], v1, b)

        v1cv = allocate(v1)
        get_vector!(M, v1cv, pts[1], v1c, b)
        @test isapprox(M, pts[1], v1, v1cv)
    end
end

@testset "Complex DeaultManifold with real and complex Cached Bases" begin
    M = ManifoldsBase.DefaultManifold(3; field = ℂ)
    p = [1.0, 2.0im, 3.0]
    X = [1.2, 2.2im, 2.3im]
    b = [Matrix{Float64}(I, 3, 3)[:, i] for i in 1:3]
    Bℝ = CachedBasis(DefaultOrthonormalBasis{ℝ}(), b)
    aℝ = get_coordinates(M, p, X, Bℝ)
    Yℝ = get_vector(M, p, aℝ, Bℝ)
    @test Yℝ ≈ X
    @test ManifoldsBase.number_of_coordinates(M, Bℝ) == 3

    bℂ = [b..., (b .* 1im)...]
    Bℂ = CachedBasis(DefaultOrthonormalBasis{ℂ}(), bℂ)
    aℂ = get_coordinates(M, p, X, Bℂ)
    Yℂ = get_vector(M, p, aℂ, Bℂ)
    @test Yℂ ≈ X
    @test ManifoldsBase.number_of_coordinates(M, Bℂ) == 6
end

@testset "Basis show methods" begin
    @test sprint(show, DefaultBasis()) == "DefaultBasis(ℝ)"
    @test sprint(show, DefaultOrthogonalBasis()) == "DefaultOrthogonalBasis(ℝ)"
    @test sprint(show, DefaultOrthonormalBasis()) == "DefaultOrthonormalBasis(ℝ)"
    @test sprint(show, DefaultOrthonormalBasis(ℂ)) == "DefaultOrthonormalBasis(ℂ)"
    @test sprint(show, GramSchmidtOrthonormalBasis(ℂ)) == "GramSchmidtOrthonormalBasis(ℂ)"
    @test sprint(show, ProjectedOrthonormalBasis(:svd)) ==
          "ProjectedOrthonormalBasis(:svd, ℝ)"
    @test sprint(show, ProjectedOrthonormalBasis(:gram_schmidt, ℂ)) ==
          "ProjectedOrthonormalBasis(:gram_schmidt, ℂ)"

    @test sprint(show, "text/plain", DiagonalizingOrthonormalBasis(Float64[1, 2, 3])) == """
                                                                               DiagonalizingOrthonormalBasis(ℝ) with eigenvalue 0 in direction:
                                                                               3-element $(sprint(show, Vector{Float64})):
                                                                                1.0
                                                                                2.0
                                                                                3.0"""

    M = DefaultManifold(2, 3)
    x = collect(reshape(1.0:6.0, (2, 3)))
    pb = get_basis(M, x, DefaultOrthonormalBasis())

    @test sprint(show, "text/plain", pb) == """
    DefaultOrthonormalBasis(ℝ) with 6 basis vectors:
     E1 =
      2×3 $(sprint(show, Matrix{Float64})):
       1.0  0.0  0.0
       0.0  0.0  0.0
     E2 =
      2×3 $(sprint(show, Matrix{Float64})):
       0.0  0.0  0.0
       1.0  0.0  0.0
     ⋮
     E5 =
      2×3 $(sprint(show, Matrix{Float64})):
       0.0  0.0  1.0
       0.0  0.0  0.0
     E6 =
      2×3 $(sprint(show, Matrix{Float64})):
       0.0  0.0  0.0
       0.0  0.0  1.0"""

    b = DiagonalizingOrthonormalBasis(get_vectors(M, x, pb)[1])
    dpb = CachedBasis(b, Float64[1, 2, 3, 4, 5, 6], get_vectors(M, x, pb))
    @test sprint(show, "text/plain", dpb) == """
    DiagonalizingOrthonormalBasis(ℝ) with eigenvalue 0 in direction:
     2×3 $(sprint(show, Matrix{Float64})):
       1.0  0.0  0.0
       0.0  0.0  0.0
    and 6 basis vectors.
    Basis vectors:
     E1 =
      2×3 $(sprint(show, Matrix{Float64})):
       1.0  0.0  0.0
       0.0  0.0  0.0
     E2 =
      2×3 $(sprint(show, Matrix{Float64})):
       0.0  0.0  0.0
       1.0  0.0  0.0
     ⋮
     E5 =
      2×3 $(sprint(show, Matrix{Float64})):
       0.0  0.0  1.0
       0.0  0.0  0.0
     E6 =
      2×3 $(sprint(show, Matrix{Float64})):
       0.0  0.0  0.0
       0.0  0.0  1.0
    Eigenvalues:
     6-element $(sprint(show, Vector{Float64})):
      1.0
      2.0
      3.0
      4.0
      5.0
      6.0"""

    M = DefaultManifold(1, 1, 1)
    x = reshape(Float64[1], (1, 1, 1))
    pb = get_basis(M, x, DefaultOrthonormalBasis())
    @test sprint(show, "text/plain", pb) == """
    DefaultOrthonormalBasis(ℝ) with 1 basis vector:
     E1 =
      1×1×1 $(sprint(show, Array{Float64,3})):
      [:, :, 1] =
       1.0"""

    dpb = CachedBasis(
        DiagonalizingOrthonormalBasis(get_vectors(M, x, pb)),
        Float64[1],
        get_vectors(M, x, pb),
    )

    @test sprint(show, "text/plain", dpb) == """
    DiagonalizingOrthonormalBasis(ℝ) with eigenvalue 0 in direction:
     1-element $(sprint(show, Vector{Array{Float64,3}})):
       [1.0]
    and 1 basis vector.
    Basis vectors:
     E1 =
      1×1×1 $(sprint(show, Array{Float64,3})):
      [:, :, 1] =
       1.0
    Eigenvalues:
     1-element $(sprint(show, Vector{Float64})):
      1.0"""
end

ext: jl


text:
module CubedSphere

export sn, cn, conformal_cubed_sphere_mapping, conformal_cubed_sphere_inverse_mapping, cartesian_to_lat_lon

using Printf
using TaylorSeries

include("complex_jacobi_elliptic.jl")
include("rancic_taylor_coefficients.jl")
include("conformal_cubed_sphere.jl")
include("cartesian_to_lat_lon.jl")

end # module


text:
module AtomShell

using Compat; import Compat.String

abstract Shell

_shell = nothing

include("install.jl")
include("process.jl")
include("window.jl")

end


text:
""" Define a schema and ACSet which allows for simulations that schedule events which occur
    after some time delay. When an event occurs, a function called an event listener is called
    on the subset of persons who were scheduled for that event, which may queue state updates or
    schedule further events. This allows simulation of processes with memory (i.e. non-Markovian models).
"""
module SchemaEvents

export TheorySchedulingIBM, AbstractSchedulingIBM, SchedulingIBM,
    add_event, schedule_event, get_scheduled, clear_schedule,
    event_tick, event_process

using Catlab
using Catlab.CategoricalAlgebra
using Catlab.CategoricalAlgebra.FinSets
using Catlab.Present
using Catlab.Theories

using ..SchemaBase

""" A schema for an individual-based model inheriting from `TheoryIBM`
    which allows for events to be scheduled for persons.
"""
@present TheorySchedulingIBM <: TheoryIBM begin
    # set of events which occur after a timed delay
    Event::Ob
    Scheduled::Ob
    scheduled_to_event::Hom(Scheduled, Event)
    scheduled_to_person::Hom(Scheduled, Person)

    Delay::AttrType
    delay::Attr(Scheduled, Delay)

    # labels of events
    EventLabel::AttrType
    eventlabel::Attr(Event, EventLabel)

    # event listeners
    EventListener::AttrType
    eventlistener::Attr(Event, EventListener)
end

""" An abstract ACSet for an individual-based model inheriting from `AbstractIBM`
    which allows for events to be scheduled for persons.
"""
@abstract_acset_type AbstractSchedulingIBM <: AbstractIBM

"""
A concrete ACSet for an individual-based model inheriting from `AbstractSchedulingIBM`
which allows for events to be scheduled for persons.
"""
@acset_type SchedulingIBM(TheorySchedulingIBM, index=[:state, :state_update, :scheduled_to_event, :scheduled_to_person], unique_index = [:statelabel, :eventlabel]) <: AbstractSchedulingIBM

"""
    add_event(model::AbstractSchedulingIBM, label::String, listeners...)

Add an event to the model, with name 'label',
"""
function add_event(model::AbstractSchedulingIBM, label::String, listeners...)
    !(label in subpart(model, :eventlabel)) || throw(ArgumentError("event with name 'label' already assigned in model"))
    listener_vec = Function[]
    for func = listeners
        push!(listener_vec, func)
    end
    add_parts!(model, :Event, 1, eventlabel = label, eventlistener = [listener_vec]);
end

"""
    schedule_event(model::AbstractSchedulingIBM, target, delay, event)

Schedule a set of persons in `target` for the `event` after some `delay`. Note that `event` should correspond to an element
in the set `EventLabel` in your model.
"""
function schedule_event(model::AbstractSchedulingIBM, target, delay, event)
    add_parts!(model, :Scheduled, length(target), scheduled_to_person = target, delay = delay, scheduled_to_event = incident(model, event, :eventlabel))
end

"""
    get_scheduled(model::AbstractSchedulingIBM, event)

Get the set of persons scheduled for `event`. Note that `event` should correspond to an element
in the set `EventLabel` in your model.
"""
function get_scheduled(model::AbstractSchedulingIBM, event)
    model[:scheduled_to_person][incident(model, event, [:scheduled_to_event, :eventlabel])]
end

"""
    clear_schedule(model::AbstractSchedulingIBM, target)

Clear the persons in `target` from any events they are scheduled for.
If 'target' is not specified, clear all scheduled events.
"""
function clear_schedule(model::AbstractSchedulingIBM, target)
    rem_parts!(model, :Scheduled, incident(model, target, [:scheduled_to_person]))
end

function clear_schedule(model::AbstractSchedulingIBM)
    rem_parts!(model, :Scheduled, parts(model, :Scheduled))
end

"""
    event_tick(model::AbstractSchedulingIBM)

Reduce all delays by 1, called at the end of a time step.
"""
function event_tick(model::AbstractSchedulingIBM)
    subpart(model, :delay) .-= 1
end

"""
    event_process(model::AbstractSchedulingIBM, t::Int)

Process events which are ready to fire.
"""
function event_process(model::AbstractSchedulingIBM, t::Int)
    # get every event ready to fire
    ready_to_fire = incident(model, 0, :delay)
    if length(ready_to_fire) > 0
        # apply the event listener for each event type that is ready to fire
        for event = parts(model, :Event)
            fire_this_event = ready_to_fire[model[:scheduled_to_event][ready_to_fire] .== event]
            fire_this_event = model[:scheduled_to_person][fire_this_event]
            # if individuals are scheduled for this event, fire all listeners
            if length(fire_this_event) > 0
                for listener = subpart(model, event, :eventlistener)
                    listener(fire_this_event, t)
                end
            end
        end
        # remove people whose events have fired
        rem_parts!(model, :Scheduled, ready_to_fire)
    end
end

"""
    simulation_loop(model::AbstractIBM, processes::Union{Function, AbstractVector{Function}}, steps::Integer)

A method for simulation models with event scheduling.
"""
function SchemaBase.simulation_loop(model::AbstractSchedulingIBM, processes::Union{Function, AbstractVector{Function}}, steps::Integer)
    if processes isa Function
        processes = [processes]
    end

    apply_attr_updates = create_attr_update(model)
    apply_state_updates = create_state_update(model)

    for t = 1:steps
        for p = processes
            p(t)
        end
        # process events
        event_process(model, t)
        event_tick(model)
        # apply updates
        apply_attr_updates()
        apply_state_updates()
    end
end

end


text:
using Weave
# using Gadfly, Cairo, Fontconfig # Dependencies for default Weave example.
# using IJulia # Dependency for output to notebook.

#=
    Output the default Weave example to different outputs.
=#

# weave(Pkg.dir("Weave","examples","gadfly_md_sample.jmd"), out_path = :pwd, doctype = "md2html") # HTML output.
# weave(Pkg.dir("Weave","examples","gadfly_md_sample.jmd"), out_path = :pwd, doctype = "pandoc") # Markdown output.
# # weave(Pkg.dir("Weave","examples","gadfly_md_sample.jmd"), out_path = :pwd, doctype = "md2pdf") # PDF output. CURRENTLY BROKEN.
# notebook(Pkg.dir("Weave","examples","gadfly_md_sample.jmd"), :pwd) # Notebook output.
# tangle(Pkg.dir("Weave","examples","gadfly_md_sample.jmd"), out_path = :pwd) # Tangled output.

#=
    Test the more specialized features that we need.
=#

# weave("economics_examples.jmd", out_path = :pwd, doctype = "md2html", template = "custom_html.tpl") # HTML output
# # weave("economics_examples.jmd", out_path = :pwd, doctype = "md2pdf") # PDF output. Currently broken.
# convert_doc("economics_examples.jmd", "economics_examples.ipynb") # Notebook output.
# tangle("economics_examples.jmd", out_path = :pwd) # Tangled output.

# Tests for equation numbering, etc.
weave("economics_examples.jmd", out_path = :pwd, doctype = "md2html")
weave("economics_examples.jmd", out_path = :pwd, doctype = "md2tex")


text:
export ZeModule

mutable struct ZeModule
    handle::ze_module_handle_t

    context::ZeContext
    device::ZeDevice

    function ZeModule(ctx::ZeContext, dev::ZeDevice, image; build_flags="", log::Bool=true)
        log_ref = if log
            log_ref = Ref{ze_module_build_log_handle_t}()
        else
            C_NULL
        end

        constants = Ref(ze_module_constants_t(
            0,
            C_NULL,
            C_NULL
        ))

        # compile the module
        GC.@preserve image build_flags constants begin
            desc_ref = Ref(ze_module_desc_t(;
                format=ZE_MODULE_FORMAT_IL_SPIRV,
                inputSize=sizeof(image),
                pInputModule=pointer(image),
                pBuildFlags=pointer(build_flags),
                pConstants=Base.unsafe_convert(Ptr{ze_module_constants_t}, constants)
            ))
            handle_ref = Ref{ze_module_handle_t}()
            res = unsafe_zeModuleCreate(ctx, dev, desc_ref, handle_ref, log_ref)
        end

        # read the log
        if log_ref !== C_NULL
            log_size_ref = Ref{Csize_t}(0)
            zeModuleBuildLogGetString(log_ref[], log_size_ref, C_NULL)
            log_buf = Vector{UInt8}(undef, log_size_ref[])
            zeModuleBuildLogGetString(log_ref[], log_size_ref, pointer(log_buf))
            zeModuleBuildLogDestroy(log_ref[])

            log = String(log_buf)[1:end-1] # strip null terminator
            if !isempty(log)
                if res == ZE_RESULT_ERROR_MODULE_BUILD_FAILURE
                    @error """Module compilation failed:
                              $log"""
                else
                    @debug """Build log:
                              $log"""
                end
            end
        end

        if res != RESULT_SUCCESS
            throw_api_error(res)
        end

        obj = new(handle_ref[], ctx, dev)
        finalizer(obj) do obj
            zeModuleDestroy(obj)
        end
        return obj
    end
end

Base.unsafe_convert(::Type{ze_module_handle_t}, mod::ZeModule) = mod.handle

## kernels

export ZeKernel

mutable struct ZeKernel
    mod::ZeModule
    handle::ze_kernel_handle_t

    function ZeKernel(mod, name)
        GC.@preserve name begin
            desc_ref = Ref(ze_kernel_desc_t(; pKernelName=pointer(name)))
            handle_ref = Ref{ze_kernel_handle_t}()
            zeKernelCreate(mod, desc_ref, handle_ref)
        end
        obj = new(mod, handle_ref[])

        finalizer(obj) do obj
            zeKernelDestroy(obj)
        end
        obj
    end
end

Base.unsafe_convert(::Type{ze_kernel_handle_t}, kernel::ZeKernel) = kernel.handle

## kernel iteration

export kernels

struct ZeModuleKernelDict <: AbstractDict{String,ZeKernel}
    mod::ZeModule
    names::Vector{String}

    function ZeModuleKernelDict(mod)
        count_ref = Ref{UInt32}(0)
        zeModuleGetKernelNames(mod, count_ref, C_NULL)
        names_ref = Vector{Ptr{Cchar}}(undef, count_ref[])
        zeModuleGetKernelNames(mod, count_ref, names_ref)
        new(mod, unsafe_string.(names_ref))
    end
end

kernels(mod::ZeModule) = ZeModuleKernelDict(mod)

function Base.iterate(dict::ZeModuleKernelDict, i=1)
    i > length(dict.names) && return nothing
    name = dict.names[i]
    kernel = ZeKernel(dict.mod, name)
    return (Pair{String,ZeKernel}(name, kernel), i+1)
end

Base.length(dict::ZeModuleKernelDict) = length(dict.names)

function Base.get(dict::ZeModuleKernelDict, name::AbstractString, def)
    in(name, dict.names) || return def
    ZeKernel(dict.mod, name)
end

## group sizes

export ZeDim, suggest_groupsize, groupsize!

"""
    ZeDim3(x)

    ZeDim3((x,))
    ZeDim3((x, y))
    ZeDim3((x, y, x))

A type used to specify dimensions, consisting of 3 integers for respectively the `x`, `y`
and `z` dimension. Unspecified dimensions default to `1`.

Often accepted as argument through the `ZeDim` type alias, allowing to pass dimensions as a
plain integer or a tuple without having to construct an explicit `ZeDim3` object.
"""
struct ZeDim3
    x::Int
    y::Int
    z::Int
end

ZeDim3(dims::Integer)             = ZeDim3(dims,    1,       1)
ZeDim3(dims::NTuple{1,<:Integer}) = ZeDim3(dims[1], 1,       1)
ZeDim3(dims::NTuple{2,<:Integer}) = ZeDim3(dims[1], dims[2], 1)
ZeDim3(dims::NTuple{3,<:Integer}) = ZeDim3(dims[1], dims[2], dims[3])

# Type alias for conveniently specifying the dimensions
# (e.g. `(len, 2)` instead of `ZeDim3((len, 2))`)
const ZeDim = Union{Integer,
                    Tuple{Integer},
                    Tuple{Integer, Integer},
                    Tuple{Integer, Integer, Integer}}

function suggest_groupsize(kernel::ZeKernel, global_sz::ZeDim)
    global_sz = ZeDim3(global_sz)
    group_sz_x = Ref{UInt32}()
    group_sz_y = Ref{UInt32}()
    group_sz_z = Ref{UInt32}()
    zeKernelSuggestGroupSize(kernel, global_sz.x, global_sz.y, global_sz.z,
                             group_sz_x, group_sz_y, group_sz_z)
    return ZeDim3(group_sz_x[], group_sz_y[], group_sz_z[])
end

function groupsize!(kernel::ZeKernel, sz::ZeDim)
    sz = ZeDim3(sz)
    zeKernelSetGroupSize(kernel, sz.x, sz.y, sz.z)
end

## arguments

export arguments

struct ZeKernelArgumentList
    kernel::ZeKernel
end

arguments(kernel::ZeKernel) = ZeKernelArgumentList(kernel)

function Base.setindex!(args::ZeKernelArgumentList, value::Any, index::Integer)
    @assert isbits(value)
    zeKernelSetArgumentValue(args.kernel, index-1, Core.sizeof(value), Base.RefValue(value))
end

## attributes

export indirect_access, indirect_access!, source_attributes

function indirect_access(kernel::ZeKernel)
    flags_ref = Ref{ze_kernel_indirect_access_flags_t}()
    zeKernelGetIndirectAccess(kernel, flags_ref)
    return flags_ref[]
end

indirect_access!(kernel::ZeKernel, flags) = zeKernelSetIndirectAccess(kernel, flags)

function source_attributes(kernel::ZeKernel)
    size_ref = Ref{UInt32}(0)
    zeKernelGetSourceAttributes(kernel, size_ref, C_NULL)

    data = Vector{UInt8}(undef, size_ref[])
    ptr_ref = Ref{Ptr{Cchar}}(pointer(data))
    zeKernelGetSourceAttributes(kernel, size_ref, ptr_ref)
    str = String(data)

    # the attribute string is null-terminated, with attributes separated by space
    return split(str[1:end-1])
end

## properties

export properties

function properties(kernel::ZeKernel)
    props_ref = Ref(ze_kernel_properties_t())
    zeKernelGetProperties(kernel, props_ref)

    props = props_ref[]
    return (
        numKernelArgs=Int(props.numKernelArgs),
        requiredGroupSize=ZeDim3(props.requiredGroupSizeX,
                                 props.requiredGroupSizeY,
                                 props.requiredGroupSizeZ),
        requiredNumSubGroups=Int(props.requiredNumSubGroups),
        requiredSubgroupSize=Int(props.requiredSubgroupSize),
        maxSubgroupSize=Int(props.maxSubgroupSize),
        maxNumSubgroups=Int(props.maxNumSubgroups),
        localMemSize=Int(props.localMemSize),
        privateMemSize=Int(props.privateMemSize),
        spillMemSize=Int(props.spillMemSize),
        kernel_uuid=Base.UUID(reinterpret(UInt128, [props.uuid.kid...])[1]),
        module_uuid=Base.UUID(reinterpret(UInt128, [props.uuid.mid...])[1]),
    )
end

## execution

export append_launch!

ze_group_count_t(dim::ZeDim3) = ze_group_count_t(dim.x, dim.y, dim.z)

function append_launch!(list, kernel, group_count::ZeDim, signal_event=nothing,
                        wait_events::ZeEvent...)
    group_count = ze_group_count_t(ZeDim3(group_count))
    zeCommandListAppendLaunchKernel(list, kernel, Ref(group_count),
                                    something(signal_event, C_NULL),
                                    length(wait_events), [wait_events...])
end


text:
# This file is a part of SimilaritySearch.jl

# Database interface
using Random

export AbstractDatabase, MatrixDatabase, DynamicMatrixDatabase, VectorDatabase, SubDatabase

"""
    abstract type AbstractDatabase end

Base type to represent databases. A database is a collection of objects that can be accessed like a
similar interface to `AbstractVector`. It is separated to allow `SimilaritySearch` methods to know what is a database and what is an object (since most object
representations will look as vectors and matrices).

The basic implementations are:
- `MatrixDatabase`: A wrapper for object-vectors stored in a `Matrix`, columns are the objects. It is static.
- `DynamicMatrixDatabase`: A dynamic representation for vectors that allows adding new vectors.
- `VectorDatabase`: A wrapper for vector-like structures. It can contain any kind of objects.
- `SubDatabase`: A sample of a given database

In particular, the storage details are not used by `VectorDatabase` and `MatrixDatabase`.
For instance, it is possible to use matrices like `Matrix`, `SMatrix` or `StrideArrays`; or even use
generated objects with `VectorDatabase` (supporting a vector-like interface).

If the storage backend support it, it is possible to use vector operations, for example:
- get the `i`-th element `obj = db[i]`, elements in the database are identified by position
- get the elements list in a list of indices `lst` as `db[lst]` (also using `view`)
- set a value at the `i`-th element `db[i] = obj`
- random sampling `rand(db)`, `rand(db, 3)`
- iterate and collect objects in the database
- get the number of elements in the database `length(db)`
- add new objects to the end of the database (not all internal containers will support it)
  - `push!(db, u)` adds a single element `u`
  - `append!(db, lst)` adds a list of objects to the end of the database
"""
abstract type AbstractDatabase end

"""
    view(db::AbstractDatabase, map)

Constructs a `SubDatabase` from `db` using the specified indexes in `map` (e.g., an array or an slice of indices)
"""
@inline Base.view(db::AbstractDatabase, map) = SubDatabase(db, map)
"""
    size(db::AbstractDatabase)

The size of the database `(length(db),)`. There is no concept of explicit dimension for `AbstractDatabase` implementations.
"""
@inline Base.size(db::AbstractDatabase) = (length(db),)

"""
    rand(db::AbstractDatabase)

Retrieves a random element from `db`
"""
@inline Random.rand(db::AbstractDatabase) = @inbounds db[rand(eachindex(db))]

"""
    rand(db::AbstractDatabase, n)

Retrieves `n` random elements from `db`, returnes a `SubDatabase` object
"""
@inline Random.rand(db::AbstractDatabase, n::Integer) = SubDatabase(db, rand(1:length(db), n))

"""
    firstindex(db::AbstractDatabase)

First index in `db` (1)
"""
@inline Base.firstindex(db::AbstractDatabase) = 1

"""
    lastindex(db::AbstractDatabase)

Last index in `db` (1)
"""
@inline Base.lastindex(db::AbstractDatabase) = length(db)

"""
    eachindex(db::AbstractDatabase)

An index iterator of `db`
"""
@inline Base.eachindex(db::AbstractDatabase) = firstindex(db):lastindex(db)

"""
    getindex(db::AbstractDatabase, lst)

A subset of `db` (using a collection of indexes `lst`)
"""
@inline Base.getindex(db::AbstractDatabase, lst::AbstractVector) = SubDatabase(db, lst)
@inline Base.getindex(db::AbstractDatabase, lst::AbstractVector{Bool}) = SubDatabase(db, findall(p -> p > 0, lst))

"""
    struct DynamicMatrixDatabase{DType,Dim} <: AbstractDatabase

A dynamic matrix with elements of type `DType` and dimension `Dim`
"""
struct DynamicMatrixDatabase{DType,Dim} <: AbstractDatabase
    data::Vector{DType}
end

"""
    struct MatrixDatabase{M<:AbstractDatabase} <: AbstractDatabase

    MatrixDatabase(matrix::AbstractMatrix)

Wraps a matrix-like object `matrix` into a `MatrixDatabase`. Please see [`AbstractDatabase`](@ref) for general usage.
"""
struct MatrixDatabase{M<:AbstractMatrix} <: AbstractDatabase
    matrix::M  # abstract matrix
end

"""
    struct VectorDatabase{V} <: AbstractDatabase

A wrapper for vector-like databases
"""
struct VectorDatabase{V} <: AbstractDatabase
    vecs::V  # abstract vector or something that looks like a vector
end

#
# Dynamic matrix-like dataset, i.e., columns are objects. Stored as a large vector to support appending items
#

"""
    eltype(db::DynamicMatrixDatabase{DType,Dim})

The type of stored elements `AbstractVector{DType}`
"""
@inline Base.eltype(db::DynamicMatrixDatabase{DType,Dim}) where {DType,Dim} = AbstractVector{DType}

"""
    DynamicMatrixDatabase(matrix::AbstractMatrix{DType})

Creates a `DynamicMatrixDatabase` from a matrix-like object. Please see [`AbstractDatabase`](@ref) for general usage.
"""
DynamicMatrixDatabase(matrix::AbstractMatrix{DType}) where DType = DynamicMatrixDatabase{DType,size(matrix, 1)}(vec(matrix))

"""
    DynamicMatrixDatabase(::Type{DType}, Dim::Integer)

Creates an empty `DynamicMatrixDatabase` such that `length(db[i]) == Dim` and `eltype(db[i]) == DType`
"""
DynamicMatrixDatabase(::Type{DType}, Dim::Integer) where DType = DynamicMatrixDatabase{DType,Dim}(Vector{DType}(undef, 0))

"""
    DynamicMatrixDatabase(V::DynamicMatrixDatabase{DType,Dim})

Creates a `DynamicMatrixDatabase` from another `DynamicMatrixDatabase`, they will share their internal data. Please see [`AbstractDatabase`](@ref) for general usage.
"""
DynamicMatrixDatabase(V::DynamicMatrixDatabase{DType,Dim}) where {DType,Dim} = DynamicMatrixDatabase{DType,Dim}(V.data)

"""
    DynamicMatrixDatabase(M::MatrixDatabase)

Creates a `DynamicMatrixDatabase` from a `MatrixDatabase`, copies internal data. Please see [`AbstractDatabase`](@ref) for general usage.
"""
DynamicMatrixDatabase(M::MatrixDatabase) = DynamicMatrixDatabase{eltype{M.data},size(M.matrix, 1)}(copy(vec(M.matrix)))

@inline function Base.getindex(db::DynamicMatrixDatabase{DType,Dim}, i::Integer) where {DType,Dim}
    ep = i * Dim
    @inbounds @view db.data[(ep - Dim + 1):ep]
end
@inline Base.setindex!(db::DynamicMatrixDatabase, value, i) = @inbounds (db.data[i] .= value)

@inline Base.length(db::DynamicMatrixDatabase{DType,Dim}) where {DType,Dim} = length(db.data) ÷ Dim
@inline function Base.push!(db::DynamicMatrixDatabase{DType,Dim}, v) where {DType,Dim}
    append!(db.data, v)
    db
end

@inline Base.append!(a::DynamicMatrixDatabase{DType,Dim}, b::DynamicMatrixDatabase{DType2,Dim})  where {DType,Dim,DType2} = append!(a.data, b.data)
@inline function Base.append!(db::DynamicMatrixDatabase{DType,Dim}, B) where {DType,Dim}
    sizehint!(db.data, length(db.data) + Dim * length(B))
    for b in B
        push!(db, b)
    end

    db
end

#
# Wrapper for matrix-like containers
#

@inline Base.eltype(db::MatrixDatabase) = AbstractVector{eltype(db.matrix)}

"""
    MatrixDatabase(V::MatrixDatabase)

Creates another `MatrixDatabase` from another `MatrixDatabase`. They will share their internal data. Please see [`AbstractDatabase`](@ref) for general usage.
"""
MatrixDatabase(V::MatrixDatabase) = MatrixDatabase(V.matrix)

@inline Base.getindex(db::MatrixDatabase, i::Integer) = view(db.matrix, :, i)
@inline Base.setindex!(db::MatrixDatabase, value, i::Integer) = @inbounds (db.matrix[:, i] .= value)
@inline Base.push!(db::MatrixDatabase, v) = error("push! is not supported for MatrixDatabase, please see DynamicMatrixDatabase")
@inline Base.append!(a::MatrixDatabase, b) = error("append! is not supported for MatrixDatabase, please see DynamicMatrixDatabase")
@inline Base.length(db::MatrixDatabase) = size(db.matrix, 2)

#
# Wrapper for array-like containers
#

@inline Base.eltype(db::VectorDatabase) = eltype(db.vecs)

"""
    VectorDatabase(vecs::T)

Creates a `VectorDatabase` from `vecs`. Please see [`AbstractDatabase`](@ref) for general usage.
"""
VectorDatabase(vecs::T) where {T<:AbstractVector} = VectorDatabase{T}(vecs)

"""
    VectorDatabase(M::AbstractMatrix)

Creates a `VectorDatabase` from a matrix-like object. It will copy columns as objects  of the new `VectorDatabase`.
Please see [`AbstractDatabase`](@ref) for general usage.
"""
VectorDatabase(M::T) where {T<:AbstractMatrix} = VectorDatabase([Vector(c) for c in eachcol(M)])

"""
    VectorDatabase(D::AbstractDatabase)

Creates a `VectorDatabase` from an `AbstractDatabase`. It copies internal data.
Please see [`AbstractDatabase`](@ref) for general usage.
"""
VectorDatabase(D::AbstractDatabase) = VectorDatabase([Vector(c) for c in D])

"""
    VectorDatabase(V::VectorDatabase)

Creates a `VectorDatabase` from another an `AbstractDatabase`. They will share their internal data.
Please see [`AbstractDatabase`](@ref) for general usage.
"""
VectorDatabase(V::VectorDatabase) = VectorDatabase(V.vecs)

"""
    VectorDatabase(; type=Vector{Float32})

Creates an empty `VectorDatabase` where each object is of type `type`.
Please see [`AbstractDatabase`](@ref) for general usage.
"""
VectorDatabase(; type=Vector{Float32}) = VectorDatabase(type[])

@inline Base.getindex(db::VectorDatabase, i::Integer) = db.vecs[i]
@inline Base.setindex!(db::VectorDatabase, value, i::Integer) = setindex!(db.vecs, value, i)
@inline Base.length(db::VectorDatabase) = length(db.vecs)
@inline Base.push!(db::VectorDatabase, v) = (push!(db.vecs, v); db)

@inline function Base.append!(db::VectorDatabase, B)
    for b in B
        push!(db, b)
    end

    db
end

# SubDatabase ~ view of the dataset
#
struct SubDatabase{DBType,RType} <: AbstractDatabase
    db::DBType
    map::RType
end

@inline Base.getindex(sdb::SubDatabase, i::Integer) = @inbounds sdb.db[sdb.map[i]]
@inline Base.length(sdb::SubDatabase) = length(sdb.map)
@inline Base.eachindex(sdb::SubDatabase) = eachindex(sdb.map)
@inline function Base.push!(sdb::SubDatabase, v)
    error("push! unsupported operation on SubDatabase")
end
@inline Base.eltype(sdb::SubDatabase) = eltype(sdb.db)
@inline Random.rand(db::SubDatabase, n::Integer) = SubDatabase(db.db, rand(db.map, n))
#
# Generic functions
#
function Base.iterate(db::AbstractDatabase, state::Int=1)
    n = length(db)
    if n == 0 || state > n
        nothing
    else
        @inbounds db[state], state + 1
    end
end

"""
    convert(::Type{AbstractDatabase}, M::AbstractDatabase)
    convert(::Type{AbstractDatabase}, M::AbstractMatrix)
    convert(::Type{AbstractDatabase}, M::Vector)
    convert(::Type{AbstractDatabase}, M::Vector{Any})
    convert(::Type{AbstractDatabase}, M::AbstractVector)
    convert(::Type{<:AbstractVector}, M::VectorDatabase)
    convert(::Type{<:AbstractVector}, M::AbstractDatabase)

Convenience functions to convert different kinds of data into a some kind of database.
"""
Base.convert(::Type{AbstractDatabase}, M::AbstractDatabase) = M
Base.convert(::Type{AbstractDatabase}, M::AbstractMatrix) = MatrixDatabase(M)
Base.convert(::Type{AbstractDatabase}, M::Vector) = VectorDatabase(M)
Base.convert(::Type{AbstractDatabase}, M::Vector{Any}) = VectorDatabase(typeof(first(M)).(M))
Base.convert(::Type{AbstractDatabase}, M::AbstractVector) = VectorDatabase(typeof(first(M)).(M))
Base.convert(::Type{<:AbstractVector}, M::VectorDatabase{T}) where T = M.data
Base.convert(::Type{<:AbstractVector}, M::AbstractDatabase) = collect(M)

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/NNS

max_forks_repo_name: sadit/SimilaritySearch


text:
n = 0
while n < 10
    global n; n += 1
    println(n)
end

for n in 1:10
    println(n)
end

for friend in ["Ted", "Robyn", "Barney","Lily","Marshall"]
    println("Hi $friend, it's great to see you!")
end

#=
let's use for loops to create some addition tables,
where the value of every entry is the sum of its row and column indices.
=#

function create_table1(m,n)
    A = fill(0, (m, n))
    for i in 1:m
        for j in 1:n
            A[i,j] = i+j
        end
    end
    return A
end

#=
Here is some syntactic sugar for the same nested `for` loop
=#

function create_table2(m,n)
    A = fill(0, (m, n))
    for i in 1:m, j in 1:n
        A[i,j] = i+j
    end
    return A
end

#=
The more "Julia" way to create this addition table would have been
with an array comprehension.
=#

function create_table3(m,n)
    C = [i+j for i in 1:m, j in 1:n]
    return C
end

@show create_table1(4,5)
@show create_table2(4,5)
@show create_table3(4,5)

using BenchmarkTools
m,n=1000,1000
@btime create_table1(m,n)
@btime create_table2(m,n)
@btime create_table3(m,n)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
const enum_sym = Symbol("@enum")
const cenum_sym = Symbol("@cenum")
const bitmask_enum_sym = Symbol("@bitmask_flag")
const auto_hash_equals_sym = Symbol("@auto_hash_equals")

isline(x) = false
isline(x::LineNumberNode) = true

function rmlines(ex)
    @match ex begin
        Expr(:macrocall, m, _...) => Expr(:macrocall, m, nothing, filter(!isline, ex.args[3:end])...)
        ::Expr                    => Expr(ex.head, filter(!isline, ex.args)...)
        a                         => a
    end
end

symbols(ex::Expr) = vcat(map(symbols, ex.args)...)
symbols(sym::Symbol) = sym
symbols(x) = []

walk(ex::Expr, inner, outer) = outer(Expr(ex.head, map(inner, ex.args)...))
walk(ex, inner, outer) = outer(ex)

postwalk(f, ex) = walk(ex, x -> postwalk(f, x), f)
prewalk(f, ex) = walk(f(ex), x -> prewalk(f, x), identity)

striplines(ex) = prewalk(rmlines, ex)

function unblock(ex::Expr)
    prewalk(ex) do ex
        ex isa Expr && ex.head == :block && length(ex.args) == 1 ? ex.args[1] : ex
    end
end

function category(ex)
    @match ex begin
        Expr(:struct, _...) || Expr(:macrocall, &auto_hash_equals_sym, _...)                                 => :struct
        Expr(:const, _...)                                                                                => :const
        Expr(:function, _...) || Expr(:(=), Expr(:call, _...) || Expr(:(::), Expr(:call, _...), _...), _...) => :function
        Expr(:macrocall, &enum_sym || &cenum_sym || &bitmask_enum_sym, _...)                                 => :enum
        :(Core.@doc $_ $docstring $ex)                                                                       => :doc
        _                                                                                => nothing
    end
end

prettify(ex) = ex |> striplines |> unblock

isblock(ex) = false
isblock(ex::Expr) = ex.head == :block

broadcast_ex(ex) = Expr(:., ex.args[1], Expr(:tuple, ex.args[2:end]...))
broadcast_ex(ex, cond::Bool) = cond ? broadcast_ex(ex) : ex
broadcast_ex(::Nothing, ::Bool) = nothing

is_broadcast(ex) = @match ex begin
    :($_.($(_...))) => true
    _ => false
end

concat_exs(x) = x
concat_exs(x, y) = Expr(:block, vcat((isblock(x) ? x.args : x), (isblock(y) ? y.args : y))...)
concat_exs(x, y, z...) = foldl(concat_exs, z; init = concat_exs(x, y))

name(sym::Symbol) = sym

function name(ex::Expr)
    @match ex begin
        Expr(:(::), a...)                                                             => name(first(a))
        Expr(:<:, a...)                                                               => name(first(a))
        Expr(:struct, _, _name, _)                                                    => name(_name)
        Expr(:call, f, _...)                                                          => name(f)
        Expr(:., subject, attr, _...)                                                 => name(subject)
        Expr(:function, sig, _...)                                                    => name(sig)
        Expr(:const, assn, _...)                                                      => name(assn)
        Expr(:(=), call, body, _...)                                                  => name(call)
        Expr(:macrocall, &enum_sym || &cenum_sym || &bitmask_enum_sym, _, decl, _...) => name(decl)
        Expr(:kw, _name, _...)                                                        => _name
        :(Core.@doc $_ $docstring $ex)                                                => name(ex)
        Expr(:macrocall, &auto_hash_equals_sym, _, ex)                                => name(ex)
        Expr(expr_type, _...)                                                         => error("Can't extract name from ", expr_type, " expression:\n", "    $ex\n")
    end
end

function name(p::Dict)
    !haskey(p, :category) && return p[:name]
    @match p[:category] begin
        :struct || :enum => name(p[:decl])
        _ => p[:name]
    end
end

type(ex::Expr) = @match ex begin
    :($_ <: $T) || :($_::$T) || :(::$T) => T
    _ => nothing
end

function deconstruct(ex::Expr)
    ex = striplines(ex)
    dict = Dict{Symbol,Any}(:category => category(ex))
    @switch ex begin
        @case Expr(:struct, is_mutable, decl, fields)
        dict[:name] = name(decl)
        dict[:decl] = decl
        dict[:is_mutable] = is_mutable
        dict[:fields] = []
        dict[:constructors] = []
        fields = @match fields begin
            Expr(:block, fields...) => fields
            _ => [fields]
        end
        for field ∈ fields
            field isa LineNumberNode && continue
            @when :function = category(field) begin
                push!(dict[:constructors], deconstruct(field))
                @otherwise
                push!(dict[:fields], field)
            end
        end

        @case Expr(:const, Expr(:(=), sym, val))
        dict[:name] = sym
        dict[:value] = val

        @case Expr(func_sym, :($call::$return_type), body) || Expr(func_sym, :($call::$return_type), body)
        dict = deconstruct(Expr(func_sym, call, body))
        dict[:return_type] = return_type

        @case Expr(:(=), call, body) || Expr(:function, call, body)
        sig_params = deconstruct(call)
        dict[:name] = sig_params[:name]
        dict[:args] = sig_params[:args]
        dict[:kwargs] = sig_params[:kwargs]
        dict[:body] = body
        dict[:short] = ex.head ≠ :function

        @case Expr(:call, f, args...)
        dict[:name] = f
        dict[:args] = []
        dict[:kwargs] = []
        for arg ∈ args
            @switch arg begin
                @case Expr(:parameters, kwargs...)
                append!(dict[:kwargs], kwargs)
                @case _arg
                push!(dict[:args], arg)
            end
        end

        @case :(Core.@doc $_ $docstring $ex)
        dict[:docstring] = docstring
        dict[:ex] = ex

        @case Expr(:macrocall, m, _, decl, args...)
        dict[:macro] = m
        dict[:name] = name(decl)
        dict[:type] = decl isa Symbol ? Int : type(decl)
        dict[:decl] = decl
        dict[:values] = @match args begin
            [::Expr] => rmlines(args[1].args)
            x        => x
        end

        @case _
        error("Matching non-exhaustive for expr $ex\n$(dump(ex))")
    end
    dict
end

function reconstruct_call(d::Dict; is_decl = true, with_typeassert=true)
    args = get(d, :args, [])
    kwargs = get(d, :kwargs, [])

    if !is_decl
        args = name.(args)
        kwargs = name.(kwargs)
    end

    call = @match (args, kwargs) begin
        (args, []) => Expr(:call, name(d), args...)
        (args, kwargs) => Expr(:call, name(d), Expr(:parameters, kwargs...), args...)
    end
    rt = get(d, :return_type, nothing)
    if isnothing(rt) || !is_decl || !with_typeassert
        call
    else
        :($call::$rt)
    end
end

function reconstruct_documented(d::Dict, ex)
    docstring = get(d, :docstring, "")
    isempty(docstring) ? ex : :(Core.@doc $docstring $ex)
end

function reconstruct(d::Dict)
    category = d[:category]
    ex = @match category begin
        :struct   => begin
            props = [get(d, :fields, []); get(d, :constructors, [])]
            Expr(:struct, get(d, :is_mutable, false), d[:decl], Expr(:block, props...))
        end
        :const    => :(const $(d[:name]) = $(d[:value]))
        :enum     => Expr(:macrocall, d[:macro], nothing, d[:decl], Expr(:block, d[:values]...))
        :function => begin
            call = reconstruct_call(d)
            get(d, :short, false) ? :($call = $(d[:body])) : Expr(:function, call, d[:body])
        end
        :doc      => :(Core.@doc $(d[:docstring]) $(d[:ex]))
        _         => error("Category $category cannot be constructed")
    end
    if category == :enum
        ex
    else
        unblock(ex)
    end
end

function to_expr(p::Dict)
    if p[:category] == :function
        if haskey(p, :relax_signature)
            val = p[:relax_signature]
            if val isa Bool
                val && relax_function_signature!(_ -> true, p)
            else
                relax_function_signature!(p[:relax_signature], p)
            end
        end
    end
    reconstruct(p)
end


text:
"""
```
Cloud
```

The `Cloud` type contains all of the relevant information for a given cloud of
particles in the SMC algorithm. Information for a single iteration is stored at
any given time (and thus the final output will be the final cloud of particles,
 of which only the particle values will be saved).

### Fields
- `particles::Matrix{Float64}`: The vector of particles. The first `1:(size(particles, 2) - 5)`
    columns hold the parameter values. The last five columns
    contain "metadata" on the particles, namely (in order from left to right)
    (1) the log-likelihood, (2) log-prior, (3) the old log-likelihood from a previous estimation
    (used for tempered updates, set to zero if not tempering), (4) the accept rate
    for the particle's mutation steps, and (5) the particle's normalized weight.
- `tempering_schedule::Vector{Float64}`: The vector of ϕ_ns (tempering factors)
- `ESS::Vector{Float64}`: The vector of effective sample sizes (resample if ESS
    falls under the threshold)
- `stage_index::Int`: The current iteration index of the algorithm
- `n_Φ::Int`: The total number of stages of in the fixed tempering schedule
    (if the algorithm is run with an adaptive ϕ schedule then this is used to
    calibrate the ϕ_prop)
- `resamples::Int`: The number of times the particle population was resampled
- `c::Float64`: The mutation step size
- `accept::Float64`: The average acceptance rate of mutation steps
- `total_sampling_time::Float64`: Total amount of time that SMC algorithm took
     to execute
"""
mutable struct Cloud
    particles::Matrix{Float64}
    tempering_schedule::Vector{Float64}
    ESS::Vector{Float64}
    stage_index::Int
    n_Φ::Int
    resamples::Int
    c::Float64
    accept::Float64
    total_sampling_time::Float64
end

"""
```
Cloud(n_params::Int, n_parts::Int)
```
Easier constructor for Cloud, which initializes the weights to be
equal, and everything else in the particle object to be empty.
"""
function Cloud(n_params::Int, n_parts::Int)
    return Cloud(Matrix{Float64}(undef, n_parts, n_params + 5),
                 zeros(1), zeros(1), 1, 0, 0, 0., 0.25, 0.)
end

"""
Find correct indices for accessing columns of cloud array.
"""
ind_para_end(N::Int)  = N-5
ind_loglh(N::Int)     = N-4
ind_logprior(N::Int)  = N-3
ind_old_loglh(N::Int) = N-2
ind_accept(N::Int)    = N-1
ind_weight(N::Int)    = N

"""
```
get_weights(c::Matrix{Float64})
```
Returns Vector{Float64}(n_parts) of weights of particles in cloud.
"""
function get_weights(c::Matrix{Float64})
    return c[:, ind_weight(size(c,2))]
end
"""
```
get_weights(c::Cloud)
```
Returns Vector{Float64}(n_parts) of weights of particles in cloud.
"""
function get_weights(c::Cloud)
    return c.particles[:, ind_weight(size(c.particles,2))]
end

"""
```
get_vals(c::Cloud; transposed::Bool = true)
```
Returns Matrix{Float64}(n_params, n_parts) of parameter values in particle cloud.
"""
function get_vals(c::Matrix{Float64}; transpose::Bool = true)
    return transpose ? Matrix{Float64}(c[:, 1:ind_para_end(size(c, 2))]') :
                                       c[:, 1:ind_para_end(size(c, 2))]
end
"""
```
get_vals(c::Matrix{Float64})
```
Returns Matrix{Float64}(n_params, n_parts) of parameter values in particle cloud.
"""
function get_vals(c::Cloud; transpose::Bool = true)
    return transpose ? Matrix{Float64}(c.particles[:, 1:ind_para_end(size(c.particles,2))]') :
                                       c.particles[:, 1:ind_para_end(size(c.particles,2))]
end

"""
```
get_loglh(c::Matrix{Float64})
get_loglh(c::Cloud)
```
Returns Vector{Float64}(n_parts) of log-likelihood of particles in cloud.
"""
function get_loglh(c::Matrix{Float64})
    return c[:, ind_loglh(size(c,2))]
end
function get_loglh(c::Cloud)
    return c.particles[:, ind_loglh(size(c.particles,2))]
end

"""
```
cloud_isempty(c::Matrix{Float64})
cloud_isempty(c::Cloud)
```
Check if cloud has no particles.
"""
function cloud_isempty(c::Matrix{Float64})
    return isempty(c)
end
function cloud_isempty(c::Cloud)
    # TODO: 'ISEMPTY' IS INCOMPATIBLE WITH INITIALIZATION
    return isempty(c.particles)
end

"""
```
get_old_loglh(c::Matrix{Float64})
get_old_loglh(c::Cloud)
```
Returns Vector{Float64}(n_parts) of old log-likelihood of particles in cloud.
"""
function get_old_loglh(c::Matrix{Float64})
    return c[:, ind_old_loglh(size(c,2))]
end
function get_old_loglh(c::Cloud)
    return c.particles[:, ind_old_loglh(size(c.particles,2))]
end

"""
```
get_logpost(c::Matrix{Float64})
get_logpost(c::Cloud)
```
Returns Vector{Float64}(n_parts) of log-posterior of particles in cloud.
"""
function get_logpost(c::Matrix{Float64})
    return c[:, ind_loglh(size(c,2))] .+ c[:, ind_logprior(size(c,2))]
end
function get_logpost(c::Cloud)
    return c.particles[:, ind_loglh(size(c.particles,2))] .+
        c.particles[:, ind_logprior(size(c.particles,2))]
end

"""
```
get_logprior(c::Matrix{Float64})
get_logprior(c::Cloud)
```
Returns Vector{Float64}(n_parts) of log-prior of particles in cloud.
"""
function get_logprior(c::Matrix{Float64})
    return c[:, ind_logprior(size(c,2))]
end
function get_logprior(c::Cloud)
    return c.particles[:, ind_logprior(size(c.particles,2))]
end

"""
```
get_accept(c::Matrix{Float64})
get_accept(c::Cloud)
```
Returns Vector{Float64}(n_parts) of old log-likelihood of particles in cloud.
"""
function get_accept(c::Matrix{Float64})
    return c[:, ind_accept(size(c,2))]
end
function get_accept(c::Cloud)
    return c.particles[:, ind_accept(size(c.particles,2))]
end

"""
```
get_likeliest_particle_value(c::Matrix{Float64})
get_likeliest_particle_value(c::Cloud)
```
Return parameter vector of particle with highest log-likelihood.
"""
function get_likeliest_particle_value(c::Matrix{Float64})
    return c[argmax(get_loglh(c)), 1:ind_para_end(size(c,2))]
end
function get_likeliest_particle_value(c::Cloud)
    return c.particles[argmax(get_loglh(c)), 1:ind_para_end(size(c.particles,2))]
end

"""
```
get_highest_posterior_particle_value(c::Matrix{Float64})
get_highest_posterior_particle_value(c::Cloud)
```
Return parameter vector of particle with highest log-posterior.
"""
function get_highest_posterior_particle_value(c::Matrix{Float64})
    return c[argmax(get_logpost(c)), 1:ind_para_end(size(c,2))]
end
function get_highest_posterior_particle_value(c::Cloud)
    return c.particles[argmax(get_logpost(c)), 1:ind_para_end(size(c.particles,2))]
end

"""
```
update_draws!(c::Cloud, draws::Matrix{Float64})
```
Update parameter draws in cloud.
"""
function update_draws!(c::Cloud, draws::Matrix{Float64})
    I, J     = size(draws)
    n_parts  = length(c)
    n_params = ind_para_end(size(c.particles, 2))
    if (I, J) == (n_parts, n_params)
        for i = 1:I, j=1:J
            c.particles[i, j] = draws[i, j]
        end
    elseif (I, J) == (n_params, n_parts)
        for i = 1:I, j=1:J
            c.particles[j, i] = draws[i, j]
        end
    else
        throw(error("update_draws!(c::Cloud, draws::Matrix): Draws are incorrectly sized!"))
    end
end

"""
```
update_weights!(c::Matrix{Float64}, incweight::Vector{Float64})
update_weights!(c::Cloud, weights::Vector{Float64})
```
Update weights in cloud.
"""
function update_weights!(c::Matrix{Float64}, incweight::Vector{Float64})
    @assert size(c, 1) == length(incweight) "Dimensional mismatch in inc. weights"
    N = ind_weight(size(c,2))
    for i=1:length(incweight)
        c[i, N] *= incweight[i]
    end
end
function update_weights!(c::Cloud, weights::Vector{Float64})
    update_weights!(c.particles, weights)
end

"""
```
set_weights!(c::Cloud, weights::Vector{Float64})
```
Set weights to specific values. Contrast to update_weights, which multiplies
existing weights by provided incremental weights.
"""
function set_weights!(c::Cloud, weights::Vector{Float64})
    @assert length(c) == length(weights) "Dimensional mismatch in set_weights"
    N = ind_weight(size(c.particles,2))
    for i=1:length(c)
        c.particles[i, N] = weights[i]
    end
end

"""
```
update_loglh!(c::Matrix{Float64}, incweight::Vector{Float64})
update_loglh!(c::Cloud, loglh::Vector{Float64})
```
Update log-likelihood in cloud.
"""
function update_loglh!(c::Matrix{Float64}, loglh::Vector{Float64})
    @assert size(c,1) == length(loglh) "Dimensional mismatch"
    N = ind_loglh(size(c,2))
    for i=1:length(loglh)
        c[i, N] = loglh[i]
    end
end
function update_loglh!(c::Cloud, loglh::Vector{Float64})
    update_loglh!(c.particles, loglh)
end

"""
```
update_logprior!(c::Matrix{Float64}, incweight::Vector{Float64})
update_logprior!(c::Cloud, logprior::Vector{Float64})
```
Update log-prior in cloud.
"""
function update_logprior!(c::Matrix{Float64}, logprior::Vector{Float64})
    @assert size(c, 1) == length(logprior) "Dimensional mismatch"
    N = ind_logprior(size(c,2))
    for i=1:length(logprior)
        c[i, N] = logprior[i]
    end
end
function update_logprior!(c::Cloud, logprior::Vector{Float64})
    update_logprior!(c.particles, logprior)
end

"""
```
update_old_loglh!(c::Matrix{Float64}, incweight::Vector{Float64})
update_old_loglh!(c::Cloud, old_loglh::Vector{Float64})
```
Update log-likelihood in cloud.
"""
function update_old_loglh!(c::Matrix{Float64}, old_loglh::Vector{Float64})
    @assert size(c, 1) == length(old_loglh) "Dimensional mismatch"
    N = ind_old_loglh(size(c,2))
    for i=1:length(old_loglh)
        c[i, N] = old_loglh[i]
    end
end
function update_old_loglh!(c::Cloud, old_loglh::Vector{Float64})
    update_old_loglh!(c.particles, old_loglh)
end

"""
```
adjust_weights!(c::Matrix{Float64}, loglh::Vector{Float64}, logprior::Vector{Float64})
adjust_weights!(c::Cloud, loglh::Vector{Float64}, logprior::Vector{Float64})
```
Update the weights based on the loglh and logprior.
Not actually used anywhere in the code.
"""
function adjust_weights!(c::Matrix{Float64}, loglh::Vector{Float64}, logprior::Vector{Float64})
    @assert size(c, 1) == length(loglh) && size(c, 1) == length(logprior) "Dimensional mismatch"
    N = ind_weight(size(c,2))
    for i=1:length(loglh)
        c[i, N] = loglh[i] + logprior[i]
    end
    c[:,N] = exp.(c[:,N])
    c[:,N] ./= sum(c[:,N])
    c[:,N] .*= length(loglh)
end

function adjust_weights!(c::Cloud, loglh::Vector{Float64}, logprior::Vector{Float64})
    adjust_weights!(c.particles, loglh, logprior)
end

"""
```
normalize_weights!(c::Matrix{Float64})
normalize_weights!(c::Cloud)
```
Normalize weights in cloud to N, the number of particles.
"""
function normalize_weights!(c::Matrix{Float64})
    sum_weights = sum(get_weights(c))
    c[:, ind_weight(size(c,2))] *= size(c,1)
    c[:, ind_weight(size(c,2))] /= sum_weights
end
function normalize_weights!(c::Cloud)
    normalize_weights!(c.particles)
end

"""
```
reset_weights!(c::Matrix{Float64})
reset_weights!(c::Cloud)
```
Uniformly reset weights of all particles to 1/n_parts.
"""
function reset_weights!(c::Matrix{Float64})
    c[:, ind_weight(size(c,2))] .= 1.0
end
function reset_weights!(c::Cloud)
    reset_weights!(c.particles)
end

"""
```
zero_bad_loglh_weights!(c::Matrix{Float64})
zero_bad_loglh_weights!(c::Cloud)
```
puts zero weight on particles with -Inf log-likelihoods
"""
function zero_bad_loglh_weights!(c::Matrix{Float64})
    n_cloud_cols = size(c, 2)
    badloglh = findall(c[:, ind_loglh(n_cloud_cols)] .== -Inf)
    c[badloglh, ind_weight(n_cloud_cols)] .= 0.
end
function zero_bad_loglh_weights!(c::Cloud)
    zero_bad_loglh_weights!(c.particles)
end

"""
```
update_mutation!(p::Vector{Float64}, para::Vector{Float64},
                          like::Float64, prior::Float64, old_like::Float64,
                          accept::Float64)
```
Update a particle's parameter vector, log-likelihood, log-prior,
old log-likelihood, and acceptance rate at the end of mutation.
"""
function update_mutation!(p::Vector{Float64}, para::Vector{Float64}, like::Float64,
                          prior::Float64, old_like::Float64, accept::Float64)
    N = length(p)
    p[1:ind_para_end(N)] = para
    p[ind_loglh(N)]      = like
    p[ind_logprior(N)]   = prior
    p[ind_old_loglh(N)]  = old_like
    p[ind_accept(N)]     = accept
end

"""
```
update_cloud!(cloud::Cloud, new_particles::Matrix{Float64})
```
Updates cloud values with those of new particles in particle array.
"""
function update_cloud!(cloud::Cloud, new_particles::Matrix{Float64})
    I, J = size(new_particles)
    if I == length(cloud)
        cloud.particles = new_particles
    elseif J == length(cloud)
        for k = 1:length(cloud)
            cloud.particles[k,:] = new_particles[:, k]
        end
    else
        throw(error("update_cloud!(c::Cloud, draws): draws are incorrect size!"))
    end
end

"""
```
update_val!(p::Vector{Float64}, val::Vector{Float64})
```
Update parameter vector of particle.
"""
function update_val!(p::Vector{Float64}, val::Vector{Float64})
    @assert ind_para_end(length(p)) == length(val) "Parameter vector length is wrong!"
    p[1:length(val)] = val
end

"""
```
update_weight!(p::Vector{Float64}, weight::Vector{Float64})
```
Update weight of particle.
"""
function update_weight!(p::Vector{Float64}, weight::Float64)
    p[ind_weight(length(p))] = weight
end

"""
```
update_acceptance_rate!(c::Cloud)
```
Update cloud's acceptance rate with the mean of its particle acceptance rates.
"""
function update_acceptance_rate!(c::Cloud)
    c.accept = mean(get_accept(c))
end

function Base.length(c::Cloud)
    return size(c.particles, 1)
end

"""
```
weighted_mean(c::Cloud)
weighted_mean(c::Matrix{Float64})
```
Compute weighted mean of particle cloud.
"""
function weighted_mean(c::Matrix{Float64})
    return get_vals(c) * get_weights(c) / sum(get_weights(c))
end
function weighted_mean(c::Cloud)
    return weighted_mean(c.particles)
end

"""
```
weighted_quantile(c::Matrix{Float64}, i::Int64)
weighted_quantile(c::Cloud, i::Int64)
```
Compute weighted quantiles of particle cloud for input parameter, indexed by i.
"""
function weighted_quantile(c::Matrix{Float64}, i::Int64)
    @assert i <= ind_end_para(size(c,2)) "Parameter index invalid."
    lb = quantile(c[:, i], Weights(get_weights(c)), .05)
    ub = quantile(c[:, i], Weights(get_weights(c)), .95)
    return lb, ub
end
function weighted_quantile(c::Cloud, i::Int64)
    return weighted_quantile(c.particles, i)
end

"""
```
weighted_std(c::Cloud)
weighted_std(c::Matrix{Float64})
```
Compute weighted standard deviation of particle cloud.
"""
function weighted_std(c::Matrix{Float64})
    return sqrt.(diag(weighted_cov(c)))
end
function weighted_std(c::Cloud)
    return weighted_std(c.particles)
end

"""
```
weighted_cov(c::Cloud)
weighted_cov(c::Matrix{Float64})
```
Compute weighted covariance of particle cloud.
"""
function weighted_cov(c::Matrix{Float64})
    return cov(get_vals(c; transpose = false),
               Weights(get_weights(c) / sum(get_weights(c))), corrected = false)
end
function weighted_cov(c::Cloud)
    return weighted_cov(c.particles)
end

"""
```
split_cloud(filename::String, n_pieces::Int)
```
splits the cloud saved in `filename` to multiple files (`n_pieces` of them).
For large clouds, the memory usage by one file may be too large. For example,
the file may exceed GitHub's 100MB memory limit.
"""
function split_cloud(filename::String, n_pieces::Int)
    cloud = load(filename, "cloud")
    w    = load(filename, "w")
    W    = load(filename, "W")
    n_part = size(cloud.particles, 1)
    n_para = size(cloud.particles, 2)

    @assert mod(n_part, n_pieces) == 0
    npart_small = Int(n_part / n_pieces)
    clouds = Vector{Cloud}(undef, n_pieces)
    ws     = Vector{Matrix{Float64}}(undef, n_pieces)
    Ws     = Vector{Matrix{Float64}}(undef, n_pieces)

    for i = 1:n_pieces
        inds = ((i-1)*npart_small+1):(i*npart_small)
        ws[i] = w[inds, :]
        Ws[i] = W[inds, :]
        clouds[i] = Cloud(n_para, npart_small)

        # Since loglh, logprior, old_loglh are all stored in cloud.particles, this updates them all too!
        SMC.update_cloud!(clouds[i], cloud.particles[inds, :])

        clouds[i].ESS = cloud.ESS
        clouds[i].c = cloud.c
        clouds[i].stage_index = cloud.stage_index
        clouds[i].total_sampling_time = cloud.total_sampling_time
        clouds[i].accept = cloud.accept
        clouds[i].n_Φ = cloud.n_Φ
        clouds[i].resamples = cloud.resamples
        clouds[i].tempering_schedule = cloud.tempering_schedule

        new_filename = replace(filename, ".jld2" => "_part$(i).jld2")
        jldopen(new_filename, true, true, true, IOStream) do file
            write(file, "cloud", clouds[i])
            write(file, "w", ws[i])
            write(file, "W", Ws[i])
        end
    end
end

"""
```
join_cloud(filename::String, n_pieces::Int; save_cloud::Bool = true)
```
joins a cloud saved in `filename` that had previously been
split into multiple files (`n_pieces` of them) and saves the newly
joined cloud when the kwarg `save_cloud` is `true`.
For large clouds, the memory usage by one file may be too large. For example,
the file may exceed GitHub's 100MB memory limit. For this reason,
it is useful to split the cloud file into multiple files and then
rejoin later.
"""
function join_cloud(filename::String, n_pieces::Int; save_cloud::Bool = true)
    clouds     = Vector{Cloud}(undef, n_pieces)
    ws         = Vector{Matrix{Float64}}(undef, n_pieces)
    Ws         = Vector{Matrix{Float64}}(undef, n_pieces)
    for i = 1:n_pieces
        small_filename = replace(filename, ".jld2" => "_part$(i).jld2")

        clouds[i]     = load(small_filename, "cloud")
        ws[i]         = load(small_filename, "w")
        Ws[i]         = load(small_filename, "W")
    end
    n_part  = sum([size(y.particles, 1) for y in clouds])
    n_para  = size(clouds[1].particles, 2)
    n_stage = size(ws[1], 2)

    cloud = Cloud(n_para, n_part)

    particles  = Matrix{Float64}(undef, 0, n_para)
    loglhs     = Vector{Float64}(undef, 0)
    old_loglhs = Vector{Float64}(undef, 0)
    logpriors  = Vector{Float64}(undef, 0)
    w          = Matrix{Float64}(undef, 0, n_stage)
    W          = Matrix{Float64}(undef, 0, n_stage)

    for i = 1:n_pieces
        particles     = vcat(particles, clouds[i].particles)
        loglhs        = vcat(loglhs, SMC.get_loglh(clouds[i]))
        old_loglhs    = vcat(old_loglhs, SMC.get_old_loglh(clouds[i]))
        logpriors     = vcat(logpriors, SMC.get_logprior(clouds[i]))
        w             = vcat(w, ws[i])
        W             = vcat(W, Ws[i])
    end

    # Since loglh, logprior, old_loglh are all stored in cloud.particles, this updates them all too!
    SMC.update_cloud!(cloud, particles)

    cloud.ESS                 = clouds[1].ESS
    cloud.c                   = clouds[1].c
    cloud.stage_index         = clouds[1].stage_index
    cloud.total_sampling_time = clouds[1].total_sampling_time
    cloud.accept              = clouds[1].accept
    cloud.n_Φ                 = clouds[1].n_Φ
    cloud.resamples           = clouds[1].resamples
    cloud.tempering_schedule  = clouds[1].tempering_schedule

    if save_cloud
        jldopen(filename, true, true, true, IOStream) do file
            write(file, "cloud", cloud)
            write(file, "w", w)
            write(file, "W", W)
        end
    end

    return cloud
end

"""
```
add_parameters_to_cloud(old_cloud_file::String, para::ParameterVector,
                        old_para_inds::BitVector; regime_switching::Bool = false)

add_parameters_to_cloud(old_cloud::Cloud, para::ParameterVector, old_para_inds::BitVector;
                        regime_switching::Bool = false) where {T <: Real}
```
extends a `Cloud` from a previous estimation to include new parameters.
This function helps construct a bridge distribution when
you want to estimate a model that extends a previous model by
adding additional parameters.

To be concrete, suppose we have two models ``\\mathcal{M}_1``
and ``\\mathcal{M}_2`` such that the parameters of the first model
are a subset of the parameters of the second model. For example,
suppose ``\\theta_1`` are the parameters for the first model,
and ``\\theta_2 = [\\theta_1, \\tilde{\\theta}]^T``,
where ``\\tilde{\\theta}`` are the new parameters for ``\\mathcal{M}_2``.
Assume that

(1) the likelihood function for ``\\mathcal{M}_1` does not depend on ``\tilde{\\theta}``,
(2) the priors for ``\theta_1`` and ``\tilde{\theta}`` are independent.

Then the posterior for ``\\theta_2`` given ``\\mathcal{M}_1`` is just
``math
\\begin{aligned}
  \\pi(\\theta_2 \\vert Y, \\mathcal{M}_1) = \\pi(\\theta_1 \\vert Y, \\mathcal{M}_1) p(\\tilde{\\theta}).
\\end{aligned}
``
Therefore, we can construct the posterior ``\\pi(\\theta_2 \\vert Y, \\mathcal{M}_1)``
by concatenating draws from the prior for ``\\tilde{\\theta}`` to the previous estimation of ``\\mathcal{M}_1``.
This function performs this concatenation and allows for regime-switching.

### Inputs
- `old_cloud` or `old_cloud_file`: this input should specify the `Cloud` from a previous estimation
    of the old model. If a different package was used to estimate the model, then the user
    can construct a new `Cloud` object. The only fields which must be set correctly
    are the `particles` and `ESS` field. The others can be set to default values (see `?Cloud` and
    the source code for the construction of the `Cloud` object).

- `para`: the parameter vector of the new model (i.e. ``\\theta_2``).
    The parameters in `para` that belonged to the old model should have the same
    settings as the ones used in previous estimation, e.g. the same prior.

- `old_para_inds`: indicates which parameters were used in the old estimation.
    If `regime_switching = true`, then this vector should specify which values
    correspond to old parameter values based on the matrix returned by
    `SMC.get_values(para)`. For example, if `old_para` is the `ParameterVector`
    used by the old estimation, then `SMC.get_values(old_para) == SMC.get_values(para)[old_para_inds]`

### Keyword Arguments
- `regime_switching`: this kwarg is needed to be able to draw from the prior correctly
    for the new parameters.
"""
function add_parameters_to_cloud(old_cloud_file::String, para::ParameterVector,
                                 old_para_inds::BitVector; regime_switching::Bool = false)
    old_cloud = load(old_cloud_file, "cloud")
    return add_parameters_to_cloud(old_cloud, para, old_para_inds;
                                   regime_switching = regime_switching)
end

function add_parameters_to_cloud(old_cloud::Cloud, para::ParameterVector{T}, old_para_inds::BitVector;
                                 regime_switching::Bool = false) where {T <: Real}

    # Sample from prior
    n_parts   = length(old_cloud)
    para_vals = regime_switching ? Matrix{T}(undef, n_parts, n_parameters_regime_switching(para)) :
        Matrix{T}(undef, n_parts, length(para))
    for i in 1:n_parts
        para_vals[i, :] = rand(para; regime_switching = regime_switching)
    end

    # Use same particles as the old cloud. Re-sampling is unnecessary
    # because old_cloud.particles has the needed weights information
    part_dim2 = size(old_cloud.particles, 2)
    old_para  = old_cloud.particles[:, 1:ind_para_end(part_dim2)]

    # Combine old parameters (drawn from their posterior from the old estimation)
    # and the new parameters (drawn from a prior)
    # It is assumed that the order of old_para matches the order of old_para_inds
    # and that the order of parameters haven't been switched around.
    # If that is the case, then the user need to write a function
    # that maps the old parameters to their correct indices.
    para_vals[:, old_para_inds] = old_para

    # Create logprior columns
    meta_info = Matrix{T}(undef, n_parts, 5) # additional 5 columns of "meta" information about particles
    for i in 1:n_parts
        # Update ParameterVector para
        update!(para, view(para_vals, i, :))

        # Compute logprior
        meta_info[i, 2] = prior(para)
    end

    # Copy loglh from old model
    meta_info[:, 1] = view(old_cloud.particles, :, ind_loglh(part_dim2))

    # Compute old loglh (just zeros since we will be creating a "new" Cloud)
    meta_info[:, 3] .= 0.

    # Add acceptance rate based on the bold cloud
    meta_info[:, 4] = view(old_cloud.particles, :, ind_accept(part_dim2))

    # Add weights for each particle based on the old cloud
    meta_info[:, 5] = view(old_cloud.particles, :, ind_weight(part_dim2))

    # Form a new Cloud
    return Cloud(hcat(para_vals, meta_info), zeros(1), old_cloud.ESS, 1, 0, 0, 0., .25, 0.)
end


text:
@testset "IO" begin
    @testset "FileExist" begin
        @test isfile("data/in.cfl")
        @test isfile("data/in.hdr")
    end

    @testset "ReadWrite" begin
        ## Read bart files
        println(pwd())
        filenameIn = "data/in"
        im = BartIO.readcfl(filenameIn)

        ## test writing to file
        filenameOut="data/out"
        BartIO.writecfl(filenameOut,im)
        ## Read back that files
        im2 = BartIO.readcfl(filenameOut)
        rm(filenameOut*".hdr")
        rm(filenameOut*".cfl")

        ## test
        @test im == im2

    end
end

@testset "BART" begin
    pathtobart = "/home/runner/work/BartIO.jl/BartIO.jl/bart"
    pathtobartpy = "/home/runner/work/BartIO.jl/BartIO.jl/bartpy"

    @testset "BART_files" begin
        @test isdir(pathtobart)
        @test isdir(pathtobartpy)

        strpy = readchomp(`python3 script.py`)
        @test strpy == "hello world"
    end

    @testset "BART_exec" begin
        #=
        python_pycall = PyCall.python

        PyCall.py"""
        import os
        os.environ['TOOLBOX_PATH'] = $pathtobart
        print(os.environ['TOOLBOX_PATH'])
        os.chdir($pathtobartpy)
        os.system($python_pycall + " setup.py install --user")
        """

        #@PyCall.pyimport bartpy.tools as bartpy #Equivalent to -> bartpy = pyimport("bartpy.tools") but does not work in module...
        bartpy = pyimport("bartpy.tools")
        bartpy.version()
        =#

        bartpy = BartIO.initBart(pathtobart,pathtobartpy)
        @test typeof(bartpy) == PyObject
        @test size(bartpy.phantom()) == (128, 128)

        bartpy2 = BartIO.initBart()
        @test size(bartpy2.phantom()) == (128, 128)
    end
end


text:
module HandleInvoke

using Mixtape
using MacroTools
using BenchmarkTools

@noinline foo(x) = x^5
@noinline bar(x) = x^10
apply(f, x1, x2::Val{T}) where {T} = f(x1, T)

@noinline function f(x)
    g = x < 5 ? foo : bar
    return invoke(g, Tuple{Int}, 2)
end

@ctx (false, false, false) struct MyMix end

allow(ctx::MyMix, m::Module) = m == HandleInvoke
show_after_inference(ctx::MyMix) = false
show_after_optimization(ctx::MyMix) = false
debug(ctx::MyMix) = false

swap(e) = e
function swap(e::Expr)
    new = MacroTools.postwalk(e) do s
        isexpr(s, :call) || return s
        s.args[1] == Base.literal_pow || return s
        return Expr(:call, apply, Base.:(*), s.args[3:end]...)
    end
    return new
end

function transform(::MyMix, b)
    display(b)
    b = Mixtape.widen_invokes!(b)
    for (v, st) in b
        replace!(b, v, swap(st))
    end
    return b
end

Mixtape.@load_call_interface()
display(@time call(MyMix(), f, 3))

end # module


text:
function read_VolumeImage(fname::AbstractString)

    Logging.info("Creating volume image from file $fname")

    if contains(fname, ".dat")
        x, y, z, s, t = read_dat(fname)
        method = "CLARA"
        header = Dict()
        units = "nAm/cm^3"
        x = x/1000
        y = y/1000
        z = z/1000
        t = t/1000
    else
        Logging.warn("Unknown file type")
    end

    header["FileName"] = fname

    coord_system = "?"

    VolumeImage(s, units, collect(x), collect(y), collect(z), collect(t), method, header, coord_system)
end


text:
# Demographics

function init_demographics()
    d = Demographics(60, 45, 1.0);
    @assert validate_demographics(d);
    return d
end

function validate_demographics(d :: Demographics)
    isValid = (d.lifeSpan > d.lastWorkAge);
    return isValid
end

## ---------  Helpers

lifespan(d :: Demographics) = d.lifeSpan;
workspan(d :: Demographics) = d.lastWorkAge;
isretired(d :: Demographics, t :: Integer) = (t > workspan(d));

# ------------


text:
@with_kw mutable struct Labels
  RESPONDER_PACKAGE_NAME::String = ""
  RESPONDER_FUNCTION_NAME::String = ""
  RESPONDER_COMMIT::Union{Nothing, String} = nothing
  RESPONDER_TREE_HASH::String = ""
  RESPONDER_PKG_SOURCE::Union{Nothing, String} = nothing
  IS_JOT_GENERATED::String = "false" # set to "true" in get_labels(res)
  user_defined_labels::AbstractDict{String, String} = Dict{String, String}()
end
StructTypes.StructType(::Type{Labels}) = StructTypes.Mutable()

function filter_user_defined_only(l::Labels)::Dict
  l.user_defined_labels
end

function Labels(d::AbstractDict{String, String})::Labels
  field_names = map(String, fieldnames(Labels))
  kwargs = Dict{Symbol, Any}()
  user_defined_labels = Dict{String, String}()
  for (k, v) in d
    if k in field_names
      kwargs[Symbol(k)] = v
    else
      user_defined_labels[k] = v
    end
  end
  kwargs[:user_defined_labels] = user_defined_labels
  Labels(; kwargs...)
end

function copy(l::Labels)::Labels
  flds = Dict(Symbol(fn) => getfield(l, fn) for fn in fieldnames(Labels))
  Labels(;flds...)
end

function get_responder_full_function_name(labels::Labels)::String
  labels.RESPONDER_PACKAGE_NAME * "." * labels.RESPONDER_FUNCTION_NAME
end

function to_aws_shorthand(l::Labels)::String
  normal_tags = join(
    ["$k=$(getfield(l, k))" for k in fieldnames(Labels) if k != :user_defined_labels],
    ",")
  if length(l.user_defined_labels) > 0
    addnl_tags = join(["$k=$v" for (k,v) in l.user_defined_labels], ",")
    normal_tags * "," * addnl_tags
  else
    normal_tags
  end
end

function to_json(l::Labels)::String
  normal_tags = [
    OrderedDict("Key" => String(k), "Value" => (isnothing(getfield(l, k)) ? "" : getfield(l, k)))
    for k in fieldnames(Labels) if k != :user_defined_labels
  ]
  addnl_tags = [OrderedDict("Key" => k, "Value" => v) for (k, v) in l.user_defined_labels]
  JSON3.write([ normal_tags ; addnl_tags ])
end

function to_docker_buildfile_format(l::Labels)::String
  normal_labels = join(
    ["$(String(k))=$(getfield(l, k))" for k in fieldnames(Labels) if k != :user_defined_labels],
    " ",
  )
  user_defined_labels = join(["$k=$v" for (k, v) in l.user_defined_labels], " ")
  normal_labels * " " * user_defined_labels
end

function add_user_defined_labels(l::Labels, new_tags::AbstractDict{String, String})::Labels
  c_l = copy(l)
  c_l.user_defined_labels = merge(c_l.user_defined_labels, new_tags)
  c_l
end


text:
function foreach_api_candidate(f, m::Module)
    exports = Set{Symbol}()
    for name in names(m)
        api = API(m, name; exported = true)
        f(api)
        push!(exports, name)
    end
    for name in names(m; all = true)
        name in exports && continue
        api = API(m, name; exported = false)
        f(api)
    end
end

function DocumentationOverview.find(m::Module; include = hasdoc)
    apis = API[]
    foreach_api_candidate(m) do api
        if include(api)
            push!(apis, api)
        end
        return
    end
    return apis
end


text:
module OnlineAdLearning

import JuMP
import MathOptInterface
const MOI = MathOptInterface

using LinearAlgebra
using Printf

using Distributions, StatsBase

export DataPoint, DataSet, Learner, TrainingResults

export Clairvoyant, RandomGuessing, ZeroReservePrice, ConstantPolicy, MIP, LPRelaxation, GradientAscent, DifferenceOfConvex

export L1Regularizer, L2Regularizer

include("data.jl")
include("learner.jl")
include("clairvoyant.jl")
include("random_guessing.jl")
include("constant_policy.jl")
include("mip.jl")
include("lp_relaxation.jl")
include("gradient_ascent.jl")
include("difference_of_convex.jl")
include("formulations.jl")
include("util.jl")

end # module

max_stars_repo_name: joehuchette/reserve-price-optimization

max_issues_repo_name: joehuchette/reserve-price-optimization

max_forks_repo_name: joehuchette/reserve-price-optimization


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

JLLWrappers.@generate_main_file_header("opusfile")
JLLWrappers.@generate_main_file("opusfile", UUID("4fb77979-f048-53ce-ae57-a45249826a33"))
end  # module opusfile_jll


text:
using Celeste
using Celeste.Model
using Celeste.SensitiveFloats

import Celeste.PSF: get_psf_at_point,
       get_psf_transform, initialize_psf_params, transform_psf_params!,
       unwrap_psf_params, wrap_psf_params,
       unconstrain_psf_params, constrain_psf_params,
       transform_psf_sensitive_float!,
       PsfOptimizer, fit_raw_psf_for_celeste,
       BivariateNormalDerivatives

using ForwardDiff

"""
Evaluate the sum of squared difference between the raw_psf and the psf
represented by psf_params, as well as the derivatives and hessians with
respect to unconstrained parameters.

Returns:
  - A sensitive float for the sum of squared differences.
"""
function evaluate_psf_fit{NumType <: Number}(
    psf_params::Vector{Vector{NumType}}, raw_psf::Matrix{Float64},
    calculate_derivs::Bool)

  K = length(psf_params)
  x_mat = PSF.get_x_matrix_from_psf(raw_psf);

  # TODO: allocate these outside?
  bvn_derivs = BivariateNormalDerivatives{NumType}(NumType);
  log_pdf = SensitiveFloats.zero_sensitive_float(PsfParams, NumType, 1);
  pdf = SensitiveFloats.zero_sensitive_float(PsfParams, NumType, 1);

  pixel_value = SensitiveFloats.zero_sensitive_float(PsfParams, NumType, K);
  squared_error = SensitiveFloats.zero_sensitive_float(PsfParams, NumType, K);

  PSF.evaluate_psf_fit!(
      psf_params, raw_psf, x_mat, bvn_derivs,
      log_pdf, pdf, pixel_value, squared_error, calculate_derivs)

  squared_error
end

function load_raw_psf(; x::Float64=500., y::Float64=500.)
  run_num = 3900
  camcol_num = 6
  field_num = 269
  b = 3

  psf_filename =
    @sprintf("%s/%s/%s/%s/psField-%06d-%d-%04d.fit",
        datadir, run_num, camcol_num, field_num,
                 run_num, camcol_num, field_num)
  psf_fits = FITSIO.FITS(psf_filename);
  raw_psf_comp = SDSSIO.read_psf(psf_fits, band_letters[b]);
  close(psf_fits)

  raw_psf_comp(x, y);
end

function test_transform_psf_params()
  K = 2
  psf_params = initialize_psf_params(K, for_test=true);
  psf_params_original = deepcopy(psf_params);
  psf_params_free = deepcopy(psf_params);
  psf_transform = PSF.get_psf_transform(psf_params);

  transform_psf_params!(psf_params, psf_params_free, psf_transform, true);
  transform_psf_params!(psf_params, psf_params_free, psf_transform, false);

  psf_params_free_2 = unconstrain_psf_params(psf_params, psf_transform)
  psf_params_2 = constrain_psf_params(psf_params_free, psf_transform)

  for k=1:K
    @test_approx_eq psf_params[k] psf_params_original[k]
    @test_approx_eq psf_params_free[k] psf_params_free_2[k]
    @test_approx_eq psf_params[k] psf_params_2[k]
  end
end

function test_psf_fit()
  raw_psf = load_raw_psf();

  # Initialize params
  K = 2
  psf_params = initialize_psf_params(K, for_test=true)
  psf_param_vec = wrap_psf_params(psf_params)

  function pixel_value_wrapper_sf{NumType <: Number}(
      psf_param_vec::Vector{NumType}, calculate_derivs::Bool)

    local psf_params = unwrap_psf_params(psf_param_vec)
    bvn_derivs = BivariateNormalDerivatives{NumType}(NumType);
    log_pdf = zero_sensitive_float(PsfParams, NumType, 1);
    pdf = zero_sensitive_float(PsfParams, NumType, 1);

    local pixel_value = zero_sensitive_float(PsfParams, NumType, K);

    local sigma_vec, sig_sf_vec, bvn_vec
    sigma_vec, sig_sf_vec, bvn_vec = PSF.get_sigma_from_params(psf_params);

    # sigma_vec = Array(Matrix{NumType}, K);
    # sig_sf_vec = Array(GalaxySigmaDerivs{NumType}, K);
    #
    # for k = 1:K
    #   sigma_vec[k] = PSF.get_bvn_cov(psf_params[k][psf_ids.e_axis],
    #                                   psf_params[k][psf_ids.e_angle],
    #                                   psf_params[k][psf_ids.e_scale])
    #   sig_sf_vec[k] = GalaxySigmaDerivs(
    #     psf_params[k][psf_ids.e_angle],
    #     psf_params[k][psf_ids.e_axis],
    #     psf_params[k][psf_ids.e_scale], sigma_vec[k], calculate_tensor=calculate_derivs);
    #
    # end

    clear!(pixel_value)
    PSF.evaluate_psf_pixel_fit!(
        x, psf_params, sigma_vec, sig_sf_vec, bvn_vec,
        bvn_derivs, log_pdf, pdf, pixel_value, calculate_derivs)

    pixel_value
  end

  function pixel_value_wrapper_value{NumType <: Number}(psf_param_vec::Vector{NumType})
    pixel_value_wrapper_sf(psf_param_vec, false).v[1]
  end

  x = Float64[1.0, 2.0]

  sigma_vec = Array(Matrix{Float64}, K);
  for k = 1:K
    sigma_vec[k] = PSF.get_bvn_cov(psf_params[k][psf_ids.e_axis],
                                    psf_params[k][psf_ids.e_angle],
                                    psf_params[k][psf_ids.e_scale])
  end

  println("Testing single pixel value")
  psf_components = PsfComponent[
    PsfComponent(psf_params[k][psf_ids.weight], psf_params[k][psf_ids.mu], sigma_vec[k])
                  for k = 1:K ];

  psf_rendered = get_psf_at_point(psf_components, rows=[ x[1] ], cols=[ x[2] ])[1];
  @test_approx_eq psf_rendered pixel_value_wrapper_value(psf_param_vec)

  pixel_value = deepcopy(pixel_value_wrapper_sf(psf_param_vec, true));

  ad_grad = ForwardDiff.gradient(pixel_value_wrapper_value, psf_param_vec);
  ad_hess = ForwardDiff.hessian(pixel_value_wrapper_value, psf_param_vec);

  @test_approx_eq ad_grad pixel_value.d[:]
  @test_approx_eq ad_hess[:] pixel_value.h[:]

  # Test the whole least squares function.
  println("Testing psf least squares")

  # Fewer pixels for quick testing.  Also, ForwardDiff.hessian runs into strange
  # problems on the whole image.
  keep_pixels = 20:30

  function evaluate_psf_fit_wrapper_sf{NumType <: Number}(
        psf_param_vec::Vector{NumType}, calculate_derivs::Bool)
    local psf_params = unwrap_psf_params(psf_param_vec)
    local squared_error =
      evaluate_psf_fit(psf_params, raw_psf[keep_pixels, keep_pixels], calculate_derivs)
    squared_error
  end

  function evaluate_psf_fit_wrapper_value{NumType <: Number}(psf_param_vec::Vector{NumType})
    local squared_error = evaluate_psf_fit_wrapper_sf(psf_param_vec, false);
    squared_error.v[1]
  end

  squared_error = deepcopy(evaluate_psf_fit_wrapper_sf(psf_param_vec, true));

  ad_grad = ForwardDiff.gradient(evaluate_psf_fit_wrapper_value, psf_param_vec);
  ad_hess = ForwardDiff.hessian(evaluate_psf_fit_wrapper_value, psf_param_vec);

  @test_approx_eq ad_grad squared_error.d[:]
  @test_approx_eq ad_hess[:] squared_error.h[:]
end

function test_transform_psf_sensitive_float()
  raw_psf = load_raw_psf();

  K = 2
  psf_params = initialize_psf_params(K, for_test=true);
  psf_transform = PSF.get_psf_transform(psf_params);
  psf_params_free = unconstrain_psf_params(psf_params, psf_transform);
  psf_params_free_vec = wrap_psf_params(psf_params_free)[:];

  # Fewer pixels for quick testing.  Also, ForwardDiff.hessian runs into strange
  # problems on the whole image.
  keep_pixels = 20:30

  function psf_fit_for_optim{NumType <: Number}(
      psf_params_free_vec::Vector{NumType}, calculate_derivs::Bool)

    local sf_free = zero_sensitive_float(PsfParams, NumType, K);
    local psf_params_free = unwrap_psf_params(psf_params_free_vec)
    local psf_params = constrain_psf_params(psf_params_free, psf_transform)
    local sf = evaluate_psf_fit(
      psf_params, raw_psf[keep_pixels, keep_pixels], calculate_derivs);
    transform_psf_sensitive_float!(
      psf_params, psf_transform, sf, sf_free, calculate_derivs)

    sf_free
  end

  function psf_fit_for_optim_val{NumType <: Number}(
      psf_params_free_vec::Vector{NumType})

    psf_fit_for_optim(psf_params_free_vec, false).v[1]
  end

  sf_free = deepcopy(psf_fit_for_optim(psf_params_free_vec, true));

  expected_value =
    evaluate_psf_fit(psf_params, raw_psf[keep_pixels, keep_pixels], false).v[1];
  ad_grad = ForwardDiff.gradient(psf_fit_for_optim_val, psf_params_free_vec);
  ad_hess = ForwardDiff.hessian(psf_fit_for_optim_val, psf_params_free_vec);

  @test_approx_eq expected_value sf_free.v[1]
  @test_approx_eq sf_free.d[:] ad_grad
  @test_approx_eq sf_free.h ad_hess
end

function test_psf_optimizer()
  raw_psf = load_raw_psf();

  K = psf_K
  psf_params = initialize_psf_params(K, for_test=false);
  psf_transform = get_psf_transform(psf_params);
  psf_optimizer = PsfOptimizer(psf_transform, K);

  nm_result = psf_optimizer.fit_psf(raw_psf, psf_params)
  psf_params_fit =
    constrain_psf_params(unwrap_psf_params(nm_result.minimum), psf_transform)

  # Could this test be tighter?
  @test 0.0 < nm_result.f_minimum < 1e-3

  celeste_psf = fit_raw_psf_for_celeste(raw_psf)[1]
  rendered_psf = get_psf_at_point(celeste_psf);

  @test_approx_eq nm_result.f_minimum sum((raw_psf - rendered_psf) .^ 2)

  # Make sure that re-using the optimizer gets the same results.
  raw_psf_10_10 = load_raw_psf(x=10., y=10.);
  celeste_psf_10_10_v1, psf_params_10_10_v1 = fit_raw_psf_for_celeste(raw_psf_10_10)
  celeste_psf_10_10_v2, psf_params_10_10_v2 =
    fit_raw_psf_for_celeste(raw_psf_10_10, psf_optimizer, psf_params)
  for k=1:K
    @test_approx_eq psf_params_10_10_v1[k] psf_params_10_10_v2[k]
    for field in fieldnames(celeste_psf_10_10_v1[k])
      @test_approx_eq getfield(celeste_psf_10_10_v1[k], field) getfield(celeste_psf_10_10_v2[k], field)
    end
  end
end

test_transform_psf_sensitive_float()
test_transform_psf_params()
test_psf_fit()
test_psf_optimizer()