text:
# Trigger several successive GC sweeps. This is more comprehensive than running just a
# single sweep, since freeable objects may need more than one sweep to be appropriately
# marked and freed.
gcscrub() = (gc(); gc(); gc(); gc())

#############
# Benchmark #
#############

type Benchmark{id}
    params::Parameters
end

params(b::Benchmark) = b.params

function loadparams!(b::Benchmark, params::Parameters, fields...)
    loadparams!(b.params, params, fields...)
    return b
end

#############
# execution #
#############
# Note that trials executed via `run` and `lineartrial` are always executed at the top-level
# scope of the module returned by `current_module()`. This is to avoid any weird quirks
# when calling the experiment from within other contexts.

sample(b::Benchmark, args...) = error("no execution method defined on type $(typeof(b))")
_run(b::Benchmark, args...; kwargs...) = error("no execution method defined on type $(typeof(b))")

function Base.run(b::Benchmark, p::Parameters = b.params; kwargs...)
    return eval(current_module(), :(BenchmarkTools._run($(b), $(p); $(kwargs...))))
end

function Base.run(group::BenchmarkGroup, args...; verbose::Bool = false, pad = "", kwargs...)
    result = similar(group)
    gcscrub() # run GC before running group, even if individual benchmarks don't manually GC
    i = 1
    for id in keys(group)
        verbose && (println(pad, "($(i)/$(length(group))) benchmarking ", repr(id), "..."); tic())
        result[id] = run(group[id], args...; verbose = verbose, pad = pad*"  ", kwargs...)
        verbose && (println(pad, "done (took ", toq(), " seconds)"); i += 1)
    end
    return result
end

function _lineartrial(b::Benchmark, p::Parameters = b.params; maxevals = RESOLUTION, kwargs...)
    params = Parameters(p; kwargs...)
    estimates = zeros(maxevals)
    completed = 0
    params.gctrial && gcscrub()
    start_time = time()
    for evals in eachindex(estimates)
        params.gcsample && gcscrub()
        params.evals = evals
        estimates[evals] = first(sample(b, params))
        completed += 1
        ((time() - start_time) > params.seconds) && break
    end
    return estimates[1:completed]
end

function lineartrial(b::Benchmark, p::Parameters = b.params; kwargs...)
    return eval(current_module(), :(BenchmarkTools._lineartrial($(b), $(p); $(kwargs...))))
end

warmup(item, verbose = true) = run(item; verbose = verbose, samples = 1, evals = 1,
                                   gctrial = false, gcsample = false)

####################
# parameter tuning #
####################

# The tuning process is as follows:
#
#   1. Using `lineartrial`, take one sample of the benchmark for each `evals` in `1:RESOLUTION`.
#
#   2. Extract the minimum sample found in this trial. Hopefully, this value will be
#      reasonably close to the the true benchmark time. At the very least, we can be certain
#      that the minimum sample overcomes discretization noise, assuming a sufficient time
#      budget for the trial (which technically should never need to be more than a couple of
#      seconds if `RESOLUTION` ∼ 1μs).
#
#   3. Using the sample from step 2 as a reasonable estimate `t`, calculate `evals`. This
#      is done by indexing into a pregenerated table that maps potential `t` values to the
#      appropriate `evals` values. The actual generation of this table is accomplished via
#      heuristic, driven mainly by empircal insight. It's built on the following postulates:
#          - Small variations in `t` should not cause large variations in `evals`
#          - `RESOLUTION` <= 1000ns
#          - `evals` <= `RESOLUTION` if `RESOLUTION` is in nanoseconds and `t` >= 1ns
#          - Increasing `t` should generally decrease `evals`
#
#   4. By default, tune the `samples` parameter as a function of `t`, `evals`, and the time
#      budget for the benchmark.

# the logistic function is useful for determining `evals` for `1 < t < RESOLUTION`
logistic(u, l, k, t, t0) = round(Int, ((u - l) / (1 + exp(-k * (t - t0)))) + l)

const EVALS = Vector{Int}(9000) # any `t > length(EVALS)` should get an `evals` of 1
for t in 1:400    (EVALS[t] = logistic(1006, 195, -0.025, t, 200)) end # EVALS[1] == 1000, EVALS[400] == 200
for t in 401:1000 (EVALS[t] = logistic(204, -16, -0.01, t, 800))   end # EVALS[401] == 200, EVALS[1000] == 10
for i in 1:8      (EVALS[((i*1000)+1):((i+1)*1000)] = 11 - i)      end # linearly decrease from EVALS[1000]

guessevals(t) = t <= length(EVALS) ? EVALS[t] : 1

function tune!(group::BenchmarkGroup; verbose::Bool = false, pad = "", kwargs...)
    gcscrub() # run GC before running group, even if individual benchmarks don't manually GC
    i = 1
    for id in keys(group)
        verbose && (println(pad, "($(i)/$(length(group))) tuning ", repr(id), "..."); tic())
        tune!(group[id]; verbose = verbose, pad = pad*"  ", kwargs...)
        verbose && (println(pad, "done (took ", toq(), " seconds)"); i += 1)
    end
    return group
end

function tune!(b::Benchmark, p::Parameters = b.params;
               verbose::Bool = false, pad = "", kwargs...)
    warmup(b, false)
    estimate = ceil(Int, minimum(lineartrial(b, p; kwargs...)))
    b.params.evals = guessevals(estimate)
    return b
end

#############################
# @benchmark/@benchmarkable #
#############################

function prunekwargs(args...)
    firstarg = first(args)
    if isa(firstarg, Expr) && firstarg.head == :parameters
        return prunekwargs(drop(args, 1)..., firstarg.args...)
    else
        core = firstarg
        params = collect(drop(args, 1))
        for ex in params
            if isa(ex, Expr) && ex.head == :(=)
                ex.head = :kw
            end
        end
        if isa(core, Expr) && core.head == :kw
            core.head = :(=)
        end
        return core, params
    end
end

function hasevals(params)
    for p in params
        if isa(p, Expr) && p.head == :kw && first(p.args) == :evals
            return true
        end
    end
    return false
end

function collectvars(ex::Expr, vars::Vector{Symbol} = Symbol[])
    if ex.head == :(=)
        lhs = first(ex.args)
        if isa(lhs, Symbol)
            push!(vars, lhs)
        elseif isa(lhs, Expr) && lhs.head == :tuple
            append!(vars, lhs.args)
        end
    elseif (ex.head == :comprehension || ex.head == :generator)
        arg = ex.args[1]
        isa(arg, Expr) && collectvars(arg, vars)
    else
        for arg in ex.args
            isa(arg, Expr) && collectvars(arg, vars)
        end
    end
    return vars
end

macro benchmark(args...)
    tmp = gensym()
    _, params = prunekwargs(args...)
    tune_expr = hasevals(params) ? :() : :(BenchmarkTools.tune!($(tmp)))
    return esc(quote
        $(tmp) = BenchmarkTools.@benchmarkable $(args...)
        BenchmarkTools.warmup($(tmp))
        $(tune_expr)
        BenchmarkTools.Base.run($(tmp))
    end)
end

macro benchmarkable(args...)
    core, params = prunekwargs(args...)

    # extract setup/teardown if present, removing them from the original expression
    setup, teardown = :(), :()
    delinds = Int[]
    for i in eachindex(params)
        ex = params[i]
        if ex.args[1] == :setup
            setup = ex.args[2]
            push!(delinds, i)
        elseif ex.args[1] == :teardown
            teardown = ex.args[2]
            push!(delinds, i)
        end
    end
    deleteat!(params, delinds)

    # extract any variable bindings shared between the core and setup expressions
    setup_vars = collectvars(setup)
    core_vars = isa(core, Expr) ? collectvars(core) : []
    out_vars = filter(var -> var in setup_vars, core_vars)

    # generate the benchmark definition
    return esc(quote
        BenchmarkTools.generate_benchmark_definition(current_module(),
                                                     $(Expr(:quote, out_vars)),
                                                     $(Expr(:quote, setup_vars)),
                                                     $(Expr(:quote, core)),
                                                     $(Expr(:quote, setup)),
                                                     $(Expr(:quote, teardown)),
                                                     BenchmarkTools.Parameters($(params...)))
    end)
end

# `eval` an expression that forcibly defines the specified benchmark at
# top-level of `eval_module`, hopefully ensuring that the "location" of the benchmark's
# definition will not be a factor for the sake of performance testing.
#
# The double-underscore-prefixed variable names are not particularly hygienic - it's
# possible for them to conflict with names used in the setup or teardown expressions.
# A more robust solution would be preferable.
function generate_benchmark_definition(eval_module, out_vars, setup_vars,
                                       core, setup, teardown, params)
    id = Expr(:quote, gensym("benchmark"))
    corefunc = gensym("core")
    samplefunc = gensym("sample")
    signature = Expr(:call, corefunc, setup_vars...)
    if length(out_vars) == 0
        returns = :()
        invocation = signature
        core_body = core
    elseif length(out_vars) == 1
        returns = :(return $(out_vars[1]))
        invocation = :($(out_vars[1]) = $(signature))
        core_body = :($(core); $(returns))
    else
        returns = :(return $(Expr(:tuple, out_vars...)))
        invocation = :($(Expr(:tuple, out_vars...)) = $(signature))
        core_body = :($(core); $(returns))
    end
    eval(eval_module, quote
        @noinline $(signature) = begin $(core_body) end
        @noinline function $(samplefunc)(__params::BenchmarkTools.Parameters)
            $(setup)
            __evals = __params.evals
            __gc_start = Base.gc_num()
            __start_time = time_ns()
            for __iter in 1:__evals
                $(invocation)
            end
            __sample_time = time_ns() - __start_time
            __gcdiff = Base.GC_Diff(Base.gc_num(), __gc_start)
            $(teardown)
            __time = max((__sample_time / __evals) - __params.overhead, 0.001)
            __gctime = max((__gcdiff.total_time / __evals) - __params.overhead, 0.0)
            __memory = Int(fld(__gcdiff.allocd, __evals))
            __allocs = Int(fld(__gcdiff.malloc + __gcdiff.realloc +
                               __gcdiff.poolalloc + __gcdiff.bigalloc,
                               __evals))
            return __time, __gctime, __memory, __allocs
        end
        function BenchmarkTools.sample(b::BenchmarkTools.Benchmark{$(id)},
                                       p::BenchmarkTools.Parameters = b.params)
            return $(samplefunc)(p)
        end
        function BenchmarkTools._run(b::BenchmarkTools.Benchmark{$(id)},
                                     p::BenchmarkTools.Parameters;
                                     verbose = false, pad = "", kwargs...)
            params = BenchmarkTools.Parameters(p; kwargs...)
            @assert params.seconds > 0.0 "time limit must be greater than 0.0"
            params.gctrial && BenchmarkTools.gcscrub()
            start_time = time()
            trial = BenchmarkTools.Trial(params)
            iters = 1
            while (time() - start_time) < params.seconds
                params.gcsample && BenchmarkTools.gcscrub()
                push!(trial, $(samplefunc)(params)...)
                iters += 1
                iters > params.samples && break
            end
            return sort!(trial)
        end
        BenchmarkTools.Benchmark{$(id)}($(params))
    end)
end


text:
#=
Mod:
- Julia version: 1.0
- Author: singularitti
- Date: 2019-07-28
=#
module Mod

using Test

using Test

using Fortran90Namelists.FortranToJulia

# Test data referenced from: http://www-classes.usc.edu/engr/ce/108/text/fbk01.htm.

@testset "Test `parse*` functions" begin
    @test parse(Float32, f"1.0E-6") === 1.0f-6
    @test parse(Float32, f"3.2767e+2") === 3.2767f2
    @test parse(Float32, f"1.89e-14") === 1.89f-14
    @test parse(Float32, f"-0.65e-2") === -0.65f-2
    @test parse(Float32, f"+1e8") === 1f8

    @test parse(Float64, f"1.0D-6") === 1.0e-6
    @test parse(Float64, f"1d-8") === 1e-8
    @test parse(Float64, f"0.") === 0.
    @test parse(Float64, f"1.00") === 1.00
    @test parse(Float64, f"+3.141593") === 3.141593
    @test parse(Float64, f"+3.1415926535d+0") === +3.1415926535
    @test parse(Float64, f"-4.78d+6") === -4.78e6
    @test parse(Float64, f"1.0d+0") === 1.0

    @test parse(Complex{Float64}, f"(5.229, -4.78)") === Complex(5.229, -4.78)
    @test parse(Complex{Float64}, f"(0.0,1.0)") === Complex(0.0, 1.0)
    @test parse(Complex{Float64}, f"(0.0,1)") === Complex(0.0, 1)
    @test parse(Complex{Float32}, f"(3.2767e+2, -0.65e-2)") === Complex{Float32}(3.2767e2, -0.65e-2)

    @test parse(Int, f"124") === 124
    @test parse(Int, f"-448") === -448
    @test parse(Int, f"0") === 0
    @test parse(Int, f"32767") === 32767
    @test parse(Int, f"2147483647") === 2147483647
    @test parse(Int, f"-9874") === -9874

    @test parse(Bool, f".true.") === true
    @test parse(Bool, f".t.") === true
    @test parse(Bool, f".false.") === false
    @test parse(Bool, f".f.") === false

    @test parse(String, f"''") === ""
    @test parse(String, f"\"./tmp234\"") === "./tmp234"
    @test parse(String, f"'david'") === "david"
end  # testset

end

ext: jl


text:
using PlanktonIndividuals, Serialization

grid = RegularRectilinearGrid(size = (1, 1, 16), spacing = (32, 32, 2),
                              topology = (Bounded, Bounded, Bounded), halo = (2, 2, 2))

model = PlanktonModel(CPU(), grid)

TP = sum((interior(model.nutrients.PO4.data, grid) .+
          interior(model.nutrients.DOP.data, grid) .+
          interior(model.nutrients.POP.data, grid)) .* grid.Δx .* grid.Δy .* grid.Δz)
TP = TP + sum(model.individuals.phytos.sp1.data.Pq .+
              model.individuals.phytos.sp1.data.Bm .* model.individuals.phytos.sp1.p.R_PC)

uvel = zeros(2,1,16,11)
vvel = zeros(1,2,16,11)
wvel = zeros(1,1,17,11)

for i in 1:11
    uvel[:,:,:,i] .= randn(2,1,16) .* 1e-4
    vvel[:,:,:,i] .= randn(1,2,16) .* 1e-4
    wvel[:,:,:,i] .= randn(1,1,17) .* 1e-4
end

sim = PlanktonSimulation(model, ΔT = 60, iterations = 10, vels=(u=uvel, v=vvel, w=wvel))

update!(sim)

TPt = sum((interior(model.nutrients.PO4.data, grid) .+
          interior(model.nutrients.DOP.data, grid) .+
          interior(model.nutrients.POP.data, grid)) .* grid.Δx .* grid.Δy .* grid.Δz)
TPt = TPt + sum(model.individuals.phytos.sp1.data.Pq .+
                model.individuals.phytos.sp1.data.Bm .* model.individuals.phytos.sp1.p.R_PC)


text:
using LightGraphs, SparseArrays, LinearAlgebra

import Base: copy
import LightGraphs: nv, has_edge, add_edge!, rem_edge!, rem_vertex!,
    rem_vertices!, add_vertex!, add_vertices!, outneighbors, inneighbors, neighbors,
    vertices, adjacency_matrix, ne, is_directed, degree, indegree, outdegree, edges,
    has_vertex, all_neighbors

export Multigraph

mutable struct Multigraph{T<:Integer} <: AbstractMultigraph{T}
    adjlist::Dict{T, Vector{T}}
    _idmax::T
    function Multigraph{T}(d::Dict{T, Vector{T}}, _idmax::T) where {T<:Integer}
        adjlist = deepcopy(d)
        vs = keys(adjlist)
        for (v, l) in adjlist
            if l ⊆ vs
                sort!(l)
            else
                error("Some vertices connected to $v is not in the multigraph!")
            end
        end
        _idmax >= 0 || (_idmax = isempty(adjlist) ? 0 : maximum(vs))
        new{T}(adjlist, _idmax)
    end
end

Multigraph(adjlist::Dict{T, Vector{T}}) where {T<:Integer} = Multigraph{T}(adjlist, isempty(adjlist) ? 0 : maximum(keys(adjlist)))
function Multigraph(adjmx::AbstractMatrix{U}) where {U<:Integer}
    m, n = size(adjmx)
    if m != n
        error("Adjacency matrices should be square!")
    end
    if !issymmetric(adjmx)
        error("Adjacency matrices should be symmetric!")
    end
    if sum(adjmx .!= 0) != sum(adjmx .> 0)
        error("All elements in adjacency matrices should be non-negative!")
    end
    adjlist = Dict(zip((1:m), [Int[] for _ = 1:m]))
    for v1 = 1:m
        for v2 = 1:m
            for i = 1:adjmx[v1, v2]
                push!(adjlist[v1], v2)
            end
        end
    end
    Multigraph{Int}(adjlist, m)
end
function Multigraph(n::T) where {T<:Integer}
    n >= 0 || error("Number of vertices should be non-negative")
    adjlist = Dict{T, Vector{T}}()
    for i = 1:n
        adjlist[i] = T[]
    end
    return Multigraph(adjlist)
end
Multigraph(g::SimpleGraph{T}) where {T<:Integer} = Multigraph(Dict(zip(T(1):nv(g), LightGraphs.SimpleGraphs.fadj(g))))

copy(mg::Multigraph{T}) where {T} = Multigraph{T}(deepcopy(mg.adjlist), mg._idmax)

nv(mg::Multigraph{T}) where {T<:Integer} = T(length(mg.adjlist))
vertices(mg::Multigraph) = collect(keys(mg.adjlist))
has_vertex(mg::Multigraph, v::Integer) = haskey(mg.adjlist, v)

function adjacency_matrix(mg::Multigraph)
    adjmx = spzeros(Int, nv(mg), nv(mg))

    ids = sort!(vertices(mg))
    for id1 in ids
        v1 = searchsortedfirst(ids, id1)
        for id2 in mg.adjlist[id1]
            v2 = searchsortedfirst(ids, id2)
            @inbounds adjmx[v1, v2] += 1
        end
    end
    return adjmx
end

function add_edge!(mg::Multigraph, me::AbstractMultipleEdge)
    s = src(me)
    d = dst(me)
    m = mul(me)
    if has_vertex(mg, s) && has_vertex(mg, d)
        for i = 1:m
            insert!(mg.adjlist[s], searchsortedfirst(mg.adjlist[s], d), d)
            insert!(mg.adjlist[d], searchsortedfirst(mg.adjlist[d], s), s)
        end
        return true
    end
    return false
end

function rem_edge!(mg::Multigraph, me::AbstractMultipleEdge)
    if has_edge(mg, me)
        s = src(me)
        d = dst(me)
        m = mul(me)
        for i = 1:m
            deleteat!(mg.adjlist[s], searchsortedfirst(mg.adjlist[s], d))
            deleteat!(mg.adjlist[d], searchsortedfirst(mg.adjlist[d], s))
        end
        return true
    else
        return false
    end
end

function rem_vertices!(mg::Multigraph{T}, vs::Vector{T}) where {T<:Integer}
    if all(has_vertex(mg, v) for v in vs)
        for v in vs
            for u in neighbors(mg, v)
                l = mg.adjlist[u]
                deleteat!(l, searchsorted(l, v))
            end
            delete!(mg.adjlist, v)
        end
        if mg._idmax in vs
            mg._idmax = maximum(keys(mg.adjlist))
        end
        return true
    end
    return false
end

function add_vertices!(mg::Multigraph{T}, n::Integer) where {T<:Integer}
    idmax = mg._idmax
    mg._idmax += n
    new_ids = collect((idmax+1):(idmax+n))
    for i in new_ids
        mg.adjlist[i] = T[]
    end
    return new_ids
end

function outneighbors(mg::Multigraph, v::Integer; count_mul::Bool = false)
    has_vertex(mg, v) || error("Vertex not found!")
    if count_mul
        return copy(mg.adjlist[v])
    else
        return sort!(collect(Set(mg.adjlist[v])))
    end
end
neighbors(mg::Multigraph, v::Integer; count_mul::Bool = false) = outneighbors(mg, v, count_mul = count_mul)
inneighbors(mg::Multigraph, v::Integer; count_mul::Bool = false) = outneighbors(mg, v, count_mul = count_mul)
all_neighbors(mg::Multigraph, v::Integer) = outneighbors(mg, v)

function mul(mg::Multigraph, s::Integer, d::Integer)
    (has_vertex(mg, s) && has_vertex(mg, d)) || error("Vertices not found!")
    return length(searchsorted(mg.adjlist[s], d))
end

is_directed(mg::Multigraph) = false
function ne(mg::Multigraph; count_mul::Bool = false)
    if count_mul
        return sum([sum(mg.adjlist[v] .>= v) for v in vertices(mg)])
    else
        return sum([sum(Set(mg.adjlist[v]) .>= v) for v in vertices(mg)])
    end
end

degree(mg::Multigraph) = Dict(zip(vertices(mg), [length(mg.adjlist[v]) for v in vertices(mg)]))
indegree(mg::Multigraph) = degree(mg)
outdegree(mg::Multigraph) = degree(mg)
degree(mg::Multigraph, v::Integer) = length(mg.adjlist[v])
indegree(mg::Multigraph, v::Integer) = degree(mg, v)
outdegree(mg::Multigraph, v::Integer) = degree(mg, v)


text:
#JSON Parser
#Adapted from http://www.mathworks.com/matlabcentral/fileexchange/23393
#Original BSD Licence, (c) 2011, François Glineur

function parse_json(strng::String)
    pos = 1
    len = length(strng)
    # String delimiters and escape characters are identified beforehand to improve speed
    # esc = regexp(str, "[\"\\\\]"); index_esc = 1; len_esc = length(esc);  #TODO Enable for speed

    function parse_object()
        parse_char('{')
        object = Dict{String, Any}()
        if next_char() != '}'
            while true
                str = parse_string()
                if isempty(str)
                    error("Name of value at position $pos cannot be empty")
                end
                parse_char(':')
                val = parse_value()
                object[str] = val
                if next_char() == '}'
                    break
                end
                parse_char(',')
            end
        end
        parse_char('}')
        return object
    end

    function parse_array()
        parse_char('[')
        object = Set()
        if next_char != ']'
            while true
                val = parse_value()
                push!(object, val)
                if next_char() == ']'
                    break
                end
                parse_char(',')
            end
        end
        parse_char(']')
        return object
    end

    function parse_char(c::Char)
        skip_whitespace()
        if pos > len || strng[pos] != c
            error("Expected $c at position $pos")
        else
            pos = pos + 1
            skip_whitespace()
        end
    end

    function next_char()
        skip_whitespace()
        if pos > len
            c = '\0'
        else
            c = strng[pos]
        end
    end

    function skip_whitespace()
        while pos <= len && isspace(strng[pos])
            pos = pos + 1
        end
    end

    function parse_string()
        if strng[pos] != '"'
            error("String starting with quotation expected at position $pos")
        else
            pos = pos + 1
        end
        str = ""
        while pos <= len
            # while index_esc <= len_esc && esc(index_esc) < pos
            #     index_esc = index_esc + 1
            # end
            # if index_esc > len_esc
            #     str = string(str, strng[pos:end])
            #     pos = len + 1
            #     break
            # else
            #     str = string(str, strng[pos:esc(index_esc)-1])
            #     pos = esc(index_esc)
            # end
            nc = strng[pos]
            if nc == '"'
                pos = pos + 1
                return string(str)
            elseif nc ==  '\\'
                if pos+1 > len
                    error_pos("End of file reached right after escape character")
                end
                pos = pos + 1
                anc = strng[pos]
                if anc == '"' || anc == '\\' || anc == '/'
                    str = string(str, strng[pos])
                    pos = pos + 1
                elseif anc ==  'b' || anc == 'f'|| anc == 'n' || anc == 'r' || anc == 't'
                    str = string(str, '\\', string[pos])
                    pos = pos + 1
                elseif anc == 'u'
                    if pos+4 > len
                        error_pos("End of file reached in escaped unicode character")
                    end
                    str = string(str, strng[pos-1:pos+4])
                    pos = pos + 5
                end
            else # should never happen
                str = string(str,strng[pos])
                pos = pos + 1
            end
        end
        error("End of file while expecting end of string")
    end

    function parse_number()
        num_regex = r"^[\w]?[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?[\w]?"
        m = match(num_regex, strng[pos:min(len,pos+20)])
        if m === nothing
            error("Error reading number at position $pos")
        end
        delta = m.offset + length(m.match)
        pos = pos + delta -1
        return float(m.match)
    end

    function  parse_value()
        nc = strng[pos]
        if nc == '"'
            val = parse_string()
            return val
        elseif nc == '['
            val = parse_array()
            return val
        elseif nc == '{'
            val = parse_object()
            return val
        elseif nc == '-' || nc == '0' || nc == '1' || nc == '2' || nc == '3' || nc == '4' || nc == '5' || nc == '6' || nc == '7' || nc == '8' || nc == '9'
            val = parse_number()
            return val
        elseif nc == 't'
            if pos+3 <= len && strng[pos:pos+3] == "true"
                val = true
                pos = pos + 4
                return val
            end
        elseif nc == 'f'
            if pos+4 <= len && strng[pos:pos+4] == "false"
                val = false
                pos = pos + 5
                return val
            end
        elseif nc == 'n'
            if pos+3 <= len && strng[pos:pos+3] == "null"
                val = []
                pos = pos + 4
                return val
            end
        end
        error("Value expected at position $pos")
    end

    if pos <= len
        nc = next_char()
        if nc == '{'
            return parse_object()
        elseif nc ==  '['
            return parse_array()
        else
            error("Outer level structure must be an object or an array")
        end
    end
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
# test aspects specific to partial factors

using IncrementalInference
using Test

@testset "testing basic partial factor functions..." begin

fg = initfg()
addVariable!(fg, :x0, ContinuousScalar)

addFactor!(fg, [:x0], Prior(Normal()) )

fc = getFactor(fg, :x0f1)
@test !isPartial(fc)

fg = initfg()

addVariable!(fg, :x1, ContinuousMultivariate(2))

addFactor!(fg, [:x1], PartialPrior(Normal(), (1,)) )

@test isPartial(getFactor(fg, :x1f1))

end


text:
#!/usr/bin/env julia
"""Lempel-Ziv complexity for a binary sequence, in naive Julia code.

- How to use it? From Julia, it's easy:

julia> using LempelZiv
julia> s = "1001111011000010"
julia> LempelZiv.lempel_ziv_complexity(s)  # 1 / 0 / 01 / 1110 / 1100 / 0010
6

- Note: there is also a Python version, if you need.

- MIT Licensed, (C) 2017 Lilian Besson (Naereen)
  https://GitHub.com/Naereen/LempelZiv.jl
"""

__author__ = "Lilian Besson (Naereen)"
__version__ = "0.0.1"

using LempelZiv
using Base.Test

# write your own tests here
@test 6 == lempel_ziv_complexity("1001111011000010")


text:
using Test
using Plots

using GCMMesh
using GCMMesh.Mesh

FT = Float64

x1min, x1max = FT(0), FT(1)
x2min, x2max = FT(0), FT(1)

limits = (x1min, x2min, x1max, x2max)

@testset "simple rectangular grid opposing face" begin

    @testset "assert correct element numbering" begin
        mesh = equispaced_rectangular_mesh(limits..., 1, 1, (true, true))
        @test_throws AssertionError Mesh.opposing_face(mesh, 0, 1)
        @test_throws AssertionError Mesh.opposing_face(mesh, 2, 1)
    end

    @testset "1×1 element quad mesh with all periodic boundries" begin
        mesh = equispaced_rectangular_mesh(limits..., 1, 1, (true, true))
        @test Mesh.opposing_face(mesh, 1, 1) == (1, 3, true)
        @test Mesh.opposing_face(mesh, 1, 2) == (1, 4, true)
        @test Mesh.opposing_face(mesh, 1, 3) == (1, 1, true)
        @test Mesh.opposing_face(mesh, 1, 4) == (1, 2, true)
    end

    @testset "1×1 element quad mesh with 1 periodic boundary" begin
        mesh = equispaced_rectangular_mesh(limits..., 1, 1, (true, false))
        @test Mesh.opposing_face(mesh, 1, 1) == (0, 1, true)
        @test Mesh.opposing_face(mesh, 1, 2) == (1, 4, true)
        @test Mesh.opposing_face(mesh, 1, 3) == (0, 3, true)
        @test Mesh.opposing_face(mesh, 1, 4) == (1, 2, true)
    end

    @testset "1×1 element quad mesh with non-periodic boundaries" begin
        mesh = equispaced_rectangular_mesh(limits..., 1, 1, (false, false))
        @test Mesh.opposing_face(mesh, 1, 1) == (0, 1, true)
        @test Mesh.opposing_face(mesh, 1, 2) == (0, 2, true)
        @test Mesh.opposing_face(mesh, 1, 3) == (0, 3, true)
        @test Mesh.opposing_face(mesh, 1, 4) == (0, 4, true)
    end

    @testset "2×2 element quad mesh with non-periodic boundaries" begin
        mesh = equispaced_rectangular_mesh(limits..., 2, 2, (false, false))
        @test Mesh.opposing_face(mesh, 1, 1) == (0, 1, true)
        @test Mesh.opposing_face(mesh, 1, 2) == (2, 4, true)
        @test Mesh.opposing_face(mesh, 1, 3) == (3, 1, true)
        @test Mesh.opposing_face(mesh, 1, 4) == (0, 4, true)
        @test Mesh.opposing_face(mesh, 2, 1) == (0, 1, true)
        @test Mesh.opposing_face(mesh, 2, 2) == (0, 2, true)
        @test Mesh.opposing_face(mesh, 2, 3) == (4, 1, true)
        @test Mesh.opposing_face(mesh, 2, 4) == (1, 2, true)
    end
end

@testset "simple rectangular grid interior faces iterator" begin
    @testset "1×1 element quad mesh with all periodic boundries" begin
        mesh = equispaced_rectangular_mesh(limits..., 1, 1, (true, true))
        @test length(Mesh.interior_faces(mesh)) == 2
        faces = Mesh.interior_faces(mesh)
        @test tuple(mesh.face_neighbors[faces[1], :]...) == (1, 4, 1, 2, -1)
        @test tuple(mesh.face_neighbors[faces[2], :]...) == (1, 1, 1, 3, -1)
    end
    @testset "1×1 element quad mesh with 1 periodic boundary" begin
        mesh = equispaced_rectangular_mesh(limits..., 1, 1, (true, false))
        @test length(Mesh.interior_faces(mesh)) == 1
        faces = Mesh.interior_faces(mesh)
        @test tuple(mesh.face_neighbors[faces[1], :]...) == (1, 4, 1, 2, -1)
    end
    @testset "1×1 element quad mesh with non-periodic boundaries" begin
        mesh = equispaced_rectangular_mesh(limits..., 1, 1, (false, false))
        @test Mesh.interior_faces(mesh) == nothing
    end
    @testset "2×2 element quad mesh with non-periodic boundaries" begin
        mesh = equispaced_rectangular_mesh(limits..., 2, 2, (false, false))
        @test length(Mesh.interior_faces(mesh)) == 4
        faces = Mesh.interior_faces(mesh)
        @test tuple(mesh.face_neighbors[faces[1], :]...) == (1, 2, 2, 4, -1)
        @test tuple(mesh.face_neighbors[faces[2], :]...) == (3, 2, 4, 4, -1)
        @test tuple(mesh.face_neighbors[faces[3], :]...) == (1, 3, 3, 1, -1)
        @test tuple(mesh.face_neighbors[faces[4], :]...) == (2, 3, 4, 1, -1)
    end
end

@testset "simple rectangular grid boundary faces iterator" begin
    @testset "1×1 element quad mesh with all periodic boundries" begin
        mesh = equispaced_rectangular_mesh(limits..., 1, 1, (true, true))
        @test Mesh.boundary_faces(mesh, 1) == nothing
        @test Mesh.boundary_faces(mesh, 2) == nothing
        @test Mesh.boundary_faces(mesh, 3) == nothing
        @test Mesh.boundary_faces(mesh, 4) == nothing
    end
    @testset "1×1 element quad mesh with 1 periodic boundary" begin
        mesh = equispaced_rectangular_mesh(limits..., 1, 1, (true, false))
        @test Mesh.boundary_faces(mesh, 1) == nothing
        @test Mesh.boundary_faces(mesh, 2) == nothing
        faces3 = Mesh.boundary_faces(mesh, 3)
        faces4 = Mesh.boundary_faces(mesh, 4)
        @test length(faces3) == 1
        @test length(faces4) == 1
        @test [Mesh.boundary_face_info(mesh, face) for face in faces3] == [(1, 1)]
        @test [Mesh.boundary_face_info(mesh, face) for face in faces4] == [(1, 3)]
    end
    @testset "1×1 element quad mesh with non-periodic boundaries" begin
        mesh = equispaced_rectangular_mesh(limits..., 1, 1, (false, false))
        faces1 = Mesh.boundary_faces(mesh, 1)
        faces2 = Mesh.boundary_faces(mesh, 2)
        faces3 = Mesh.boundary_faces(mesh, 3)
        faces4 = Mesh.boundary_faces(mesh, 4)

        @test length(faces1) == 1
        @test length(faces2) == 1
        @test length(faces3) == 1
        @test length(faces4) == 1
        @test [Mesh.boundary_face_info(mesh, face) for face in faces1] == [(1, 4)]
        @test [Mesh.boundary_face_info(mesh, face) for face in faces2] == [(1, 2)]
        @test [Mesh.boundary_face_info(mesh, face) for face in faces3] == [(1, 1)]
        @test [Mesh.boundary_face_info(mesh, face) for face in faces4] == [(1, 3)]
        @test Mesh.boundary_tag(mesh, :west) == 1
        @test Mesh.boundary_tag(mesh, :east) == 2
        @test Mesh.boundary_tag(mesh, :south) == 3
        @test Mesh.boundary_tag(mesh, :north) == 4
        @test Mesh.boundary_tag(mesh, :northeast) == nothing # no such boundary for rectangle
    end
    @testset "2×3 element quad mesh with non-periodic boundaries" begin
        mesh = equispaced_rectangular_mesh(limits..., 2, 3, (false, false))
        faces1 = Mesh.boundary_faces(mesh, 1)
        faces2 = Mesh.boundary_faces(mesh, 2)
        faces3 = Mesh.boundary_faces(mesh, 3)
        faces4 = Mesh.boundary_faces(mesh, 4)
        @test length(faces1) == 3
        @test length(faces2) == 3
        @test length(faces3) == 2
        @test length(faces4) == 2

        @test [Mesh.boundary_face_info(mesh, face) for face in faces1] == [(1, 4), (3, 4), (5, 4)]
        @test [Mesh.boundary_face_info(mesh, face) for face in faces2] == [(2, 2), (4, 2), (6, 2)]
        @test [Mesh.boundary_face_info(mesh, face) for face in faces3] == [(1, 1), (2, 1)]
        @test [Mesh.boundary_face_info(mesh, face) for face in faces4] == [(5, 3), (6, 3)]
    end
end

@testset "simple rectangular grid vertex iterator" begin
    @testset "1×1 element quad mesh with all periodic boundries" begin
        mesh = equispaced_rectangular_mesh(limits..., 1, 1, (true, true))
        uverts = mesh.unique_verts
        @test length(uverts) == 1
        @test Mesh.unique_vertex_connectivity(mesh, 1) == [(1, 1), (1, 2), (1, 3), (1, 4)]
        @test length(Mesh.unique_vertex_connectivity(mesh, 1)) == 4
    end
    @testset "1×1 element quad mesh with 1 periodic boundary" begin
        mesh = equispaced_rectangular_mesh(limits..., 1, 1, (true, false))
        uverts = mesh.unique_verts
        @test length(uverts) == 2
        vc1 = Mesh.unique_vertex_connectivity(mesh, 1)
        vc2 = Mesh.unique_vertex_connectivity(mesh, 2)
        @test length(vc1) == 2
        @test vc1 == [(1, 1), (1, 2)]
        @test length(vc2) == 2
        @test vc2 == [(1, 3), (1, 4)]
    end
    @testset "1×1 element quad mesh with non-periodic boundaries" begin
        mesh = equispaced_rectangular_mesh(limits..., 1, 1, (false, false))
        uverts = mesh.unique_verts
        @test length(uverts) == 4
        vc1 = Mesh.unique_vertex_connectivity(mesh, 1)
        vc2 = Mesh.unique_vertex_connectivity(mesh, 2)
        vc3 = Mesh.unique_vertex_connectivity(mesh, 3)
        vc4 = Mesh.unique_vertex_connectivity(mesh, 4)
        @test length(vc1) == 1
        @test collect(vc1) == [(1, 1)]
        @test length(vc2) == 1
        @test collect(vc2) == [(1, 2)]
        @test length(vc3) == 1
        @test collect(vc3) == [(1, 4)]
        @test length(vc4) == 1
        @test collect(vc4) == [(1, 3)]
    end
    @testset "2×3 element quad mesh with non-periodic boundaries" begin
        mesh = equispaced_rectangular_mesh(limits..., 2, 3, (false, false))
        uverts = mesh.unique_verts
        @test length(uverts) == 3 * 4
        vc1 = Mesh.unique_vertex_connectivity(mesh, 1)
        @test length(vc1) == 1
        @test vc1 == [(1, 1)]
    end
    @testset "2×3 element quad mesh with periodic boundaries" begin
        mesh = equispaced_rectangular_mesh(limits..., 2, 3, (true, true))
        uverts = mesh.unique_verts
        @test length(uverts) == 2 * 3
        vc1 = Mesh.unique_vertex_connectivity(mesh, 1)
        vc2 = Mesh.unique_vertex_connectivity(mesh, 2)
        vc3 = Mesh.unique_vertex_connectivity(mesh, 3)
        vc4 = Mesh.unique_vertex_connectivity(mesh, 4)
        vc5 = Mesh.unique_vertex_connectivity(mesh, 5)
        vc6 = Mesh.unique_vertex_connectivity(mesh, 6)
        @test sort(vc1) == sort([(1, 1), (2, 2), (5, 4), (6, 3)])
        @test sort(vc2) == sort([(2, 1), (1, 2), (5, 3), (6, 4)])
        @test sort(vc3) == sort([(3, 1), (4, 2), (1, 4), (2, 3)])
        @test sort(vc4) == sort([(4, 1), (3, 2), (2, 4), (1, 3)])
        @test sort(vc5) == sort([(5, 1), (6, 2), (3, 4), (4, 3)])
        @test sort(vc6) == sort([(6, 1), (5, 2), (4, 4), (3, 3)])
    end
end

@testset "simple rectangular grid coordinates" begin
    @testset "1×1 element quad mesh with all periodic boundries" begin
        x1min, x1max = FT(0), FT(1)
        x2min, x2max = FT(0), FT(1)
        limits = (x1min, x2min, x1max, x2max)
        mesh = equispaced_rectangular_mesh(limits..., 1, 1, (true, true))
        c1, c2, c3, c4 = Mesh.vertex_coordinates(mesh, 1)
        @test c1 == (0.0, 0.0)
        @test c2 == (1.0, 0.0)
        @test c3 == (1.0, 1.0)
        @test c4 == (0.0, 1.0)

        x1min, x1max = FT(-1), FT(1)
        x2min, x2max = FT(-1), FT(1)
        limits = (x1min, x2min, x1max, x2max)
        mesh = equispaced_rectangular_mesh(limits..., 1, 1, (false, false))
        c1, c2, c3, c4 = Mesh.vertex_coordinates(mesh, 1)
        @test c1 == (-1.0, -1.0)
        @test c2 == (1.0, -1.0)
        @test c3 == (1.0, 1.0)
        @test c4 == (-1.0, 1.0)
    end
    @testset "2×4 element quad mesh with non-periodic boundaries" begin
        x1min, x1max = FT(0), FT(1)
        x2min, x2max = FT(0), FT(1)
        limits = (x1min, x2min, x1max, x2max)
        mesh = equispaced_rectangular_mesh(limits..., 2, 4, (false, false))
        c1, c2, c3, c4 = Mesh.vertex_coordinates(mesh, 1)
        @test c1 == (0.0, 0.0)
        @test c2 == (0.5, 0.0)
        @test c3 == (0.5, 0.25)
        @test c4 == (0.0, 0.25)

        c1, c2, c3, c4 = Mesh.vertex_coordinates(mesh, 8)
        @test c1 == (0.5, 0.75)
        @test c2 == (1.0, 0.75)
        @test c3 == (1.0, 1.0)
        @test c4 == (0.5, 1.0)
    end
    @testset "check coordinate type accuracy" begin
        x1min, x1max = big(0.0), big(1.0)
        x2min, x2max = big(0.0), big(1.0)
        limits = (x1min, x2min, x1max, x2max)
        mesh = equispaced_rectangular_mesh(limits..., 3, 1, (false, false))
        c1, c2, c3, c4 = Mesh.vertex_coordinates(mesh, 1)
        @test eltype(c2) == BigFloat
        @test c2[1] ≈ big(1.0) / big(3.0) rtol = eps(BigFloat)
        @test c2[2] == 0.0
    end
end

max_stars_repo_name: CliMA/GCMMesh

max_issues_repo_name: CliMA/GCMMesh

max_forks_repo_name: CliMA/GCMMesh


text:
using Documenter, DynamicDecisionNetworks

makedocs(;
    modules=[DynamicDecisionNetworks],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/zsunberg/DynamicDecisionNetworks.jl/blob/{commit}{path}#L{line}",
    sitename="DynamicDecisionNetworks.jl",
    authors="Zachary Sunberg <sunbergzach@gmail.com>",
    assets=String[],
)

deploydocs(;
    repo="github.com/zsunberg/DynamicDecisionNetworks.jl",
)


text:
module multimaterial_beam_examples
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using FinEtools.MeshExportModule

function multimaterial_beam_algo()
    println("""
    Multi-material beam. Rubber-like and metal-like halves,
    clamped, with shear traction at free end.
    """)
    E1 = 0.29e3;
    nu1 = 0.49;
    E2 = 0.4e4;
    nu2 = 0.3;
    W = 4.1;
    L = 12.;
    t = 6.5;
    nl = 2; nt = 1; nw = 1; ref = 9;
    p =   200.0/W/t;
    #  Loading in the Z direction
    loadv = [0;0;p]; dir = 3;
    tolerance  = t/1000;

    fens,fes  = H20block(L,W,t, nl*ref,nw*ref,nt*ref)

    # Clamped end of the beam
    l1  = selectnode(fens; box = [0 0 -Inf Inf -Inf Inf], inflate  =  tolerance)
    e1 = FDataDict("node_list"=>l1, "component"=>1, "displacement"=>0.0)
    e2 = FDataDict("node_list"=>l1, "component"=>2, "displacement"=>0.0)
    e3 = FDataDict("node_list"=>l1, "component"=>3, "displacement"=>0.0)

    # Traction on the opposite edge
    boundaryfes  =   meshboundary(fes);
    Toplist   = selectelem(fens,boundaryfes, box =  [L L -Inf Inf -Inf Inf], inflate =   tolerance);
    el1femm  =   FEMMBase(IntegDomain(subset(boundaryfes,Toplist), GaussRule(2, 2)))
    flux1 = FDataDict("femm"=>el1femm, "traction_vector"=>loadv)

    r1list   = selectelem(fens,fes, box =  [0 L/2. -Inf Inf -Inf Inf], inflate =   tolerance);
    r2list   = selectelem(fens,fes, box =  [L/2. L -Inf Inf -Inf Inf], inflate =   tolerance);

    # Model reduction type
    MR = DeforModelRed3D

    # Make region 1
    region1 = FDataDict("femm"=>FEMMDeforLinear(MR,
    IntegDomain(subset(fes,r1list), GaussRule(3,2)),
    MatDeforElastIso(MR, 0.0, E1, nu1, 0.0)))

    # Make region 2
    region2 = FDataDict("femm"=>FEMMDeforLinear(MR,
    IntegDomain(subset(fes,r2list), GaussRule(3,2)),
    MatDeforElastIso(MR, 0.0, E2, nu2, 0.0)))

    # Make model data
    modeldata =  FDataDict(
    "fens"=> fens, "regions"=>  [region1, region2],
    "essential_bcs"=>[e1, e2, e3], "traction_bcs"=>  [flux1])

    # Call the solver
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)
    geom = modeldata["geom"]
    u = modeldata["u"]

    # Write out mesh with displacements
    modeldata["postprocessing"] = FDataDict("file"=>"multimaterial_beam")
    modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)

    # Write out mesh with stresses
    modeldata["postprocessing"] = FDataDict("file"=>"multimaterial_beam_xy",
    "quantity"=> :Cauchy, "component"=> :xy)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)

    # Write out mesh with stresses
    modeldata["postprocessing"] = FDataDict("file"=>"multimaterial_beam_xz",
    "quantity"=> :Cauchy, "component"=> :xz)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)

    # Write out mesh with von Mises stresses
    modeldata["postprocessing"] = FDataDict("file"=>"multimaterial_beam_vm",
    "quantity"=> :vm)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)

    # Write out mesh with von Mises stresses, elementwise
    modeldata["postprocessing"] = FDataDict("file"=>"multimaterial_beam_vm-ew",
    "quantity"=> :vm)
    modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)

end # multimaterial_beam_algo

function allrun()
    println("#####################################################")
    println("# multimaterial_beam_algo ")
    multimaterial_beam_algo()
    return true
end # function allrun

end # module multimaterial_beam_examples


text:
@rule GammaMixture{N}(:out, Marginalisation) (q_switch::Any, q_a::NTuple{N, Any}, q_b::NTuple{N, GammaDistributionsFamily}) where { N } = begin
    πs = probvec(q_switch)
    return GammaShapeRate(sum(πs .* mean.(q_a)), sum(πs .* mean.(q_b)))
end


text:
abstract type EventLikelihood{T<:Real} <: AbstractLikelihood{T} end

include("poisson.jl")
include("negativebinomial.jl")

""" Return the labels in a vector of vectors for multiple outputs"""
function treat_labels!(y::AbstractVector{<:Int}, likelihood::EventLikelihood)
    return y, 1, likelihood
end

function treat_labels!(::AbstractVector{<:Real}, ::EventLikelihood)
    return error("For event count target(s) should be integers")
end


text:
__precompile__()

module DistributedStatistics

using Reexport
@reexport using DistributedArrays, StatsBase

# Extend support of fit(Histogram, vs)
function StatsBase.append!(h::AbstractHistogram{T,N}, vs::NTuple{N,DistributedArrays.DVector}) where {T,N}
    hs = asyncmap(procs(vs[1])) do p
        remotecall_fetch(p) do
            hl = append!(h, localpart.(vs))
        end
    end
    return merge!(hs...)
end
function StatsBase.append!(h::AbstractHistogram{T,N}, vs::NTuple{N,DistributedArrays.DVector}, wv::DistributedArrays.DVector) where {T,N}
    hs = asyncmap(procs(vs[1])) do p
        remotecall_fetch(p) do
            hl = append!(h, localpart.(vs), localpart(wv))
        end
    end
    return merge!(hs...)
end

end


text:
"""
This module will be deleted shortly, when we beginning using a rasterized PSF.
In the future, PSF code will appear only in `SDSSIO.jl` (to load the eigen-PSF)
and `model/psf.jl` (functions for computing parts of the log probability
relating to the PSF).
"""
module PSF

using Celeste
using Celeste: Const, @aliasscope, @unroll_loop
using ..Model
using ..Model: GalaxySigmaDerivs, BivariateNormalDerivatives, get_bvn_derivs!,
               eval_bvn_pdf!, transform_bvn_derivs!, BvnComponent, get_bvn_cov
using ..Transform
using ..SensitiveFloats.SensitiveFloat
import ..SensitiveFloats

import Optim
import WCS

using ForwardDiff
using StaticArrays

const ID_MAT_2D = eye(Float64, 2)

"""
A data type to store functions related to optimizing a PSF fit. Initialize
using the transform and number of components, and then call fit_psf
to perform a fit to a specificed psf and initial parameters.
"""
mutable struct PsfOptimizer
    psf_transform::DataTransform
    ftol::Float64
    grtol::Float64
    num_iters::Int
    verbose::Bool
    K::Int

    # Variable that will be allocated in optimization:
    bvn_derivs::BivariateNormalDerivatives{Float64}

    log_pdf::SensitiveFloat{Float64}
    pdf::SensitiveFloat{Float64}
    pixel_value::SensitiveFloat{Float64}
    squared_error::SensitiveFloat{Float64}
    sf_free::SensitiveFloat{Float64}

    psf_params_free_vec_cache::Vector{Float64}

    function PsfOptimizer(psf_transform::DataTransform, K::Int;
                          verbose::Bool=false,
                          ftol::Float64=1e-9,
                          grtol::Float64=1e-9)
        num_iters = 50

        bvn_derivs = BivariateNormalDerivatives{Float64}()

        log_pdf = SensitiveFloat{Float64}(length(PsfParams), 1, true, true)
        pdf = SensitiveFloat{Float64}(length(PsfParams), 1, true, true)
        pixel_value = SensitiveFloat{Float64}(length(PsfParams), K, true, true)
        squared_error = SensitiveFloat{Float64}(length(PsfParams), K, true, true)
        sf_free = SensitiveFloat{Float64}(length(PsfParams), K, true, true)

        psf_params_free_vec_cache = fill(NaN, K * length(PsfParams))

        new(psf_transform, ftol, grtol, num_iters, verbose, K,
            bvn_derivs, log_pdf, pdf, pixel_value, squared_error,
            sf_free, psf_params_free_vec_cache)
    end
end

############################################################################################

function psf_fit_for_optim(psf_optimizer::PsfOptimizer, raw_psf, x_mat, psf_params_free_vec::Vector)
    if psf_params_free_vec == psf_optimizer.psf_params_free_vec_cache
        return psf_optimizer.sf_free
    else
        copy!(psf_optimizer.psf_params_free_vec_cache, psf_params_free_vec)
    end

    psf_params_free = unwrap_psf_params(psf_params_free_vec)
    psf_params = constrain_psf_params(psf_params_free, psf_optimizer.psf_transform)

    # Update squared_error in place.
    evaluate_psf_fit!(psf_params, raw_psf, x_mat, psf_optimizer.bvn_derivs,
                      psf_optimizer.log_pdf, psf_optimizer.pdf,
                      psf_optimizer.pixel_value, psf_optimizer.squared_error, true)

    # Update sf_free in place.
    transform_psf_sensitive_float!(psf_params, psf_optimizer.psf_transform,
                                   psf_optimizer.squared_error,
                                   psf_optimizer.sf_free, true)

    psf_optimizer.sf_free
end

function fit_psf(psf_optimizer::PsfOptimizer, raw_psf::Matrix{Float64}, initial_params::Vector{Vector{Float64}})
    x_mat = get_x_matrix_from_psf(raw_psf)
    psf_params_free = unconstrain_psf_params(initial_params, psf_optimizer.psf_transform)
    psf_params_free_vec = vec(wrap_psf_params(psf_params_free))

    function psf_fit_value(psf_params_free_vec::Vector)
        psf_fit_for_optim(psf_optimizer, raw_psf, x_mat, psf_params_free_vec).v[]
    end

    function psf_fit_grad!(psf_params_free_vec::Vector{Float64}, grad::Vector{Float64})
        grad[:] = psf_fit_for_optim(psf_optimizer, raw_psf, x_mat, psf_params_free_vec).d[:]
    end

    function psf_fit_hess!(psf_params_free_vec::Vector{Float64}, hess::Matrix{Float64})
        hess[:] = psf_fit_for_optim(psf_optimizer, raw_psf, x_mat, psf_params_free_vec).h
        hess[:] = 0.5 * (hess + hess')
    end

    tr_method = Optim.NewtonTrustRegion(initial_delta=10.0, delta_hat=1e9, eta=0.1,
                                        rho_lower=0.2, rho_upper=0.75)

    options = Optim.Options(;
                            x_tol = 0.0, # Don't allow convergence in params
                            f_tol = psf_optimizer.ftol,
                            g_tol = psf_optimizer.grtol,
                            iterations = psf_optimizer.num_iters,
                            store_trace = psf_optimizer.verbose,
                            show_trace = false,
                            extended_trace = psf_optimizer.verbose)

    return Optim.optimize(psf_fit_value, psf_fit_grad!, psf_fit_hess!,
                          psf_params_free_vec, tr_method, options)
end

############################################################################################

"""
Return an image of a Celeste GMM PSF evaluated at rows, cols.

Args:
- psf_array: The PSF to be evaluated as an array of PsfComponent
- rows: The rows in the image (usually in pixel coordinates)
- cols: The column in the image (usually in pixel coordinates)

Returns:
 - The PSF values at rows and cols.    The default size is the same as
     that returned by get_psf_at_point applied to FITS header values.

Note that the point in the image at which the PSF is evaluated --
that is, the center of the image returned by this function -- is
already implicit in the value of psf_array.
"""
function get_psf_at_point(psf_array::Vector{PsfComponent};
                          rows = -25:25, cols = -25:25)
    function get_psf_value(psf::PsfComponent, row::Float64, col::Float64)
         x = @SVector(Float64[row, col]) - psf.xiBar
         exp_term = exp(-0.5 * dot(x, psf.tauBarInv * x) - 0.5 * psf.tauBarLd)
         psf.alphaBar * exp_term / (2 * pi)
    end

    return Float64[
     sum(get_psf_value(psf, float(row), float(col)) for psf in psf_array)
         for row in rows, col in cols ]
end

"""
Get the PSF located at a particular world location in an image.

Args:
 - world_loc: A location in world coordinates.
 - img: An Image

Returns:
 - An array of PsfComponent objects that represents the PSF as a mixture
     of Gaussians.
"""
function get_source_psf(world_loc::Vector{Float64}, img::Image, psf_K::Int)
    # Some stamps or simulated data have no raw psf information.    In that case,
    # just use the psf from the image.
    if size(img.raw_psf_comp.rrows) == (0, 0)
        # Also return a vector of empty psf params
        return img.psf, fill(fill(NaN, length(PsfParams)), psf_K)
    else
        pixel_loc = WCS.world_to_pix(img.wcs, world_loc)
        psfstamp = Model.eval_psf(img.raw_psf_comp, pixel_loc[1], pixel_loc[2])
        return PSF.fit_raw_psf_for_celeste(psfstamp, psf_K)
    end
end

"""
Get the PSF located at a particular world location in an image.

Args:
 - world_loc: A location in world coordinates.
 - img: An Image

Returns:
 - An array of PsfComponent objects that represents the PSF as a mixture
     of Gaussians.
"""
function get_source_psf(world_loc::Vector{Float64},
                        img::Image,
                        psf_optimizer::PsfOptimizer,
                        initial_psf_params::Vector{Vector{Float64}})
    # Some stamps or simulated data have no raw psf information. In that case,
    # just use the psf from the image.
    if size(img.raw_psf_comp.rrows) == (0, 0)
        return img.psf
    else
        pixel_loc = WCS.world_to_pix(img.wcs, world_loc)
        psfstamp = Model.eval_psf(img.raw_psf_comp, pixel_loc[1], pixel_loc[2])
        return PSF.fit_raw_psf_for_celeste(psfstamp, psf_optimizer, initial_psf_params)
    end
end

"""
Get initial values for a set of PSF parameters.

Args:
    - K: The number of components.
    - for_test: If true, return a set of slightly asymmetric parameters so that
                            derivatives are interesting for testing.

Returns:
    - A vector of PSF parameters, one for each component.
"""
function initialize_psf_params(K::Int; for_test::Bool=false)
    psf_params = Vector{Vector{Float64}}(K)
    for k=1:K
        if for_test
            # Choose asymmetric values for testing.
            psf_params[k] = zeros(length(PsfParams))
            psf_params[k][psf_ids.mu] = [0.1, 0.2]
            psf_params[k][psf_ids.gal_ab] = 0.8
            psf_params[k][psf_ids.gal_angle] = pi / 4
            psf_params[k][psf_ids.gal_scale] = sqrt(2 * k)
            psf_params[k][psf_ids.weight] = 1 / K + k / 10
        else
            psf_params[k] = zeros(length(PsfParams))
            psf_params[k][psf_ids.mu] = [0.0, 0.0]
            psf_params[k][psf_ids.gal_ab] = 0.95
            psf_params[k][psf_ids.gal_angle] = 0.0
            psf_params[k][psf_ids.gal_scale] = sqrt(2 * k)
            psf_params[k][psf_ids.weight] = 1 / K
        end
    end

    psf_params
end

"""
Given PSF parameters, return a DataTransform to transfrom them to and from
an unconstrained parameterization.

Args:
    - psf_params: The parameters to be transformed
    - scale: A vector as long as PsfParams to define a linear rescaling.

Returns:
    - A DataTransform object.
"""
function get_psf_transform(
        psf_params::Vector{Vector{Float64}};
        scale::Vector{Float64}=ones(length(PsfParams)))
    K = length(psf_params)
    bounds = Vector{ParamBounds}(length(psf_params))
    # Note that, for numerical reasons, the bounds must be on the scale
    # of reasonably meaningful changes.
    for k in 1:K
        bounds[k] = ParamBounds()
        bounds[k][:mu] = ParamBox[ ParamBox(-5.0, 5.0, scale[psf_ids.mu[1]]),
                                   ParamBox(-5.0, 5.0, scale[psf_ids.mu[2]]) ]
        bounds[k][:gal_ab] = ParamBox[ ParamBox(0.1, 1.0, scale[psf_ids.gal_ab] ) ]
        bounds[k][:gal_angle] =
            ParamBox[ ParamBox(-4 * pi, 4 * pi, scale[psf_ids.gal_angle] ) ]
        bounds[k][:gal_scale] =
            ParamBox[ ParamBox(0.05, 10.0, scale[psf_ids.gal_scale] ) ]

        # Note that the weights do not need to sum to one.
        bounds[k][:weight] = ParamBox[ ParamBox(0.05, 2.0, scale[psf_ids.weight] ) ]
    end
    DataTransform(bounds, collect(1:K), K)
end

"""
Transform psf_params to and from psf_params_free in place.

Args:
    - psf_params, psf_params_free: The constrained and unconstrained parameters,
            respectively.
    - psf_transform: The DataTransform to be applied.
    - to_unconstrained: If true, then update psf_params_free with the unconstrained
            representation of psf_params.    If false, update psf_params with the
            constrained version of psf_params_free.

Returns:
    - Updates psf_params or psf_params_free in place.
"""
function transform_psf_params!{NumType<:Number}(
        psf_params::Vector{Vector{NumType}}, psf_params_free::Vector{Vector{NumType}},
        psf_transform::DataTransform, to_unconstrained::Bool)
    for k=1:length(psf_params)
        for (param, constraint_vec) in psf_transform.bounds[k]
            for ind in 1:length(getfield(psf_ids, param))
                param_ind = getfield(psf_ids, param)[ind]
                constraint = constraint_vec[ind]
                to_unconstrained ?
                    psf_params_free[k][param_ind] =
                        Transform.unbox_parameter(psf_params[k][param_ind], constraint):
                    psf_params[k][param_ind] =
                        Transform.box_parameter(psf_params_free[k][param_ind], constraint)
            end
        end
    end

    true # return type
end

"""
Allocate memory for and return a constrained parameter set.
"""
function constrain_psf_params{NumType<:Number}(
        psf_params_free::Vector{Vector{NumType}}, psf_transform::DataTransform)
    K = length(psf_params_free)
    psf_params = Vector{Vector{NumType}}(K)
    for k=1:K
        psf_params[k] = zeros(NumType, length(PsfParams))
    end

    transform_psf_params!(psf_params, psf_params_free, psf_transform, false)

    psf_params
end

"""
Allocate memory for and return an unconstrained parameter set.
"""
function unconstrain_psf_params{NumType<:Number}(
        psf_params::Vector{Vector{NumType}}, psf_transform::DataTransform)
    K = length(psf_params)
    psf_params_free = Vector{Vector{NumType}}(K)
    for k=1:K
        psf_params_free[k] = zeros(NumType, length(PsfParams))
    end

    transform_psf_params!(psf_params, psf_params_free, psf_transform, true)

    psf_params_free
end

"""
Convert a single vector of psf parameters to a vector of vectors.
"""
function unwrap_psf_params{NumType <: Number}(psf_param_vec::Vector{NumType})
    @assert length(psf_param_vec) % length(PsfParams) == 0
    K = round(Int, length(psf_param_vec) / length(PsfParams))
    psf_param_mat = reshape(psf_param_vec, length(PsfParams), K)
    psf_params = Vector{Vector{NumType}}(K)
    for k = 1:K
        psf_params[k] = psf_param_mat[:, k]
    end
    psf_params
end

"""
Convert a vector of vectors of psf parameters to a single vector.
"""
function wrap_psf_params{NumType <: Number}(psf_params::Vector{Vector{NumType}})
    psf_params_mat = zeros(NumType, length(PsfParams), length(psf_params))
    for k=1:length(psf_params)
        psf_params_mat[:, k] = psf_params[k]
    end
    psf_params_mat[:]
end

"""
Return a sensitive float representing the value of the psf at pixel x
with all its associated derivatives (with respect to the constrained
parameterization).    Note that "fit" is a bad name -- it is just the pixel value.
TODO: fix the name.

Args:
    - x: The 2d location at which to evaluate the pdf
    - calculate_gradient: If false, only update the value of pixel_value
    - other values: Pre-allocated memory for intermediate calculations

Returns:
    - Updates pixel_value in place (and all the other placeholder values as well)
"""
function evaluate_psf_pixel_fit!{NumType <: Number}(
        x::SVector{2,Float64}, psf_params::Vector{Vector{NumType}},
        sig_sf_vec::Vector{GalaxySigmaDerivs{NumType}},
        bvn_vec::Vector{BvnComponent{NumType}},
        bvn_derivs::BivariateNormalDerivatives{NumType},
        log_pdf::SensitiveFloat{NumType},
        pdf::SensitiveFloat{NumType},
        pixel_value::SensitiveFloat{NumType},
        calculate_gradient::Bool)
    SensitiveFloats.zero!(pixel_value)

    K = length(psf_params)
    sigma_ids = (psf_ids.gal_ab, psf_ids.gal_angle, psf_ids.gal_scale)
    @inbounds for k = 1:K
        # I will put in the weights later so that the log pdf sensitive float
        # is accurate.
        bvn = bvn_vec[k]
        eval_bvn_pdf!(bvn_derivs, bvn, x)
        get_bvn_derivs!(bvn_derivs, bvn, true, true)
        transform_bvn_derivs!(bvn_derivs, sig_sf_vec[k], I, true)

        SensitiveFloats.zero!(log_pdf)
        SensitiveFloats.zero!(pdf)

        # This is redundant, but it's what eval_bvn_pdf returns.
        log_pdf.v[] = log(bvn_derivs.f_pre[1])

        if calculate_gradient
            for ind=1:2
                log_pdf.d[psf_ids.mu[ind]] = bvn_derivs.bvn_u_d[ind]
            end
            for ind=1:3
                log_pdf.d[sigma_ids[ind]] = bvn_derivs.bvn_s_d[ind]
            end
            log_pdf.d[psf_ids.weight] = 0

            for ind1 = 1:2, ind2 = 1:2
                log_pdf.h[psf_ids.mu[ind1], psf_ids.mu[ind2]] =
                    bvn_derivs.bvn_uu_h[ind1, ind2]
            end
            for mu_ind = 1:2, sig_ind = 1:3
                log_pdf.h[psf_ids.mu[mu_ind], sigma_ids[sig_ind]] =
                log_pdf.h[sigma_ids[sig_ind], psf_ids.mu[mu_ind]] =
                    bvn_derivs.bvn_us_h[mu_ind, sig_ind]
            end
            for ind1 = 1:3, ind2 = 1:3
                log_pdf.h[sigma_ids[ind1], sigma_ids[ind2]] =
                    bvn_derivs.bvn_ss_h[ind1, ind2]
            end
        end

        pdf_val = exp(log_pdf.v[])
        pdf.v[] = pdf_val

        if calculate_gradient
            for ind1 = 1:length(PsfParams)
                if ind1 == psf_ids.weight
                    pdf.d[ind1] = pdf_val
                else
                    pdf.d[ind1] = psf_params[k][psf_ids.weight] * pdf_val * log_pdf.d[ind1]
                end

                for ind2 = 1:ind1
                    pdf.h[ind1, ind2] = pdf.h[ind2, ind1] =
                        psf_params[k][psf_ids.weight] * pdf_val *
                        (log_pdf.h[ind1, ind2] + log_pdf.d[ind1] * log_pdf.d[ind2])
                end
            end

            # Weight hessian terms.
            for ind1 = 1:length(PsfParams)
                pdf.h[psf_ids.weight, ind1] = pdf.h[ind1, psf_ids.weight] =
                    pdf_val * log_pdf.d[ind1]
            end

        end

        pdf.v[] *= psf_params[k][psf_ids.weight]

        SensitiveFloats.add_sources_sf!(pixel_value, pdf, k)
    end

    true # Set return type
end

"""
Convert PSF parameters to covariance matrices and derivatives and BvnComponents.
"""
function get_sigma_from_params{NumType <: Number}(psf_params::Vector{Vector{NumType}})
    K = length(psf_params)
    sigma_vec = Vector{SMatrix{2,2,NumType,4}}(K)
    sig_sf_vec = Vector{GalaxySigmaDerivs{NumType}}(K)
    bvn_vec = Vector{BvnComponent{NumType}}(K)
    for k = 1:K
        sigma_vec[k] = get_bvn_cov(psf_params[k][psf_ids.gal_ab],
            psf_params[k][psf_ids.gal_angle],
            psf_params[k][psf_ids.gal_scale])
        sig_sf_vec[k] = GalaxySigmaDerivs(
            psf_params[k][psf_ids.gal_angle],
            psf_params[k][psf_ids.gal_ab],
            psf_params[k][psf_ids.gal_scale], sigma_vec[k], 1.0, true)

        bvn_vec[k] =
            BvnComponent(SVector{2,NumType}(psf_params[k][psf_ids.mu]), sigma_vec[k], 1.0)
    end
    sigma_vec, sig_sf_vec, bvn_vec
end

"""
evaluate_psf_fit but with pre-allocated memory for intermediate calculations.
"""
function evaluate_psf_fit!{NumType <: Number}(
        psf_params::Vector{Vector{NumType}}, raw_psf::Matrix{Float64},
        x_mat::Matrix{SVector{2,Float64}},
        bvn_derivs::BivariateNormalDerivatives{NumType},
        log_pdf::SensitiveFloat{NumType},
        pdf::SensitiveFloat{NumType},
        pixel_value::SensitiveFloat{NumType},
        squared_error::SensitiveFloat{NumType},
        calculate_gradient::Bool)
    K = length(psf_params)
    sigma_vec, sig_sf_vec, bvn_vec = get_sigma_from_params(psf_params)
    SensitiveFloats.zero!(squared_error)

    @inbounds for x_ind in 1:length(x_mat)
        SensitiveFloats.zero!(pixel_value)
        evaluate_psf_pixel_fit!(
                x_mat[x_ind], psf_params, sig_sf_vec, bvn_vec,
                bvn_derivs, log_pdf, pdf, pixel_value, calculate_gradient)

        diff = (pixel_value.v[] - raw_psf[x_ind])
        squared_error.v[] +=    diff ^ 2
        if calculate_gradient
            for ind1 = 1:length(squared_error.d)
                squared_error.d[ind1] += 2 * diff * pixel_value.d[ind1]
                for ind2 = 1:ind1
                    squared_error.h[ind1, ind2] +=
                        2 * (diff * pixel_value.h[ind1, ind2] +
                                 pixel_value.d[ind1] * pixel_value.d[ind2]')
                    squared_error.h[ind2, ind1] = squared_error.h[ind1, ind2]
                end
            end
        end # if calculate_gradient
    end

    squared_error
end

"""
Transform a sf, a SensitiveFloat that has derivatives with respect to
constrianed parameters, into sf_free, a SensitiveFloat that has derivatives
with respect to unconstrained parameteres.

Args:
    - psf_params: The constrianed parameters
    - psf_transform: The data transform to be applied
    - sf: The SensitiveFloat with derivatives with respect to the constrained params
    - sf_free: Updated in place.    The SensitiveFloat with derivatives with respect
                         to the unconstrained parameters.
    - calculate_gradient: If false, only calculate the value.

Returns:
    - Updates sf_free in place.
"""
function transform_psf_sensitive_float!{NumType<:Number}(
        psf_params::Vector{Vector{NumType}}, psf_transform::DataTransform,
        sf::SensitiveFloat{NumType}, sf_free::SensitiveFloat{NumType},
        calculate_gradient::Bool)
    sf_free.v[] = sf.v[]
    if calculate_gradient
        K = length(psf_params)

        # These are the hessians of each individual parameter's transform.    We
        # can represent it this way since each parameter's transform only depends on
        # its own value and not on others.

        # This is the diagonal of the Jacobian transform.
        jacobian_diag = zeros(length(PsfParams) * K)

        # These are hte Hessians of the each parameter's transform.
        hessian_values = zeros(length(PsfParams) * K)

        for k = 1:K
            offset = length(PsfParams) * (k - 1)
            for ind = 1:2
                mu_ind = psf_ids.mu[ind]
                jac, hess = Transform.box_derivatives(
                    psf_params[k][mu_ind], psf_transform.bounds[k][:mu][ind])
                jacobian_diag[offset + mu_ind] = jac
                hessian_values[offset + mu_ind] = hess
            end

            # The rest are one-dimensional.
            for field in setdiff(fieldnames(PsfParams), [ :mu ])
                ind = getfield(psf_ids, field)
                jac, hess = Transform.box_derivatives(
                    psf_params[k][ind], psf_transform.bounds[k][field][1])
                jacobian_diag[offset + ind] = jac
                hessian_values[offset + ind] = hess
            end
        end

        # Apply the transformations.
        for ind1 = 1:length(sf_free.d)
            sf_free.d[ind1] = jacobian_diag[ind1] * sf.d[ind1]

            for ind2 = 1:ind1
                # Calculate the Hessian
                sf_free.h[ind1, ind2] = sf_free.h[ind2, ind1] =
                    (jacobian_diag[ind1] * jacobian_diag[ind2]) * sf.h[ind1, ind2]
                if ind1 == ind2
                    sf_free.h[ind1, ind2] +=    hessian_values[ind1] * sf.d[ind1]
                end
            end
        end
    end

    true # return type
end

"""
Given a psf matrix, return a matrix of the same size, where each element of
the matrix is a 2-length vector of the [x1, x2] location of that matrix cell.
The locations are chosen so that the scale is pixels, but the center of the
psf is [0, 0].
"""
function get_x_matrix_from_psf(psf::Matrix{Float64})
    psf_center1, psf_center2 = (size(psf, 1) - 1) / 2 + 1, (size(psf, 2) - 1) / 2 + 1
    return [ @SVector Float64[i - psf_center1, j - psf_center2] for i = 1:size(psf, 1), j = 1:size(psf, 2) ]
end

"""
Fit a mixture of 2d Gaussians to a PSF image (evaluated at a single point).

Args:
    - raw_psf: An matrix image of the point spread function.
    - psf_optimizer: A PsfOptimizer object for the fit.

Returns:
    - A vector of PsfComponents fit to the raw_psf.
"""
function fit_raw_psf_for_celeste{P<:PsfOptimizer}(
        raw_psf::Array{Float64, 2}, psf_optimizer::P,
        initial_psf_params::Vector{Vector{Float64}})
    K = length(initial_psf_params)
    @assert K == psf_optimizer.K
    optim_result = fit_psf(psf_optimizer, raw_psf, initial_psf_params)
    psf_params_fit =
        constrain_psf_params(
            unwrap_psf_params(Optim.minimizer(optim_result)), psf_optimizer.psf_transform)

    sigma_vec = get_sigma_from_params(psf_params_fit)[1]
    celeste_psf = Vector{PsfComponent}(K)
    for k=1:K
        mu = psf_params_fit[k][psf_ids.mu]
        weight = psf_params_fit[k][psf_ids.weight]
        celeste_psf[k] = PsfComponent(weight, SVector{2,eltype(mu)}(mu), SMatrix{2,2,eltype(mu),4}(sigma_vec[k]))
    end

    celeste_psf, psf_params_fit
end

"""
Fit a mixture of 2d Gaussians to a PSF image (evaluated at a single point).

Args:
 - raw_psf: An matrix image of the point spread function.
 - K: The number of components to fit.

 Returns:
     - A vector of PsfComponents fit to the raw_psf.
"""
function fit_raw_psf_for_celeste(raw_psf::Array{Float64, 2}, K::Integer; ftol=1e-9)
    psf_params = PSF.initialize_psf_params(K, for_test=false)
    psf_transform = PSF.get_psf_transform(psf_params)
    psf_optimizer = PsfOptimizer(psf_transform, K)
    psf_optimizer.ftol = ftol
    fit_raw_psf_for_celeste(raw_psf, psf_optimizer, psf_params)
end

function trim_psf(raw_psf::Array{Float64, 2}; trim_percent=0.999)
    h_mid = cld(size(raw_psf, 1), 2)
    w_mid = cld(size(raw_psf, 2), 2)

    width = 1
    function get_trimmed_psf_image()
        raw_psf[(h_mid - width):(h_mid + width),
                (w_mid - width):(w_mid + width)]
    end

    psf_tot = sum(abs, raw_psf)
    while sum(abs, get_trimmed_psf_image()) < trim_percent * psf_tot
        width += 1
    end

    deepcopy(get_trimmed_psf_image())
end

end # module


text:
export  excelFluxReport

function excelFluxReport( mpSim::MPsim, timeGrid::Vector{T}, fluxType::Symbol,
    nodes::AbstractString...; filename::AbstractString = "fluxReport",
    overwrite::Bool = true, timeFactor::Real = 12.0 ) where T <: Real

    filename, overwrite = setupFile( filename, overwrite )
    tStart = now()
    fluxReport = nodeFluxReport( mpSim, timeGrid, fluxType, nodes... )
    reportGenerationTime = ( now() - tStart ).value / 1000.0

    if isempty( fluxReport ) || all( node -> isempty( fluxReport[ node ] ),
        collect( keys( fluxReport ) ) )
        return
    end  # if isempty( fluxReport ) || ...

    writeNodeFluxData( fluxReport, nodes, fluxType, mpSim, timeFactor,
        reportGenerationTime, filename, overwrite )

end  # excelFluxReport( mpSim, timeGrid, fluxType, nodes, filename, overwrite,
     # timeFactor )

excelFluxReport( mpSim::MPsim, timeRes::Real, fluxType::Symbol,
    nodes::AbstractString...; filename::AbstractString = "fluxReport",
    overwrite::Bool = true, timeFactor::Real = 12.0 ) =
    excelFluxReport( mpSim, MP.generateTimeGrid( mpSim, timeRes ), fluxType,
    nodes...; filename = filename, overwrite = overwrite,
    timeFactor = timeFactor )

include( "private/nodeflux.jl" )


text:
module IntoTheJungleCodeModLanternHeartSpawner

using ..Ahorn, Maple

@mapdef Entity "IntoTheJungleCodeMod/LanternHeartSpawner" LanternHeartSpawner(x::Integer, y::Integer)

const placements = Ahorn.PlacementDict(
    "Lantern Heart Spawner (Into The Jungle code mod)" => Ahorn.EntityPlacement(
        LanternHeartSpawner
    )
)

const sprite = "collectables/heartGem/0/00.png"

function Ahorn.selection(entity::LanternHeartSpawner)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::LanternHeartSpawner, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, 0)

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
## Simple non-linear toy model proposed by:
## Netto, M., Gimeno, L. and Mendes, M., 1978. On the optimal and suboptimal
## nonlinear filtering problem for discrete-time systems. IEEE Transactions on
## Automatic Control, 23(6), pp.1062-1067.

module Netto

using SequentialMonteCarlo
using RNGPool
import SMCExamples.Particles.Float64Particle

struct NettoTheta
  σ²::Float64
  δ²::Float64
end

function makeNettoModel(theta::NettoTheta, ys::Vector{Float64})
  n::Int64 = length(ys)
  σ::Float64 = sqrt(theta.σ²)
  invδ²over2::Float64 = 0.5 / theta.δ²
  logncG::Float64 = -0.5 * log(2 * π * theta.δ²)
  @inline function lG(p::Int64, particle::Float64Particle, ::Nothing)
    x::Float64 = particle.x
    tmp::Float64 = x*x/20.0
    v::Float64 = tmp - ys[p]
    return logncG - v * invδ²over2 * v
  end
  @inline function M!(newParticle::Float64Particle, rng::RNG, p::Int64,
    particle::Float64Particle, ::Nothing)
    if p == 1
      newParticle.x = σ*randn(rng)
    else
      x::Float64 = particle.x
      newParticle.x = x/2.0 + 25*x/(1.0+x*x) + 8*cos(1.2*p) + σ*randn(rng)
    end
  end
  model::SMCModel = SMCModel(M!, lG, n, Float64Particle, Nothing)
  return model
end

function simulateNettoModel(theta::NettoTheta, n::Int64)
  model = makeNettoModel(theta, Vector{Float64}(undef, 0))
  ys = Vector{Float64}(undef, n)
  xParticle = Float64Particle()
  rng = getRNG()
  for p in 1:n
    model.M!(xParticle, rng, p, xParticle, nothing)
    ys[p] = xParticle.x*xParticle.x/20.0 + sqrt(theta.σ²)*randn(rng)
  end
  return ys
end

function defaultNettoModel(n::Int64)
  theta = NettoTheta(0.9, 0.6)
  ys = simulateNettoModel(theta, n)

  model = makeNettoModel(theta, ys)

  return model, theta, ys
end

end


text:
using DataStructures

d = Dict('A' => [(1,  1, 'B'), (0, -1, 'C')],
         'B' => [(1, -1, 'A'), (1,  1, 'D')],
         'C' => [(0, -1, 'B'), (0, -1, 'E')],
         'D' => [(1,  1, 'A'), (0,  1, 'B')],
         'E' => [(1, -1, 'F'), (1, -1, 'C')],
         'F' => [(1,  1, 'D'), (1,  1, 'A')])

t = DefaultDict(() -> 0)
s, p = 'A', 0
for _ in 1:12481997
    write, move, s = d[s][t[p]+1]
    t[p] = write
    p += move
end
println(sum(values(t)))


text:
function test_alg_morphism_char_0(A, B, AtoB)
  @test iszero(AtoB(zero(A)))
  @test iszero(AtoB\(zero(B)))
  @test isone(AtoB(one(A)))
  @test isone(AtoB\(one(B)))

  for i = 1:5
    c = rand(A, -10:10)
    d = rand(A, -10:10)
    @test AtoB\(AtoB(c)) == c
    @test AtoB\(AtoB(d)) == d
    @test AtoB(c + d) == AtoB(c) + AtoB(d)
    @test AtoB(c*d) == AtoB(c)*AtoB(d)

    e = rand(B, -10:10)
    f = rand(B, -10:10)
    @test AtoB(AtoB\(e)) == e
    @test AtoB(AtoB\(f)) == f
    @test AtoB\(e + f) == (AtoB\(e)) + (AtoB\(f))
    @test AtoB\(e*f) == (AtoB\(e)) * (AtoB\(f))
  end
end

function test_alg_morphism_char_p(A, B, AtoB)
  @test iszero(AtoB(zero(A)))
  @test iszero(AtoB\(zero(B)))
  @test isone(AtoB(one(A)))
  @test isone(AtoB\(one(B)))

  for i = 1:5
    c = rand(A)
    d = rand(A)
    @test AtoB\(AtoB(c)) == c
    @test AtoB\(AtoB(d)) == d
    @test AtoB(c + d) == AtoB(c) + AtoB(d)
    @test AtoB(c*d) == AtoB(c)*AtoB(d)

    e = rand(B)
    f = rand(B)
    @test AtoB(AtoB\(e)) == e
    @test AtoB(AtoB\(f)) == f
    @test AtoB\(e + f) == AtoB\(e) + AtoB\(f)
    @test AtoB\(e*f) == (AtoB\(e))*(AtoB\(f))
  end
end

@testset "Associative algebras" begin

  @testset "Change of ring" begin

    # Restrict from number field to Q
    Qx, x = FlintQQ["x"]
    f = x^3 - 2
    K, a = number_field(f, "a")

    A = AlgAss(MatrixAlgebra(K, 2))
    B, BtoA = Hecke.restrict_scalars(A, FlintQQ)
    @test base_ring(B) == FlintQQ
    @test dim(B) == dim(A)*degree(K)

    test_alg_morphism_char_0(B, A, BtoA)

    # Extend to K again
    C, CtoB = Hecke._as_algebra_over_center(B)
    @test is_isomorphic(K, base_ring(C))
    @test dim(C) == dim(A)

    test_alg_morphism_char_0(C, B, CtoB)

    # Restrict from number field to number field
    g = x^9 - 15x^6 - 87x^3 - 125
    L, b = number_field(g, "b")
    KtoL = hom(K, L, -2//45*b^7 + 7//9*b^4 + 109//45*b)

    A = AlgAss(MatrixAlgebra(L, 2))
    B, BtoA = Hecke.restrict_scalars(A, KtoL)

    @test base_ring(B) == K
    @test dim(B) == dim(A)*div(degree(L), degree(K))

    test_alg_morphism_char_0(B, A, BtoA)

    # Restrict from F_q to F_p
    Fp = GF(7)
    Fq, a = FiniteField(7, 3, "a")

    A = AlgAss(MatrixAlgebra(Fq, 2))
    B, BtoA = Hecke.restrict_scalars(A, Fp)
    @test base_ring(B) == Fp
    @test dim(B) == dim(A)*degree(K)

    test_alg_morphism_char_p(B, A, BtoA)

    # Extend to Fq again
    C, CtoB = Hecke._as_algebra_over_center(B)
    @test characteristic(base_ring(C)) == characteristic(Fq)
    @test degree(base_ring(C)) == degree(Fq)
    @test dim(C) == dim(A)

    test_alg_morphism_char_p(C, B, CtoB)

    # Extend from F_p^m to F_p^n
    Fqx, x = Fq["x"]
    f = x^2 + 5x + 2
    A = AlgAss(f)
    B, BtoA = Hecke._as_algebra_over_center(A)
    @test characteristic(base_ring(B)) == characteristic(Fq)
    @test degree(base_ring(B)) == degree(f)*degree(Fq)
    @test dim(B) == 1

    test_alg_morphism_char_p(B, A, BtoA)

    Fp = GF(3)
    mt = Array{elem_type(Fp), 3}(undef, 2, 2, 2)
    mt[1, 1, 1] = Fp(0)
    mt[1, 1, 2] = Fp(2)
    mt[1, 2, 1] = Fp(2)
    mt[1, 2, 2] = Fp(1)
    mt[2, 1, 1] = Fp(2)
    mt[2, 1, 2] = Fp(1)
    mt[2, 2, 1] = Fp(1)
    mt[2, 2, 2] = Fp(1)
    A = AlgAss(Fp, mt)
    B, BtoA = Hecke._as_algebra_over_center(A)
    @test characteristic(base_ring(B)) == characteristic(Fp)
    @test degree(base_ring(B)) == dim(A)
    @test dim(B) == 1

    test_alg_morphism_char_p(B, A, BtoA)

  end

  # n = dim(A)^2 = dim(B)^2
  function test_mat_alg_morphism(AtoB::Hecke.AbsAlgAssMor, n::Int)
    A = domain(AtoB)
    B = codomain(AtoB)

    test_alg_morphism_char_p(A, B, AtoB)

    sum_of_diag = AtoB\B[1]
    for i = 2:n
      sum_of_diag += AtoB\B[(i - 1)*n + i]
    end
    @test isone(sum_of_diag)

    # B[i + (j - 1)*n]*B[k + (l - 1)*n] == B[i + (l - 1)*n], if j == k, and 0, otherwise
    for j = 1:n
      jN = (j - 1)*n
      for i = 1:n
        ji = jN + i
        for l = 1:n
          lN = (l - 1)*n
          for k = 1:n
            if j == k
              @test AtoB\(B[ji]*B[k + lN]) == AtoB\B[i + lN]
            else
              @test iszero(AtoB\(B[ji]*B[k + lN]))
            end
          end
        end
      end
    end
  end

  @testset "Matrix Algebra" begin
    Fp = GF(7)
    Fq, a = FiniteField(7, 2, "a")

    A = AlgAss(MatrixAlgebra(Fq, 3))
    B, AtoB = Hecke._as_matrix_algebra(A)
    @test dim(B) == dim(A)

    test_mat_alg_morphism(AtoB, 3)

    G = SymmetricGroup(4)
    A = AlgAss(AlgGrp(Fp, G))[1]
    Adec = Hecke.decompose(A)

    i = 2
    B = Adec[1][1]
    while dim(B) != 9
      B = Adec[i][1]
      i += 1
    end

    C, BtoC = Hecke._as_matrix_algebra(B)

    test_mat_alg_morphism(BtoC, 3)

    i = 2
    B = Adec[1][1]
    while dim(B) != 4
      B = Adec[i][1]
      i += 1
    end

    C, BtoC = Hecke._as_matrix_algebra(B)

    test_mat_alg_morphism(BtoC, 2)
  end

  @testset "Splitting at infinite place" begin
    G = small_group(8, 4)
    Qx, x = PolynomialRing(FlintQQ, "x")
    K, a = NumberField(x - 1, "a")
    A = Hecke.AlgGrp(K, G)
    H = first(c[1] for c in Hecke.decompose(A) if dim(c[1]) == 4)
    P = infinite_places(K)[1]
    @test !is_split(H, P)

    K, a = NumberField(x - 1, "a")
    A = Hecke.AlgGrp(K, G)
    H = first(c[1] for c in Hecke.decompose(A) if dim(c[1]) == 1)
    P = infinite_places(K)[1]
    @test is_split(H, P)

    K, a = NumberField(x^2 - 2, "a")
    HH = Hecke.quaternion_algebra2(2, 3)
    A = AlgAss(K, map(K, HH.mult_table))
    Ps = real_places(K)
    @test is_split(A, Ps[1])
    @test is_split(A, Ps[2])
  end

end


text:
struct EdgeIter{G<:AGraphOrDiGraph}
    g::G
end

#TODO implement iterate protocolo for real
iterate(it::EdgeIter{<:AGraphOrDiGraph}, i=start(it)) = done(it, i) ? nothing : next(it, i)

@inline function start(it::EdgeIter{G}) where G<:AGraph
    s = _start(it)
    while !_done(it, s)
        e, t = _next(it, s)
        if src(e) <= dst(e)
            return (false, e, t)
        end
        s = t
    end
    (true, edge(it.g, 1, 2), s)
end

@inline function next(it::EdgeIter{G}, st) where G<:AGraph
    _, e, s = st
    while !_done(it, s)
        w, t = _next(it, s)
        if src(w) <= dst(w)
            return e, (false, w, t)
        end
        s = t
    end
    e, (true, e, s)
end

@inline done(it::EdgeIter{G}, s) where {G<:AGraph} = s[1]

@inline start(it::EdgeIter{G}) where {G<:ADiGraph} = _start(it)
@inline next(it::EdgeIter{G}, s) where {G<:ADiGraph} = _next(it, s)
@inline done(it::EdgeIter{G}, s) where {G<:ADiGraph} = _done(it, s)

@inline function _start(it::EdgeIter)
    i = 1
    xs_state = 1
    while i <= nv(it.g)
        eit = out_edges(it.g, i)
        xs_state = Base.start(eit)
        !Base.done(eit, xs_state) && break
        i += 1
    end
    return i, xs_state
end

@inline function _next(it::EdgeIter, state)
    i, xs_state = state
    e, xs_state = Base.next(out_edges(it.g, i), xs_state)
    while Base.done(out_edges(it.g, i), xs_state)
        i += 1
        i > nv(it.g) && break
        xs_state = Base.start(out_edges(it.g, i))
    end
    return e, (i, xs_state)
end

@inline _done(it::EdgeIter, state) = state[1] > nv(it.g)

eltype(it::EdgeIter) = edgetype(it.g)
length(it::EdgeIter) = ne(it.g)
IteratorSize(::Type{EdgeIter}) = HasLength()

"""
    edges(g)

Returns an iterator to the edges of a graph `g`.
The returned iterator is invalidated by changes to `g`.
"""
edges(g::AGraphOrDiGraph) = EdgeIter(g)

## Simpler but slower with julia 0.5. Try again with future versions
# edges(g::ADiGraph) = nv(g) == 0 ? #julia issue #18852
#                 (edge(g, u, v) for u=1:1 for v=1:0) :
#                 (e for u=1:nv(g) for e in out_edges(g, u))
#
# edges(g::AGraph) = nv(g) == 0 ? #julia issue #18852
#                 (edge(g, u, v) for u=1:1 for v=1:0) :
#                 (e for u=1:nv(g) for e in out_edges(g, u) if src(e) <= dst(e))
#
# edges(g::ADiGraph, vs::AbstractVector) = length(vs) == 0 ? #julia issue #18852
#                 (edge(g, u, v) for u=1:1 for v=1:0) :
#                 (e for u in vs for e in out_edges(g, u) if dst(e) ∈ vs)
#
# edges(g::AGraph, vs::AbstractVector) = length(vs) == 0 ? #julia issue #18852
#                 (edge(g, u, v) for u=1:1 for v=1:0) :
#                 (e for u in vs for e in out_edges(g, u) if src(e) <= dst(e) && dst(e) ∈ vs)


text:
using JuMP
using ParameterJuMP
using Polyhedra

mutable struct Conitope{T, VT <: AbstractVector{T}, D<:Polyhedra.FullDim}
    d::D
    points::Vector{VT}
    factory::Union{Nothing, JuMP.OptimizerFactory}
    model::Union{Nothing, JuMP.Model}
    z::Union{Nothing, Vector{ParameterJuMP.ParameterRef}}
    t_0::Union{Nothing, JuMP.VariableRef}
    function Conitope{T, VT, D}(
        d::Polyhedra.FullDim, points::Polyhedra.PointIt,
        factory::Union{Nothing, JuMP.OptimizerFactory}=nothing) where {T, VT, D}
        new{T, VT, D}(Polyhedra.FullDim_convert(D, d),
                      Polyhedra.lazy_collect(points), factory, nothing, nothing, nothing)
    end
end
function Conitope(d::Polyhedra.FullDim, points::Polyhedra.PointIt, args...)
    return Conitope{Polyhedra.coefficient_type(eltype(points)), eltype(points),
                    typeof(d)}(d, points, args...)
end

Polyhedra.FullDim(v::Conitope) = v.d

function _build_model(brp::Conitope)
    n = length(brp.points)
    brp.model = ParameterJuMP.ModelWithParams(brp.factory)
    @variable(brp.model, t[1:n] ≥ 0)
    @constraint(brp.model, sum(t) ≤ 1)
    return sum(t[i] .* brp.points[i] for i in 1:n)
end
function _parametrized_model(brp::Conitope, v)
    hull = _build_model(brp)
    brp.z = ParameterJuMP.add_parameters(brp.model, collect(v))
    @constraint(brp.model, brp.z .≤ hull)
end
function _ratio_model(brp::Conitope, v)
    hull = _build_model(brp)
    brp.t_0 = @variable(brp.model)
    @constraint(brp.model, brp.t_0 .* v .≤ hull)
end

function _ei(i::Integer, d::Integer, T::Type)
    e = zeros(T, d)
    e[i] = one(T)
    return e
end
function Polyhedra.polyhedron(
    c::Conitope{T},
    lib::Polyhedra.Library=Polyhedra.default_library(Polyhedra.FullDim(c),
                                                     T)) where T
    n = fulldim(c)
    p = polyhedron(vrep([Ray(-_ei(i, n, T)) for i in 1:n]) + vrep(c.points),
                   lib)
    return p ∩ (Matrix(-LinearAlgebra.I, n, n) * p)
end


text:
## Gridstructure accessors
function get_buses(OPFpbs::OPFProblems, scenario::String)
  return keys(OPFpbs[scenario].ds.bus)
end

function get_links(OPFpbs::OPFProblems, scenario::String)
  return keys(OPFpbs[scenario].ds.link)
end

function get_node_elems(bus::String, OPFpbs::OPFProblems, scenario::String)
  return keys(OPFpbs[scenario].ds.bus[bus])
end

function get_link_elems(link::Link, OPFpbs::OPFProblems, scenario::String)
  return keys(OPFpbs[scenario].ds.link[link])
end

function get_elem_type(elem::String, bus::String, OPFpbs::OPFProblems, scenario::String)
  return typeof(OPFpbs[scenario].ds.bus[bus][elem])
end

function get_elem_type(elem::String, link::Link, OPFpbs::OPFProblems, scenario::String)
  return typeof(OPFpbs[scenario].ds.link[link][elem])
end

function get_links_in(bus::String, gs::GridStructure)
  return haskey(gs.node_linksin, bus) ? gs.node_linksin[bus] : SortedSet{Any}()
end

function get_links_out(bus::String, gs::GridStructure)
  return haskey(gs.node_linksout, bus) ? gs.node_linksout[bus] : SortedSet{Any}()
end

function get_Var(elem_type::String, bus::String, OPFpbs::OPFProblems, scenario::String)
  return OPFpbs[scenario].gs.node_vars[bus][elem_type]
end

function get_Var(elem_type::String, link::Link, OPFpbs::OPFProblems, scenario::String)
  return OPFpbs[scenario].gs.link_vars[link][elem_type]
end

## Modeling Choice accessors
function get_elem_formulation(elem_id::String, bus::String, OPFpbs::OPFProblems, scenario::String)
  return OPFpbs[scenario].mc.node_formulations[bus][elem_id]
end

function get_elem_formulation(elem_id::Type, link::Link, OPFpbs::OPFProblems, scenario::String)
  return OPFpbs[scenario].mc.node_formulations[link][elem_id]
end


text:
mutable struct MaternCovarianceCache
    last_theta::AbstractArray{Float64, 1}
    D::AbstractArray{Float64, 2}
    D2::AbstractArray{Float64, 2}
    K::AbstractArray{Float64, 2}
end

"""
    MaternIsoKernel <: AbstractGPKernel

Matérn kernel with uniform length scale across all dimensions, ``l``.
Parameter ``\\nu`` (nu) is passed in constructor. Currently, only values of
``\\nu=1``, ``\\nu=3`` and ``\\nu=5`` are supported.

```math
\\begin{aligned}
K_{\\nu=1}(r) &= \\sigma_f^2e^{-\\sqrt{r}}\\\\
K_{\\nu=3}(r) &= \\sigma_f^2(1 + \\sqrt{3r})e^{-\\sqrt{3r}}\\\\
K_{\\nu=5}(r) &= \\sigma_f^2(1 + \\sqrt{3r} + \\frac{5}{3}r)e^{-\\sqrt{5r}}\\\\
r_{ij} &= \\sum_{k=1}^d\\frac{(x_{ik}-z_{jk})^2}{l^2}
\\end{aligned}
```

``r_{ij}`` are computed by [`scaled_squared_distance`](@ref)

# Hyperparameters
Hyperparameters vector for this kernel must contain two elements, in the following order:
1. ``\\sigma_f``: the signal standard deviation
2. ``l``: the length scale
"""
struct MaternIsoKernel <: AbstractGPKernel
    nu::Int
    cache::MaternCovarianceCache
end

function check_nu(nu::Int)
    if nu != 1 && nu != 3 && nu != 5
        error("Unsupported ν = $(nu)")
    end
end

function MaternIsoKernel(nu::Int)
    check_nu(nu)
    MaternIsoKernel(nu, MaternCovarianceCache(zeros(0), zeros(0,0), zeros(0,0), zeros(0,0)))
end

"""
    ExponentialIsoKernel

Alias for [`MaternIsoKernel`](@ref)(1)
"""
ExponentialIsoKernel() = MaternIsoKernel(1)

"""
    MaternArdKernel <: AbstractGPKernel

Matérn kernel with distinct length scale for each dimention, ``l_k``.
Parameter ``\\nu`` (nu) is passed in constructor. Currently, only values of
``\\nu=1``, ``\\nu=3`` and ``\\nu=5`` are supported.

```math
\\begin{aligned}
K_{\\nu=1}(r) &= \\sigma_f^2e^{-\\sqrt{r}}\\\\
K_{\\nu=3}(r) &= \\sigma_f^2(1 + \\sqrt{3r})e^{-\\sqrt{3r}}\\\\
K_{\\nu=5}(r) &= \\sigma_f^2(1 + \\sqrt{3r} + \\frac{5}{3}r)e^{-\\sqrt{5r}}\\\\
r_{ij} &= \\sum_{k=1}^d\\frac{(x_{ik}-z_{jk})^2}{l_k^2}
\\end{aligned}
```

``r_{ij}`` are computed by [`scaled_squared_distance`](@ref)

# Hyperparameters
The length of hyperparameters array for this kernel depends on the dimensionality
of the data. Assuming each data point is a vector in a ``d``-dimensional space,
this kernel needs ``d+1`` hyperparameters, in the following order:
1. ``\\sigma_f``: the signal standard deviation
2. ``l_1, \\ldots, l_d``: the length scales for each dimension
"""
struct MaternArdKernel <: AbstractGPKernel
    nu::Int
    cache::MaternCovarianceCache
end

"""
    ExponentialArdKernel

Alias for [`MaternArdKernel`](@ref)(1)
"""
ExponentialArdKernel() = MaternArdKernel(1)

function MaternArdKernel(nu::Int)
    check_nu(nu)
    MaternArdKernel(nu, MaternCovarianceCache(zeros(0), zeros(0,0), zeros(0,0), zeros(0,0)))
end

function get_hyperparameters_size(kernel::MaternIsoKernel, training_data::AbstractArray{Float64, 2})
    return 2
end

function get_hyperparameters_size(kernel::MaternArdKernel, training_data::AbstractArray{Float64, 2})
    return size(training_data, 2) + 1
end

function matern_poly(nu, D, D2)
    if nu == 1
        return 1
    elseif nu == 3
        return D .* sqrt(3) .+ 1
    elseif nu == 5
        return D2 .* (5/3) .+ D * sqrt(5) .+ 1
    end
end

function matern_iso_poly_grad(nu, D, D2)
    if nu == 1
        return D
    elseif nu == 3
        return 3 * D2
    elseif nu == 5
        return D2 .* (5/3) .* (D .* sqrt(5) .+ 1)
    end
end

function matern_ard_poly_grad(nu, D)
    if nu == 1
        ret = -1 ./ D
        ret[findall(D .== 0.0)] .= 0.0
        ret
    elseif nu == 3
        -sqrt(3)
    elseif nu == 5
        -(sqrt(5) .* D .+ 1) .* (sqrt(5) / 3)
    end
end

function matern_covariance_common(nu::Int, cache::MaternCovarianceCache, log_theta::AbstractArray{Float64, 1},
        x1::AbstractArray{Float64, 2}, x2::AbstractArray{Float64, 2})
    sigma_f = exp(log_theta[1] * 2)
    D2 = scaled_squared_distance(log_theta[2:end], x1, x2)
    D2[findall(D2 .< 0)] .= 0.0  # could be negative due to numerical noise
    D = sqrt.(D2)
    K = sigma_f .* matern_poly(nu, D, D2) .* exp.(-sqrt(nu) .* D)
end

function covariance(ker::MaternIsoKernel, log_theta::AbstractArray{Float64, 1},
        x1::AbstractArray{Float64, 2}, x2::AbstractArray{Float64, 2})
    matern_covariance_common(ker.nu, ker.cache, log_theta, x1, x2)
end

function covariance(ker::MaternArdKernel, log_theta::AbstractArray{Float64, 1},
        x1::AbstractArray{Float64, 2}, x2::AbstractArray{Float64, 2})
    matern_covariance_common(ker.nu, ker.cache, log_theta, x1, x2)
end

function matern_covariance_common_training(nu::Int, cache::MaternCovarianceCache,
    log_theta::AbstractArray{Float64, 1}, x::AbstractArray{Float64, 2})
    if log_theta != cache.last_theta ||
            size(cache.D, 1) != size(x, 1) ||
            size(cache.D2, 1) != size(x, 1) ||
            size(cache.K, 1) != size(x, 1)
        sigma_f = exp(log_theta[1] * 2)
        D2 = scaled_squared_distance(log_theta[2:end], x, x)
        D2[findall(D2 .< 0)] .= 0.0  # could be negative due to numerical noise
        D = sqrt.(D2)
        K = sigma_f .* matern_poly(nu, D, D2) .* exp.(-sqrt(nu) .* D)
        cache.last_theta = copy(log_theta)
        cache.D = D
        cache.D2 = D2
        cache.K = K
    end
    cache.K
end

function covariance_training(ker::MaternIsoKernel, log_theta::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2})
    matern_covariance_common_training(ker.nu, ker.cache, log_theta, x)
end

function covariance_training(ker::MaternArdKernel, log_theta::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2})
    matern_covariance_common_training(ker.nu, ker.cache, log_theta, x)
end

function covariance_diagonal(ker::MaternIsoKernel, log_theta::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2})
    fill(exp(log_theta[1] * 2), (size(x, 1), 1))
end

function covariance_diagonal(ker::MaternArdKernel, log_theta::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2})
    fill(exp(log_theta[1] * 2), (size(x, 1), 1))
end

function matern_covariance_grad_common(nu::Int, cache::MaternCovarianceCache, log_theta::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2}, R::AbstractArray{Float64, 2})
    sigma_f = exp(log_theta[1] * 2)
    if log_theta == cache.last_theta &&
            size(cache.D, 1) == size(x, 1) &&
            size(cache.D2, 1) != size(x, 1) &&
            size(cache.K, 1) != size(x, 1)
        K = cache.K
        D = cache.D
        D2 = cache.D2
    else
        D2 = scaled_squared_distance(log_theta[2:end], x, x)
        D2[findall(D2 .< 0)] .= 0.0  # could be negative due to numerical noise
        D = sqrt.(D2)
        K = sigma_f .* matern_poly(nu, D, D2) .* exp.(-sqrt(nu) .* D)
        cache.last_theta = copy(log_theta)
        cache.D = D
        cache.D2 = D2
        cache.K = K
    end
    KR = K .* R
    d_sigma_f = 2 * sum(KR)
    exp_nu_D_R = exp.(-sqrt(nu) .* D) .* R
    if length(log_theta) == 2
        d_ell = sum(matern_iso_poly_grad(nu, D, D2) .* exp_nu_D_R)
    else
        d_ell = scaled_squared_distance_grad(log_theta[2:end], x, x,
            sqrt(nu) / 2 .* exp_nu_D_R .* matern_ard_poly_grad(nu, D))
    end
    return [d_sigma_f; sigma_f * d_ell]
end

function covariance_grad(ker::MaternIsoKernel, log_theta::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2}, R::AbstractArray{Float64, 2})
    matern_covariance_grad_common(ker.nu, ker.cache, log_theta, x, R)
end

function covariance_grad(ker::MaternArdKernel, log_theta::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2}, R::AbstractArray{Float64, 2})
    matern_covariance_grad_common(ker.nu, ker.cache, log_theta, x, R)
end