text:
#!/usr/bin/env julia --project
include("utilities/boilerplate.jl")
using ClimaCoupler

########
# Set up parameters and initial conditions
########
include("parameters_initialconditions.jl")

FT = Float64

# Collects tracer values/fluxes for conservation checks
mutable struct LogThatFlux{F}
    A::F
    B::F
end

########
# Set up domain
########
# A: low atmos level (troposphere)
domainA = SphericalShell(radius = FT(1), height = FT(0.2))
# B: high atmos level (stratosphere)
domainB = SphericalShell(radius = domainA.radius + domainA.height, height = domainA.height)

gridA = DiscretizedDomain(
    domainA;
    elements = (vertical = 1, horizontal = 8),
    polynomial_order = (vertical = 0, horizontal = 3),
    overintegration_order = (vertical = 1, horizontal = 1),
)

gridB = DiscretizedDomain(
    domainB;
    elements = (vertical = 1, horizontal = 8),
    polynomial_order = (vertical = 0, horizontal = 3),
    overintegration_order = (vertical = 1, horizontal = 1),
)

########
# Set up model physics
########
physics = Physics(
    orientation = SphericalOrientation(),
    advection = NonLinearAdvection(),
    coriolis = DeepShellCoriolis{Float64}(Ω = parameters.Ω),
    gravity = Buoyancy{Float64}(α = parameters.α, g = 0.0),
    eos = BarotropicFluid{Float64}(ρₒ = parameters.ρₒ, cₛ = parameters.cₛ),
)

########
# Set up boundary conditions
########
bcsA = (
    bottom = (ρu = Impenetrable(FreeSlip()), ρθ = Insulating()),
    # top = (ρu = Impenetrable(FreeSlip()), ρθ = Insulating()), # for non-coupled testing
    top = (ρu = Impenetrable(FreeSlip()), ρθ = CoupledSecondaryBoundary()),
)
bcsB = (
    bottom = (ρu = Impenetrable(FreeSlip()), ρθ = CoupledPrimaryBoundary()),
    # bottom = (ρu = Impenetrable(FreeSlip()), ρθ = Insulating()), #for non-coupled testing
    top = (ρu = Impenetrable(FreeSlip()), ρθ = Insulating()),
)

########
# Set up model
########
modelA = ModelSetup(
    physics = physics,
    boundary_conditions = bcsA,
    initial_conditions = (ρ = ρ₀ᶜᵃʳᵗ, ρu = ρu⃗₀ᶜᵃʳᵗ, ρθ = ρθ₀ᶜᵃʳᵗA),
    numerics = (flux = RoeNumericalFlux(), flux_second_order = PenaltyNumFluxDiffusive()),
    parameters = parameters,
)

modelB = ModelSetup(
    physics = physics,
    boundary_conditions = bcsB,
    initial_conditions = (ρ = ρ₀ᶜᵃʳᵗ, ρu = ρu⃗₀ᶜᵃʳᵗ, ρθ = ρθ₀ᶜᵃʳᵗB),
    numerics = (flux = RoeNumericalFlux(), flux_second_order = PenaltyNumFluxDiffusive()),
    parameters = parameters,
)

########
# Set up time steppers (could be done automatically in simulation)
########
Δt = min_node_distance(gridA.numerical) / parameters.cₛ * 0.25
total_steps = 10000
start_time = 0
end_time = Δt * total_steps#30 * 24 * 3600
method = SSPRK22Heuns
callbacksA = (
    Info(),
    CFL(),
    #   VTKState(
    #     iteration = 250,#Int(floor(6*3600/Δt)),
    #     filepath = "./out/A/"),
)
callbacksB = (
    Info(),
    CFL(),
    #   VTKState(
    #     iteration = 250,#Int(floor(6*3600/Δt)),
    #     filepath = "./out/B/"),
)

########
# Set up simulation
########
nstepsA = 1
nstepsB = 1

epss = sqrt(eps(FT))
boundary_mask(param_set, xc, yc, zc) = @. abs((xc^2 + yc^2 + zc^2)^0.5 - domainB.radius) < epss

simA = CplSimulation(
    modelA;
    grid = gridA,
    timestepper = (method = method, timestep = Δt / nstepsA),
    time = (start = start_time, finish = end_time),
    nsteps = nstepsA,
    boundary_z = boundary_mask,
    callbacks = callbacksA,
)
simB = CplSimulation(
    modelB;
    grid = gridB,
    timestepper = (method = method, timestep = Δt / nstepsB),
    time = (start = start_time, finish = end_time),
    nsteps = nstepsB,
    boundary_z = boundary_mask,
    callbacks = callbacksB,
)

## Create a Coupler State object for holding import/export fields.
coupler = CplState()
coupler_register!(coupler, :EnergyA, deepcopy(simA.state.ρθ[simA.boundary]), simA.grid, DateTime(0), u"J") # value on top of domainA for calculating upward flux into domainB
coupler_register!(coupler, :EnergyFluxB, deepcopy(simB.state.F_ρθ_accum[simB.boundary]), simB.grid, DateTime(0), u"J") # downward flux

compA = (pre_step = preA, component_model = simA, post_step = postA)
compB = (pre_step = preB, component_model = simB, post_step = postB)
component_list = (domainA = compA, domainB = compB)
cpl_solver = CplSolver(
    component_list = component_list,
    coupler = coupler,
    coupling_dt = Δt,
    t0 = FT(start_time),
    fluxlog = LogThatFlux(zeros(total_steps), zeros(total_steps)),
)

########
# Run the simulation
########
evolve!(cpl_solver, total_steps)

########
# Check conservation
########
using Plots
fluxA = cpl_solver.fluxlog.A
fluxB = cpl_solver.fluxlog.B
fluxT = fluxA .+ fluxB
time = collect(1:1:total_steps)
rel_error = [((fluxT .- fluxT[1]) / fluxT[1])]
plot(time .* cpl_solver.dt, rel_error)
#plot(time .* cpl_solver.dt, [fluxA .- fluxA[1] fluxB .- fluxB[1]])

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
using MultivectorAlgebras
using Base.Test: @test, @test_throws

# Multivector1: same typed arguments
@test begin
    l = Multivector1(1.2, 3.4)
    r = Multivector1{Float64}(1.2, 3.4)
    l == r
end

# Multivector1: mixed typed arguments
@test begin
    l = Multivector1(1, 2.3)
    r = Multivector1{Float64}(1.0, 2.3)
    l == r
end

# Multivector1: single argument
@test begin
    l = Multivector1(2.3)
    r = Multivector1{Float64}(2.3, 0.0)
    l == r
end

# real
@test begin
    l = real(Multivector1(1, 2))
    r = 1
    l == r
end

# unreal: Multivector1
@test begin
    l = unreal(Multivector1(3, 4))
    r = 4
    l == r
end

# unreal: Real
@test begin
    l = unreal(3.14)
    r = 0.0
    l == r
end

# zero: zero-like
@test begin
    l = zero(Multivector1(1, 2))
    r = Multivector1(0, 0)
    l == r
end

# zero: type argument
@test begin
    l = zero(Multivector1{Float64})
    r = Multivector1(0.0, 0.0)
    l == r
end

# one: one-like
@test begin
    l = one(Multivector1(1, 2))
    r = Multivector1(1, 0)
    l == r
end

# one: type argument
@test begin
    l = one(Multivector1{Float64})
    r = Multivector1(1.0, 0.0)
    l == r
end

# conj
@test begin
    l = conj(Multivector1(5, 6))
    r = Multivector1(5, -6)
    l == r
end

# Distributivity of conj
@test begin
    a = Multivector1(1, 2)
    b = Multivector1(3, 4)
    l = conj(a * b)
    r = conj(a) * conj(b)
    l == r
end

# cloak
@test begin
    l = cloak(Multivector1(5, 6))
    r = Multivector1(-5, 6)
    l == r
end

# Anti-distributivity of cloak
@test begin
    a = Multivector1(1, 2)
    b = Multivector1(3, 4)
    l = -cloak(a * b)
    r = cloak(a) * cloak(b)
    l == r
end

# cloak: equivalent to -dagger
@test begin
    a = Multivector1(7, 8)
    l = cloak(a)
    r = -dagger(a)
    l == r
end

# dagger
@test begin
    l = dagger(Multivector1(1, 2))
    r = Multivector1(1, -2)
    l == r
end

# Distributivity of dagger
@test begin
    a = Multivector1(1, 2)
    b = Multivector1(3, 4)
    l = dagger(a * b)
    r = dagger(a) * dagger(b)
    l == r
end

# dagger: equivalent to conj
@test begin
    a = Multivector1(9, 8)
    l = dagger(a)
    r = conj(a)
    l == r
end

# star
@test begin
    l = star(Multivector1(2, 3))
    r = Multivector1(3, 2)
    l == r
end

# Involutivity of conj
@test begin
    l = conj(conj(Multivector1(1, 2)))
    r = Multivector1(1, 2)
    l == r
end

# Involutivity of cloak
@test begin
    l = cloak(cloak(Multivector1(1, 2)))
    r = Multivector1(1, 2)
    l == r
end

# Involutivity of dagger
@test begin
    l = dagger(dagger(Multivector1(1, 2)))
    r = Multivector1(1, 2)
    l == r
end

# Involutivity of star
@test begin
    l = star(star(Multivector1(1, 2)))
    r = Multivector1(1, 2)
    l == r
end

# +: unary
@test begin
    l = +(Multivector1(1, 2))
    r = Multivector1(1, 2)
    l == r
end

# +
@test begin
    l = Multivector1(1, 2) + Multivector1(3, 4)
    r = Multivector1(4, 6)
    l == r
end

# +: addition of real
@test begin
    l = 2 + Multivector1(8, 9)
    r = Multivector1(10, 9)
    l == r
end

# +: addition of real
@test begin
    l = Multivector1(8, 9) + 1
    r = Multivector1(9, 9)
    l == r
end

# zero is addition identity
@test begin
    l = Multivector1(5, 6) + zero(Multivector1(5, 6))
    r = Multivector1(5, 6)
    l == r
end

# Commutativity of addition
@test begin
    l = Multivector1(1, 2) + Multivector1(3, 4)
    r = Multivector1(3, 4) + Multivector1(1, 2)
    l == r
end

# Associativity of addition
@test begin
    l = (Multivector1(1, 2) + Multivector1(3, 4)) + Multivector1(5, 6)
    r = Multivector1(1, 2) + (Multivector1(3, 4) + Multivector1(5, 6))
    l == r
end

# -
@test begin
    l = Multivector1(1, 2) - Multivector1(3, 4)
    r = Multivector1(-2, -2)
    l == r
end

# -: subtraction from real
@test begin
    l = 2 - Multivector1(8, 9)
    r = Multivector1(-6, -9)
    l == r
end

# -: subtraction of real
@test begin
    l = Multivector1(8, 9) - 1
    r = Multivector1(7, 9)
    l == r
end

# -: minus
@test begin
    l = -Multivector1(7, 8)
    r = Multivector1(-7, -8)
    l == r
end

# Anti-commutativity of subtraction
@test begin
    l = Multivector1(1, 2) - Multivector1(3, 4)
    r = -(Multivector1(3, 4) - Multivector1(1, 2))
    l == r
end

# Associativity of subtraction
@test begin
    l = (Multivector1(1, 2) - Multivector1(3, 4)) - Multivector1(5, 6)
    r = Multivector1(1, 2) - (Multivector1(3, 4) + Multivector1(5, 6))
    l == r
end

# *
@test begin
    l = Multivector1(2, 3) * Multivector1(4, 5)
    r = Multivector1(8, 22)
    l == r
end

# *: multiplication by real
@test begin
    l = 5 * Multivector1(2, 3)
    r = Multivector1(10, 15)
    l == r
end

# *: multiplication by real
@test begin
    l = Multivector1(3, 4) * 3
    r = Multivector1(9, 12)
    l == r
end

# one is multiplication identity
@test begin
    l = Multivector1(2, 3) * one(Multivector1(2, 3))
    r = Multivector1(2, 3)
    l == r
end

# Commutativity of multiplication
@test begin
    l = Multivector1(1, 2) * Multivector1(3, 4)
    r = Multivector1(3, 4) * Multivector1(1, 2)
    l == r
end

# Associativity of multiplication
@test begin
    l = (Multivector1(1, 2) * Multivector1(3, 4)) * Multivector1(5, 6)
    r = Multivector1(1, 2) * (Multivector1(3, 4) * Multivector1(5, 6))
    l == r
end

# Nilpotency
@test begin
    l = Multivector1(0, 8) * Multivector1(0, 8)
    r = Multivector1(0, 0)
    l == r
end

# Positivity
@test begin
    l = abs2(Multivector1(6, 7)) > 0
    r = true
    l == r
end

# Composition
@test begin
    a = Multivector1(1, 2)
    b = Multivector1(3, 4)
    l = abs2(a) * abs2(b)
    r = abs2(a * b)
    l == r
end

# \
@test begin
    a = Multivector1(1.0, 2.0)
    b = Multivector1(3.0, 4.0)
    c = a * b
    l = a \ c
    r = b
    l == r
end

# \: zero divisor error
@test_throws ErrorException begin
    Multivector1(0, 1) \ Multivector1(2, 3)
end

# \: zero divisor error
@test_throws ErrorException begin
    Multivector1(0, 1) \ 2
end

# /
@test begin
    a = Multivector1(1.0, 2.0)
    b = Multivector1(3.0, 4.0)
    c = a * b
    l = c / b
    r = a
    l == r
end

# /: zero divisor error
@test_throws ErrorException begin
    Multivector1(1, 2) / Multivector1(0, 3)
end

# /: zero divisor error
@test_throws ErrorException begin
    2 / Multivector1(0, 3)
end

# /: division by real
@test begin
    l = Multivector1(1.2, 3.4) / 5
    r = Multivector1(1.2 / 5, 3.4 / 5)
    l == r
end

# /: division of real
@test begin
    l = 2 / Multivector1(3 // 4, 5 // 6)
    r = 2 * inv(Multivector1(3 // 4, 5 // 6))
    l == r
end

# /: zero denominator
@test_throws ErrorException begin
    Multivector1(1, 2) / 0
end

# \: division by real
@test begin
    l = 5 \ Multivector1(1.2, 3.4)
    r = Multivector1(5 \ 1.2, 5 \ 3.4)
    l == r
end

# \: division of real
@test begin
    l = Multivector1(3 // 4, 5 // 6) \ 2
    r = inv(Multivector1(3 // 4, 5 // 6)) * 2
    l == r
end

# \: zero denominator
@test_throws ErrorException begin
    0 \ Multivector1(1, 2)
end

# inv
@test begin
    a = Multivector1(1 // 2, 3 // 4)
    b = Multivector1(5 // 6, 7 // 8)
    l = a / b
    r = inv(b / a)
    l == r
end

# inv: zero divisor error
@test_throws ErrorException begin
    inv(Multivector1(0, 1))
end

# isreal: false case
@test begin
    l = isreal(Multivector1(1, 2))
    r = false
    l == r
end

# isreal: true case
@test begin
    l = isreal(Multivector1(3, 0))
    r = true
    l == r
end

# asarray
@test begin
    l = asarray(Multivector1(4, 5))
    r = [4, 5]
    l == r
end

# iszero: false case
@test begin
    l = iszero(Multivector1(1, 2))
    r = false
    l == r
end

# iszero: true case
@test begin
    l = iszero(Multivector1(0, 0))
    r = true
    l == r
end

# abs: negative case
@test begin
    l = abs(Multivector1(-2, 3))
    r = 2
    l == r
end

# abs: positive case
@test begin
    l = abs(Multivector1(2, 3))
    r = 2
    l == r
end

# abs2
@test begin
    l = abs2(Multivector1(4, 5))
    r = 16
    l == r
end

# iszerodivisor: false case
@test begin
    l = iszerodivisor(Multivector1(2, 3))
    r = false
    l == r
end

# iszerodivisor: true case
@test begin
    l = iszerodivisor(Multivector1(0, 3))
    r = true
    l == r
end

# selfstar
@test begin
    a = 1.2
    b = 3.4
    l = selfstar(Multivector1(a, b))
    r = Multivector1((a + b) / 2, (a + b) / 2)
    l == r
end

# Involutivity of selfstar
@test begin
    z = Multivector1(1.2, 3.4)
    l = selfstar(selfstar(z))
    r = selfstar(z)
    l == r
end

# antiselfstar
@test begin
    a = 1.2
    b = 3.4
    l = antiselfstar(Multivector1(a, b))
    r = Multivector1((a - b) / 2, (b - a) / 2)
    l == r
end

# Involutivity of antiselfstar
@test begin
    z = Multivector1(1.2, 3.4)
    l = antiselfstar(antiselfstar(z))
    r = antiselfstar(z)
    l == r
end

# show
@test begin
    io = IOBuffer()
    show(io, Multivector1(1, 2))
    l = String(take!(io))
    r = "[1: 1, A: 2]"
    l == r
end

# random
@test begin
    l = random(Multivector1{Int})
    r = random(Multivector1{Int})
    l != r
end

ext: jl


text:
## Setting Struct

"Store the settings required for performing the EMD."
struct EMDSetting
    "Number of siftings"
    num_siftings::Int64
    "S number"
    s_num::Int64
    "Number of IMFS"
    m::Int64
    function EMDSetting(n::Int64, num_siftings::Int64, s_num::Int64, m::Int64)
        if (s_num == 0 && num_siftings == 0) || (num_siftings < 0 || s_num < 0)
            throw(DomainError("Invalid EMD Settings."))
        end

        if m <= 0
            @warn "Invalid number of IMFs, will set to a default number."
            m = num_imfs(n)
        end

        return new(num_siftings, s_num, m)
    end
end

## Main Function

"""
    emd(input::Vector{Float64}, s::EMDSetting)

Compute EMD of the input signal with given settings.
"""
function emd(input_::Vector{Float64}, s::EMDSetting)
    input = deepcopy(input_)
    residue = deepcopy(input_)
    n = length(residue)

    output = zeros(n, s.m)

    for j = 1:s.m-1
        if j != 1
            input = deepcopy(residue)
        end

        sift!(input, s)

        residue -= input
        output[:,j] += input
    end

    output[:,s.m] += residue

    return output
end

## Helper functions

"""
    sift!(input::Vector{Float64}, s::EMDSetting)

Return the IMF that satisfies the given settings. In particular, it returns
an IMF that satisifes the S number criterion that is found within the set
number of siftings.
"""
function sift!(input::Vector{Float64}, s::EMDSetting)
    n = length(input)
    max_spline = zeros(n)
    min_spline = zeros(n)
    max_x = zeros(n)
    max_y = zeros(n)
    min_x = zeros(n)
    min_y = zeros(n)
    sift_counter = 0
    s_counter = 0

    num_max = -1
    num_min = -1
    num_zc = -1
    prev_num_max = -1
    prev_num_min = -1
    prev_num_zc = -1

    while (s.num_siftings == 0 || sift_counter < s.num_siftings)
        sift_counter += 1
        prev_num_max = num_max
        prev_num_min = num_min
        prev_num_zc = num_zc
        num_max, num_min, num_zc = find_extrema!(input, max_x, max_y,
                                                    min_x, min_y)

        if (s.s_num != 0)
            max_diff = num_max - prev_num_max
            min_diff = num_min - prev_num_min
            zc_diff = num_zc - prev_num_zc
            if abs(max_diff) + abs(min_diff) + abs(zc_diff) <= 1
                s_counter += 1
                if s_counter >= s.s_num
                    num_diff = num_min + num_max - 4 - num_zc
                    if abs(num_diff) <= 1
                        break;
                    end
                end
            end
        else
            s_counter = 0
        end
    end

    max_spline = evaluate_spline(max_x, max_y, num_max)
    min_spline = evaluate_spline(min_x, min_y, num_min)

    for i in eachindex(input)
        global input[i] -= 0.5(max_spline[i] + min_spline[i])
    end
end

"""
    find_extrema!(x::Vector{Float64}, max_x::Vector{Float64}, max_y::Vector{Float64},
                  min_x::Vector{Float64}, min_y::Vector{Float64})

Return the number of maxima, minima, and zero crossings of x after modifying max_x,
max_y, min_x, min_y to contain the maxima and minima of x.
"""
function find_extrema!(x::Vector{Float64},
                          max_x::Vector{Float64}, max_y::Vector{Float64},
                          min_x::Vector{Float64}, min_y::Vector{Float64})
    N = length(x)
    num_max = 1
    num_min = 1
    num_zc = 1

    if N == 0
        return num_max-1, num_min-1, num_zc-1
    end

    global max_x[1] = 1
    global max_y[1] = x[1]
    num_max += 1
    global min_x[1] = 1
    global min_y[1] = x[1]
    num_min += 1

    if N == 1
        return num_max-1, num_min-1, num_zc-1
    end

    prev_slope = 0
    prev_sign = (x[1] < 0) ? -1 : ((x[1] > 0) ? 1 : 0)
    flat_counter = 0

    for i=1:N-1
        if x[i+1] > x[i]
            if prev_slope == -1
                global min_x[num_min] = convert(Float64, i) - convert(Float64, flat_counter/2)
                global min_y[num_min] = x[i]
                num_min += 1
            end
            if prev_sign == -1 && x[i+1] > 0
                num_zc += 1
                prev_sign = 1
            elseif prev_sign == 0 && x[i+1] > 0
                prev_sign = 1
            end
            prev_slope = 1
            flat_counter = 0
        elseif x[i+1] < x[i]
            if prev_slope == 1
                global max_x[num_max] = convert(Float64, i) - convert(Float64, flat_counter/2)
                global max_y[num_max] = x[i]
                num_max +=1
            end
            if prev_sign == 1 && x[i+1] < 0
                num_zc += 1
                prev_sign = -1
            elseif prev_sign == 0 && x[i+1] < 0
                prev_sign = -1
            end
            prev_slope = -1
            flat_counter = 0
        else
            flat_counter += 1
        end
    end

    global max_x[num_max] = N
    global max_y[num_max] = x[N]
    num_max += 1

    global min_x[num_min] = N
    global min_y[num_min] = x[N]
    num_min += 1

    if num_max >= 4
        max_el = linear_extrapolate(max_x[2], max_y[2], max_x[3], max_y[3], 0)
        if max_el > max_y[1]
            global max_y[1] = max_el
        end
        max_er = linear_extrapolate(max_x[num_max-3], max_y[num_max-3], max_x[num_max-2], max_y[num_max-2], N)
        if max_er > max_y[num_max]
            global max_y[num_max-1] = max_er
        end
    end
    if num_min >= 4
        min_el = linear_extrapolate(min_x[2], min_y[2], min_x[3], min_y[3], 0)
        if min_el < min_y[1]
            global min_y[1] = min_el
        end
        min_er = linear_extrapolate(min_x[num_min-3], min_y[num_min-3], min_x[num_min-2], min_y[num_min-2], N)
        if min_er < min_y[num_min]
            global min_y[num_min-1] = min_er
        end
    end

    return num_max-1, num_min-1, num_zc-1
end

"""
    linear_extrapolate(x_0::Float64, y_0::Float64, x_1::Float64, y_1::Float64, x::Int64)

Return the linear extrapolation of x based on x_0, x_1, y_0, y_1.
"""
function linear_extrapolate(x_0::Float64, y_0::Float64, x_1::Float64, y_1::Float64, x::Int64)
    if x_0 == x_1
        throw(DomainError(x_1, "x_1 can not equal x_0"))
    end

    return y_0 + (y_1 - y_0)*(x - x_0)/(x_1 - x_0)
end

"""
    evaluate_spline(x::Vector{Float64}, y::Vector{Float64}, n::Int64)

Return spline generated by the first n elements of (x,y).
"""
function evaluate_spline(x::Vector{Float64}, y::Vector{Float64}, n::Int64)
    max_j = x[n]
    N = length(x)

    if n <= 1
        throw(DomainError(n, "Not enough points for a spline."))
    end

    spline_y = zeros(N)

	if n == 2
        spl = Spline1D(x[1:n], y[1:n]; k=1)
        for j=1:Int(max_j)
			spline_y[j] = spl(j)
		end

		return spline_y
	elseif n == 3
        spl = Spline1D(x[1:n], y[1:n]; k=2)
        for j=1:Int(max_j)
			spline_y[j] = spl(j)
		end

		return spline_y
	end

    sys_size::Int32 = n-2
	diag = zeros(sys_size)
	sup_diag = zeros(sys_size-1)
	sub_diag = zeros(sys_size-1)
	g = zeros(sys_size)

	h_0 = x[2] - x[1]
	h_1 = x[3] - x[2]
	h_nm1 = x[n] - x[n-1]
	h_nm2 = x[n-1] - x[n-2]

	diag[1] = h_0 + 2*h_1
	sup_diag[1] = h_1 - h_0
	g[1] = 3/(h_0 + h_1)*((y[3] - y[2]) - (h_1/h_0)*(y[2] - y[1]))

	for i=3:n-2
		h_i = x[i+1] - x[i]
		h_im1 = x[i] - x[i-1]

		sub_diag[i-2] = h_im1
		diag[i-1] = 2*(h_im1 + h_i)
		sup_diag[i-1] = h_i
		g[i-1] = 3*((y[i+1] - y[i])/h_i - (y[i] - y[i-1])/h_im1)
	end

	sub_diag[sys_size-1] = h_nm2 - h_nm1
	diag[sys_size] = 2*h_nm2 + h_nm1
	g[sys_size] = 3/(h_nm1 + h_nm2)*((h_nm2/h_nm1)*(y[n]-y[n-1]) - (y[n-1] - y[n-2]))

	LAPACK.gtsv!(sub_diag, diag, sup_diag, g)

	c = zeros(n)

	for i=2:n-1
		c[i] = g[i-1]
	end

	c[1] = c[2] + (h_0/h_1)*(c[2]-c[3])
	c[n] = c[n-1] + (h_nm1/h_nm2)*(c[n-1]-c[n-2])

	i = 1
    for j=1:Int(max_j)
		if (j > x[i+1])
			i += 1
            if i >= n
                throw(DomainError(i, "exceeded spline boundary."))
            end
		end

		dx = j - x[i]

		if dx == 0
            spline_y[j] = y[i]
			continue
		end

		h_i = x[i+1] - x[i]
        a_i = y[i]
		b_i = (y[i+1] - y[i])/h_i - (h_i/3)*(c[i+1] + 2*c[i])
		c_i = c[i]
		d_i = (c[i+1] - c[i])/(3*h_i)
        spline_y[j] = a_i + dx*(b_i + dx*(c_i + dx*d_i))
	end

	return spline_y
end


text:
using OptEvapCool
using StatsPlots
using Printf: @sprintf
using LaTeXStrings
using Plots.PlotMeasures

function harmonic_verification_trap(duration = 2, input_dir = "")
    # Physical parameters
    Np = 1e7 # Initial atom count
    T₀ = 1e-6 # Initial temperature
    species = Rb87

    # Numerical parameters
    Nt = ceil(Int64, 1e5) # Test particles
    F = Np / Nt
    Nc = 3 # Target number of test particles per cell

    # Trapping frequencies
    m = species.m
    ωx = 2π * 150
    ωy = 2π * 150
    ωz = 2π * 15

    # Cloud initialisation
    positions = harmonic_boltzmann_positions(Nt, m, T₀, ωx, ωy, ωz)
    velocities = boltzmann_velocities(Nt, m, T₀)

    # Three-body loss & background loss
    K = 4.3e-29 * 1e-12
    τ_bg = 180

    # Maximum timestep
    max_dt = 0.05 * 2π / max(ωx, ωy, ωz)

    # HARMONIC SIMULATION
    field = HarmonicField(ωx, ωy, ωz)
    acc! = acceleration(field)
    pot! = potential(field)

    function accel(p, s, t, o)
        fill!(o, 0.0)
        acc!(p, s, t, o)
        return nothing
    end

    function poten(p, s, t, o)
        fill!(o, 0.0)
        pot!(p, s, t, o)
        return nothing
    end

    # Evaporation
    η = 10
    Uₜ = η * kB * T₀
    evap = energy_evap(Uₜ, poten)

    sensor = GlobalSensor()
    measure = measurer(sensor, 0.001, ωx, ωy, ωz)

    conditions = SimulationConditions(species, F, positions, velocities,
        accel, poten, evap = evap, τ_bg = τ_bg, K = K)

    if input_dir == ""
        final_cloud = evolve(conditions, duration;
            Nc = Nc, max_dt = max_dt, measure = measure)
    else
        sensor = loadsensor("$input_dir/sensor-data.csv")
    end

    # Save sensor data (before plotting, in case plotting fails)
    ft = filetime()
    dir = "./results/$ft-harmonic-verification"
    mkpath(dir)

    savecsv(sensor, "$dir/sensor-data.csv")

    # THEORY
    # Pethick model
    peth_t, peth_N, peth_T, peth_Γ = harmonic_theory(
        species, ωx, ωy, ωz, T₀, Np, Uₜ, duration, τ_bg
    )
    # Purdue model
    pur_t, pur_N, pur_T, pur_Γ = purdue_theory(
        species, ωx, ωy, ωz, T₀, Np, Uₜ, duration, τ_bg, K
    )

    # PLOTTING
    default(fontfamily="Computer Modern",
        linewidth=2, framestyle=:box, label=nothing, grid=false)
    scalefontsizes()
    scalefontsizes(1.5)

    window_time = min(4e-2, duration/2)
    time = sensor.time
    window_size = OptEvapCool.window_time_size(time, window_time)
    rolling_time = rollmean(time, window_size)

    # Colour palette
    col1 = RGB(0.0588, 0, 0.5882)
    col2 = RGB(0.9000, 0, 0.3490)
    col3 = RGB(1, 0.6510, 0)
    col4 = RGB(0.0824, 0.8392, 0)

    # Temperature
    rolling_temp = 2 / (3 * kB) * rollmean(sensor.ke, window_size)

    temp_order = -6
    temp_yformatter(y) = @sprintf("%.1f",y/(10.0^temp_order))

    temperature_plt = plot(rolling_time, rolling_temp,
        xlabel = "Time (s)",
        ylabel = L"\textrm{Temperature\ \ }(\mu K)",
        label = false,
        linecolor = col1,
        ls = :solid,
        #ylims = (0, Inf),
        minorticks = 5,
        yformatter = temp_yformatter,
        dpi = 300)
    #plot!(peth_t, peth_T, label = false, linecolor = col2, ls = :dash)
    plot!(pur_t, pur_T, label = false, linecolor = col3, ls = :dash)

    # Number
    Np_series = sensor.Nt .* sensor.F

    num_order = 6
    num_yformatter(y) = @sprintf("%.1f",y/(10.0^num_order))

    number_plt = plot(time, Np_series,
        xlabel = "Time (s)",
        ylabel = L"\textrm{Number\ \ }({}\times10^{%$num_order})",
        label = false,
        linecolor = col1,
        ls = :solid,
        #ylims = (0, Inf),
        minorticks = 5,
        yformatter = num_yformatter,
        dpi = 300)
    #plot!(peth_t, peth_N, label = false, linecolor = col2, ls = :dash)
    plot!(pur_t, pur_N, label = false, linecolor = col3, ls = :dash)

    # Collrate
    rolling_timesteps = (
        time[window_size:end] - time[1:end - window_size + 1]
    )
    rolling_Nt = rollmean(sensor.Nt, window_size)
    rolling_Γ = (2 * rolling(sum, sensor.coll, window_size) ./
        (rolling_Nt .* rolling_timesteps)
    )

    collrate_plt = plot(rolling_time, rolling_Γ,
        xlabel = "Time (s)",
        ylabel = "Collision rate (Hz)",
        label = false,
        linecolor = col1,
        ylims = (250, 350),
        ls = :solid,
        dpi = 300)
    #plot!(peth_t, peth_Γ, label = false, linecolor = col2, ls = :dash)
    plot!(pur_t, pur_Γ, label = false, linecolor = col3, ls = :dash)

    # Phase space density
    rolling_n0 = rollmean(sensor.n0 .* sensor.F, window_size)

    rolling_psd = rolling_n0 .* (
        (2*π*h̄^2 ./ (m * kB .* rolling_temp)) .^1.5
    )

    psd_plt = plot(rolling_time, rolling_psd,
        xlabel = "Time (s)",
        ylabel = "Phase space density",
        label = false,
        linecolor = col1,
        ls = :solid,
        ylims = (0, Inf),
        minorticks = 5,
        dpi = 300)

    # Save plots
    savefig(temperature_plt, "$dir/temp.png")
    savefig(collrate_plt, "$dir/collrate.png")
    savefig(number_plt, "$dir/number.png")
    savefig(psd_plt, "$dir/psd.png")

    return nothing
end


text:
"""
```
find_threshold(Entropy(),counts,edges)
```
An algorithm for finding the threshold value for a gray-level histogram using
the entropy of the histogram.

# Output

Returns the point in the `AbstractRange` which corresponds to the threshold bin in the histogram.

# Details

This algorithm uses the entropy of a gray level histogram to produce a threshold
value.

Let ``f_1, f_2, \\ldots, f_I`` be the frequencies in the various bins of the
histogram and ``I`` the number of bins. With ``N = \\sum_{i=1}^{I}f_i``, let
``p_i = \\frac{f_i}{N}`` (``i = 1, \\ldots, I``) denote the probability
distribution of gray levels. From this distribution one derives two additional
distributions. The first defined for discrete values ``1`` to ``s`` and the
other, from ``s+1`` to ``I``. These distributions are

```math
A: \\frac{p_1}{P_s}, \\frac{p_2}{P_s}, \\ldots, \\frac{p_s}{P_s}
\\quad \\text{and} \\quad
B: \\frac{p_{s+1}}{1-P_s}, \\ldots, \\frac{p_n}{1-P_s}
\\quad \\text{where} \\quad
P_s = \\sum_{i=1}^{s}p_i.
```
The entropies associated with each distribution are as follows:

```math
H(A) = \\ln(P_s) + \\frac{H_s}{P_s}
```
```math
H(B) = \\ln(1-P_s) + \\frac{H_n-H_s}{1-P_s}
```
```math
\\quad \\text{where} \\quad
H_s = -\\sum_{i=1}^{s}p_i\\ln{p_i}
\\quad \\text{and} \\quad
H_n = -\\sum_{i=1}^{I}p_i\\ln{p_i}.
```

Combining these two entropy functions we have

```math
\\psi(s) = \\ln(P_s(1-P_s)) + \\frac{H_s}{P_s} + \\frac{H_n-H_s}{1-P_s}.
```
Finding the discrete value ``s`` which maximises the function ``\\psi(s)`` produces
the sought-after threshold value (i.e. the bin which determines the threshold).

See Section 4 of [1] for more details on the derivation of the entropy.

# Options

## Choices for `counts`
You can specify an `AbstractArray` which should be a 1D array of frequencies for
a histogram. You should submit the corresponding `edges` range for the bins of
the histogram. The function will throw an error if it detects that the `edges`
and `counts` have different lengths.

## Choices for `edges`
You can specify an `AbstractRange` which should be the corresponding range for
the bins of the histogram array passed into `counts`.

# Example

```julia

using TestImages, Images

img = testimage("cameraman")
# building a histogram with 256 bins
edges, counts = build_histogram(img, 256)
#=
  The `counts` array stores at index 0 the frequencies that were below the
  first bin edge. Since we are seeking a threshold over the interval
  partitioned by `edges` we need to discard the first bin in `counts`
  so that the dimensions of `edges` and `counts` match.
=#
find_threshold(EntropyThresholding(), counts[1:end], edges)
```

# References
[1] J. N. Kapur, P. K. Sahoo, and A. K. C. Wong, “A new method for gray-level picture thresholding using the entropy of the histogram,” *Computer Vision, Graphics, and Image Processing*, vol. 29, no. 1, p. 140, Jan. 1985.[doi:10.1016/s0734-189x(85)90156-2](https://doi.org/10.1016/s0734-189x%2885%2990156-2)
"""
function find_threshold(::Entropy, counts::AbstractArray, edges::AbstractRange)
    if length(edges) != length(counts)
        error("the lengths of edges and counts must match")
    end
    nbins = length(counts)
    pdf = counts ./ sum(counts)
    Ψₘₐₓ = 0
    Hₛ = H(pdf)
    Hₙ = -1*Hₛ[nbins]
    Pₛ = cumsum(pdf)
    # Pₛ might exceed 1 for float-point numerical stability issue (issue #35)
    Pₛ = _maybe_inplace_clamp01(Pₛ)

    # Here we use a trivial findmax implementation so as to avoid additional memory allocation.
    # Generally nbins are small number (e.g., 256) so fancy findmax strategy might give less
    # performance boost.
    threshold_bin = 1
    @inbounds for i in 1:nbins
        Pᵢ = Pₛ[i]
        Hᵢ = -1*Hₛ[i]
        Ψᵢ = log(Pᵢ*(1-Pᵢ)) + (Hᵢ/Pᵢ) + ((Hₙ-Hᵢ)/(1-Pᵢ))
        if Ψᵢ > Ψₘₐₓ
            Ψₘₐₓ = Ψᵢ
            threshold_bin = i
        end
    end
    edges[threshold_bin]
end

if VERSION > v"1.5"
    function H(pdf::AbstractArray{T}) where T
        f(x) = x == 0 ? zero(T) : convert(T, x*log(x)) # To avoid NaN when x == 0
        cumsum(f(x) for x in pdf)
    end
else
    function H(pdf::AbstractArray{T}) where T
        f(x) = x == 0 ? zero(T) : convert(T, x*log(x)) # To avoid NaN when x == 0
        cumsum(f.(pdf))
    end
end

# We can't assure that all array types support setindex!, so here we only optimize
# for the `Vector`, which is `counts[1:end]`. `SubArray` may or may not be writable.
_maybe_inplace_clamp01(A::AbstractArray{T}) where T = clamp.(A, zero(T), one(T))
_maybe_inplace_clamp01(A::Vector{T}) where T = clamp!(A, zero(T), one(T))


text:
# http://huntception.com/puzzle/1_2_3/solution/
@testset "1, 2, 3" begin
    f = best_feature(["season", "saveup", "ecowas", "ignore", "sluice", "hosni", "inbed", "barbeau", "museum", "tobiah", "unsew", "dolce", "anaphia", "teenage"])
    @test f.description == "has 3 unique consonants"
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

struct __enum_Version <: ProtoEnum
    V1_0::Int32
    V1_1::Int32
    __enum_Version() = new(0,1)
end #struct __enum_Version
const Version = __enum_Version()

struct __enum_Team <: ProtoEnum
    Self::Int32
    Opponent::Int32
    Nobody::Int32
    __enum_Team() = new(0,1,2)
end #struct __enum_Team
const Team = __enum_Team()

struct __enum_ControlType <: ProtoEnum
    Continue::Int32
    Reset::Int32
    __enum_ControlType() = new(0,1)
end #struct __enum_ControlType
const ControlType = __enum_ControlType()

mutable struct Vector2 <: ProtoType
    x::Float32
    y::Float32
    Vector2(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Vector2

mutable struct TeamInfo <: ProtoType
    team_name::AbstractString
    version::Int32
    TeamInfo(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct TeamInfo

mutable struct Ball <: ProtoType
    position::Vector2
    Ball(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Ball

mutable struct Wheel <: ProtoType
    left_speed::Float32
    right_speed::Float32
    Wheel(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Wheel

mutable struct Robot <: ProtoType
    position::Vector2
    rotation::Float32
    wheel::Wheel
    Robot(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Robot

mutable struct Field <: ProtoType
    self_robots::Base.Vector{Robot}
    opponent_robots::Base.Vector{Robot}
    ball::Ball
    tick::Int32
    Field(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Field

mutable struct Placement <: ProtoType
    robots::Base.Vector{Robot}
    ball::Ball
    Placement(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Placement

mutable struct ControlInfo <: ProtoType
    command::Int32
    ControlInfo(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ControlInfo

export Version, Team, ControlType, Vector2, TeamInfo, Ball, Wheel, Robot, Field, Placement, ControlInfo


text:
@testset "Analysis tools" begin
    # minimalistic test
    @testset "Relative Energy regain" begin
        x = [1f0 -2; 4 5; 7 8 ]
        @test DeconvOptim.relative_energy_regain(x, x .* 1) == (Float32[0.0, 0.333, 0.5, 0.601], Float32[1.0, 1.0, 1.0, 1.0])
        @test DeconvOptim.relative_energy_regain(x, x .* 0.5) == (Float32[0.0, 0.333, 0.5, 0.601], Float32[0.75, 0.75, 0.75, 0.75])
    end

    # minimalistic test
    @testset "Normalized Cross Correlation" begin
        x = [1f0 -2; 4 5; 7 8 ]

        @test DeconvOptim.normalized_cross_correlation(x, x) == 1.0f0
    end

end


text:
using Test, DataFrames, CSV
df = CSV.read("TOKCC - Open.xlsx - DECEMBER 2021.csv", DataFrame)
Score2 = TOKCC.score.(df.KB, df.Reps, df.BW)
@test all(isapprox.(df.Score, Score2; atol=0.1))

max_stars_repo_name: wjholden/TOKCC

max_issues_repo_name: wjholden/TOKCC

max_forks_repo_name: wjholden/TOKCC


text:
function solve_lyapunov(A,Q,guess = I;iter_lim=100,epsilon=1e-14)
    iter = 1
    err = 1
    while iter <= iter_lim && err > epsilon
        new_guess = A*guess*A' + Q
        err = maximum(abs.(new_guess - guess))
        guess = new_guess
        iter += 1
    end
    return(guess)
end

using LinearAlgebra

function solve_lyapunov_vec(A,Q)
    aux = (I - kron(A,A))
    vec_ans = \(aux,vec(Q))
    res = reshape(vec_ans,size(Q))
    return res
end

using ToeplitzMatrices

function build_variance(g,h,Q,T) #only works for a single series
    sig_x = solve_lyapunov_vec(h,Q)
    y_autocov(t) = g*h*sig_x*h'^(t)*g'
    vec_st = map(x->first(y_autocov(x)),1:(T))
    res = Toeplitz(vec_st,vec_st)
    return res
end

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
function printObj3DTree(scene::Scene; world::Union{Nothing, Object3D}=nothing)
    printObjInfos(world)
    printTree(scene.treeAccVelo)
end

function printAllVisuElements(scene::Scene)
    printTree(scene.allVisuElements)
end

function printTree(tree)
    for obj in tree
        printObjInfos(obj)
    end
end

printObjInfos(obj::Nothing) = nothing
function printObjInfos(obj::Object3D)
    println("parent   = ", Modia3D.fullName(obj.parent))
    println("Object3D = ", Modia3D.fullName(obj))
    println("children:")
    for child in obj.children
        println("  child = ", Modia3D.fullName(child))
    end
    println("joint   = ", obj.joint)
    println(" ")
end

function printScene(scene)
    println("\nScene")
    tree = scene.options.useOptimizedStructure ? scene.treeAccVelo : scene.tree
    if hasParent(tree[1])
        println("  ", tree[1].parent.path)
    end
    for obj in tree
        println("  ", obj.path)
        if hasParent(obj)
            println("    parent = ", obj.parent.path)
        end
        println("    r_rel = ", obj.r_rel)
        print("    joint")
        if obj.joint.path != ""
            println(" (= ", obj.joint.path, ")")
        else
            println()
        end
        println("      kind       = ", obj.joint.kind)
        println("      ndof       = ", obj.joint.ndof)
        println("      canCollide = ", obj.joint.canCollide)
        if obj.joint.kind == RevoluteKind
            println("      phi        = ", obj.joint.specific.phi)
            println("      w          = ", obj.joint.specific.w)
        elseif obj.joint.kind == PrismaticKind
            println("      s          = ", obj.joint.specific.s)
            println("      v          = ", obj.joint.specific.v)
        elseif obj.joint.kind == FreeMotionKind || obj.joint.kind == AbsoluteFreeMotionKind
            println("      r          = ", obj.joint.specific.r)
            println("      rot        = ", obj.joint.specific.rot)
            println("      v          = ", obj.joint.specific.v)
            println("      w          = ", obj.joint.specific.w)
        end

        hasMassProperties =  scene.options.useOptimizedStructure && objectHasMass(obj) ||
                             !scene.options.useOptimizedStructure && featureHasMass(obj)
        if hasMassProperties
            massProperties::Shapes.InternalMassProperties = scene.options.useOptimizedStructure ? obj.massProperties : obj.feature.massProperties
            println("    massProperties")
            println("      m   = ", massProperties.m)
            println("      rCM = ", massProperties.rCM)
        end
    end
    println()
end


text:
#------------------------------------------------------------------------------

function truncate_matrix(M::MatElem{<: Generic.AbsSeriesElem}, prec::Int)
    return map(e -> truncate(e, prec), M)
end

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

function matrix_set_precision!(M::MatElem{<: Generic.AbsSeriesElem}, prec::Int)
    map(e -> set_precision!(e, prec), M)
end

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

function ps_matrix_const_term(M::MatElem{<: Generic.AbsSeriesElem})
    return map(e -> coeff(e, 0), M)
end

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

"""
    _matrix_inv_newton_iteration(M, Minv)

Performs a single step of Newton iteration for inverting `M` with `Minv` being a partial result
"""
function _matrix_inv_newton_iteration(M::MatElem{T}, Minv::MatElem{T}) where T <: Generic.AbsSeriesElem{<: Generic.FieldElem}
    return 2 * Minv - Minv * M * Minv
end

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

"""
    ps_matrix_inv(M, prec)

Input:
- `M` - a square matrix with entries in a univariate power series ring
      it is assumed that M(0) is invertible and all entries having the same precision
- `prec` - an integer, precision, if `-1` then defaults to precision of `M`

Output:
- the inverse of `M` computed up to `prec`
"""
function ps_matrix_inv(M::MatElem{<: Generic.AbsSeriesElem{<: Generic.FieldElem}}, prec::Int=-1)
    const_term = ps_matrix_const_term(M)
    prec = (prec == -1) ? precision(M[1, 1]) : prec
    power_series_ring = base_ring(parent(M))
    result = map(a -> power_series_ring(a), Base.inv(const_term))
    cur_prec = 1
    while cur_prec < prec
        result = _matrix_inv_newton_iteration(M, result)
        cur_prec *= 2
    end
    return result
end

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

"""
    ps_diff(ps)

Input:
- `ps` - (absolute capped) unvariate power series

Output:
- the derivative of `ps`
"""
function ps_diff(ps::Generic.AbsSeriesElem{<: Generic.RingElem})
    result = zero(parent(ps))
    set_precision!(result, precision(ps))
    for exp in 1:(precision(ps) - 1)
        setcoeff!(result, exp - 1, coeff(ps, exp) * exp)
    end
    return result
end

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

"""
    ps_integrate(ps)

Input:
- `ps` - (absolute capped) unvariate power series
Output:
- the integral of `ps` without constant term
"""
function ps_integrate(ps::Generic.AbsSeriesElem{<: Generic.FieldElem})
    result = zero(parent(ps))
    set_precision!(result, precision(ps) + 1)
    for exp in 0:(precision(ps) - 1)
        setcoeff!(result, exp + 1, coeff(ps, exp) // (exp + 1))
    end
    return result
end

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

"""
    ps_matrix_log(M)

Input:
- M - a square matrix with entries in a univariate power series ring
      it is assumed that `M(0)` is the identity
Output:
- the natural log of `M`
"""
function ps_matrix_log(M::MatElem{<: Generic.AbsSeriesElem{<: Generic.FieldElem}})
    const_term = ps_matrix_const_term(M)
    if const_term != one(parent(const_term))
       throw(Base.DomainError("Constant term must be the identity matrix"))
    end
    return map(
        e -> ps_integrate(e),
        map(e -> ps_diff(e), M - one(parent(M))) *
        ps_matrix_inv(M)
    )
end

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

function _matrix_homlinear_de_newton_iteration(
        A::MatElem{<: Generic.AbsSeriesElem{T}},
        Y::MatElem{<: Generic.AbsSeriesElem{T}},
        Z::MatElem{<: Generic.AbsSeriesElem{T}},
        cur_prec::Int
    ) where T <: Generic.FieldElem
    Yprime = map(ps_diff, Y)
    Z = Z + truncate_matrix(Z * (one(parent(A)) - Y * Z), cur_prec)
    Y = Y - truncate_matrix(Y * map(ps_integrate, Z * (Yprime - truncate_matrix(A, 2 * cur_prec - 1) * Y)), 2 * cur_prec)
    return (Y, Z)
end

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

"""
    ps_matrix_homlinear_de(A, Y0, prec)

Input:
- `A` - a square matrix with entries in a univariate power series ring
- `Y0` - a square invertible matrix over the base field
Output:
- matrix `Y` such that `Y' = AY` up to precision of `A - 1` and `Y(0) = Y0`
"""
function ps_matrix_homlinear_de(
        A::MatElem{<: Generic.AbsSeriesElem{T}},
        Y0::MatElem{<: T},
        prec::Int=-1
    ) where T <: Generic.FieldElem
    prec = (prec == -1) ? precision(A[1, 1]) : prec
    ps_ring = base_ring(parent(A))
    cur_prec = 1
    Y = (one(parent(A)) + gen(ps_ring) * truncate_matrix(A, cur_prec)) * map(e -> truncate(ps_ring(e), cur_prec), Y0)
    Z = map(e -> truncate(ps_ring(e), cur_prec), Base.inv(Y0))
    while cur_prec < prec
        matrix_set_precision!(Y, 2 * cur_prec)
        matrix_set_precision!(Z, cur_prec)
        Y, Z = _matrix_homlinear_de_newton_iteration(A, Y, Z, cur_prec)
        cur_prec *= 2
    end
    return Y, Z
end

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

function _variation_of_constants(
        A::MatElem{<: Generic.AbsSeriesElem{T}},
        B::MatElem{<: Generic.AbsSeriesElem{T}},
        Yh::MatElem{<: Generic.AbsSeriesElem{T}},
        Zh::MatElem{<: Generic.AbsSeriesElem{T}},
        Y0::MatElem{<: T},
        prec::Int
    ) where T <: Generic.FieldElem
    Zh += truncate_matrix(Zh * (one(parent(A)) - Yh * Zh), prec)
    Y_particular = truncate_matrix(Yh * map(ps_integrate, Zh * B), prec)
    return Y_particular + Yh * map(e -> parent(A[1, 1])(e), Y0)
end

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

"""
    ps_matrix_linear_de(A, B, Y0, prec)

Input:
- `A`, `B` - square matrices with entries in a univariate power series ring
- `Y0` - a matrix over the base field with the rows number the same as `A`

Output:
- matrix `Y` such that `Y' = AY + B` up to precision of `A - 1` and `Y(0) = Y0`
"""
function ps_matrix_linear_de(
        A::MatElem{<: Generic.AbsSeriesElem{T}},
        B::MatElem{<: Generic.AbsSeriesElem{T}},
        Y0::MatElem{<: T},
        prec::Int=-1
    ) where T <: Generic.FieldElem
    prec = (prec == -1) ? precision(A[1, 1]) : prec
    n = nrows(A)
    identity = one(AbstractAlgebra.MatrixSpace(base_ring(parent(Y0)), n, n))
    Yh, Zh = ps_matrix_homlinear_de(A, identity, prec)
    matrix_set_precision!(Zh, prec)
    return _variation_of_constants(A, B, Yh, Zh, Y0, prec)
end

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

"""
    ps_ode_solution(equations, ic, inputs, prec)

Input:
- `equations` - a system of the form ``A(x, u, mu)x' - B(x, u, mu) = 0``,
                where A is a generically nonsingular square matrix. *Assumption*: `A` is nonzero at zero
- `ic` - initial conditions for x's (dictionary)
- `inputs` - power series for inputs represented as arrays (dictionary)
- `prec` - precision of the solution

Output:
- power series solution of the system
"""
function ps_ode_solution(
        equations::Array{P, 1},
        ic::Dict{P, T},
        inputs::Dict{P, Array{T, 1}},
        prec::Int
    ) where {T <: Generic.FieldElem, P <: MPolyElem{T}}
    n = length(equations)
    ring = parent(equations[1])
    S = AbstractAlgebra.MatrixSpace(ring, n, n)
    Sv = AbstractAlgebra.MatrixSpace(ring, n, 1)
    Svconst = AbstractAlgebra.MatrixSpace(base_ring(ring), n, 1)
    eqs = Sv(equations)

    x_vars = filter(v -> ("$(v)_dot" in map(string, gens(ring))), gens(ring))
    x_vars = [x for x in x_vars]
    x_dot_vars = [str_to_var(var_to_str(x) * "_dot", ring) for x in x_vars]

    Jac_dots = S([derivative(p, xd) for p in equations, xd in x_dot_vars])
    Jac_xs = S([derivative(p, x) for p in equations, x in x_vars])

    ps_ring, t = PowerSeriesRing(base_ring(ring), prec, "t"; model=:capped_absolute)
    solution = Dict()
    for (u, coeffs) in inputs
        solution[u] = sum(coeffs[i] * t^(i - 1) for i in 1:length(coeffs))
    end
    for x in x_vars
        solution[x] = ps_ring(ic[x])
    end
    for xd in x_dot_vars
        solution[xd] = ps_ring(0)
        set_precision!(solution[xd], 1)
    end

    cur_prec = 1
    while cur_prec < prec
        @debug "\t Computing power series solution, currently at precision $cur_prec"
        new_prec = min(prec, 2 * cur_prec)
        for i in 1:length(x_vars)
            set_precision!(solution[x_vars[i]], new_prec)
            set_precision!(solution[x_dot_vars[i]], new_prec)
        end
        eval_point = [solution[v] for v in gens(ring)]
        map(ps -> set_precision!(ps, 2 * cur_prec), eval_point)
        eqs_eval = map(p -> evaluate(p, eval_point), eqs)
        J_eval = map(p -> evaluate(p, eval_point), Jac_xs)
        Jd_eval = map(p -> evaluate(p, eval_point), Jac_dots)
        Jd_inv = ps_matrix_inv(Jd_eval, cur_prec)

        X_err = ps_matrix_linear_de(-Jd_inv * J_eval, -Jd_inv * eqs_eval, zero(Svconst))
        for i in 1:length(x_vars)
            solution[x_vars[i]] = solution[x_vars[i]] + X_err[i]
            solution[x_dot_vars[i]] = ps_diff(solution[x_vars[i]])
        end
        cur_prec = new_prec
    end

    return solution
end

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

function ps_ode_solution(
        equations::Array{P, 1},
        ic::Dict{P, Int},
        inputs::Dict{P, Array{Int, 1}},
        prec::Int
    ) where P <: MPolyElem{<: Generic.FieldElem}
    bring = base_ring(parent(equations[1]))
    ps_ode_solution(
        equations,
        Dict(x => bring(v) for (x, v) in ic),
        Dict(u => map(v -> bring(v), vv) for (u, vv) in inputs),
        prec
    )
end


text:
@testset "native" begin

include("definitions/native.jl")

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

@testset "Compilation" begin
    kernel() = nothing

    output = native_code_execution(kernel, (); validate=false)
    @test occursin("kernel", output[2])
    @test isempty(output[3])
    @test isempty(output[4])

    @testset "Undefined Functions" begin
        function undef_fn()
            ccall("extern somefunc", llvmcall, Cvoid, ())
            nothing
        end

        output = native_code_execution(undef_fn, (); validate=false)
        @test length(output[3]) == 1
        @test output[3][1] == "somefunc"
    end

    @testset "Undefined Globals" begin
        @generated function makegbl(::Val{name}, ::Type{T}, ::Val{isext}) where {name,T,isext}
            JuliaContext() do ctx
                T_gbl = convert(LLVMType, T, ctx)
                T_ptr = convert(LLVMType, Ptr{T}, ctx)
                llvm_f, _ = create_function(T_ptr)
                mod = LLVM.parent(llvm_f)
                gvar = GlobalVariable(mod, T_gbl, string(name))
                isext && extinit!(gvar, true)
                Builder(ctx) do builder
                    entry = BasicBlock(llvm_f, "entry", ctx)
                    position!(builder, entry)
                    result = ptrtoint!(builder, gvar, T_ptr)
                    ret!(builder, result)
                end
                call_function(llvm_f, Ptr{T})
            end
        end
        function undef_gbl()
            ext_ptr = makegbl(Val(:someglobal), Int64, Val(true))
            Base.unsafe_store!(ext_ptr, 1)
            ptr = makegbl(Val(:otherglobal), Float32, Val(false))
            Base.unsafe_store!(ptr, 2f0)
            nothing
        end

        output = native_code_execution(undef_gbl, ())
        @test length(output[4]) == 2
        @test output[4][1].name == "someglobal"
        @test eltype(output[4][1].type) isa LLVM.IntegerType
        @test output[4][1].external
        @test output[4][2].name == "otherglobal"
        @test eltype(output[4][2].type) isa LLVM.LLVMFloat
        @test !output[4][2].external
    end
end

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

@testset "IR" begin

@testset "basic reflection" begin
    valid_kernel() = return
    invalid_kernel() = 1

    ir = sprint(io->native_code_llvm(io, valid_kernel, Tuple{}; optimize=false, dump_module=true))

    # module should contain our function + a generic call wrapper
    @test occursin(r"define\ .* void\ @.*julia_valid_kernel.*\(\)"x, ir)
    @test !occursin("define %jl_value_t* @jlcall_", ir)

    # there should be no debug metadata
    @test !occursin("!dbg", ir)

    @test native_code_llvm(devnull, invalid_kernel, Tuple{}) == nothing
    @test_throws KernelError native_code_llvm(devnull, invalid_kernel, Tuple{}; kernel=true) == nothing
end

@testset "unbound typevars" begin
    invalid_kernel() where {unbound} = return
    @test_throws KernelError native_code_llvm(devnull, invalid_kernel, Tuple{})
end

@testset "child functions" begin
    # we often test using `@noinline sink` child functions, so test whether these survive
    @noinline child(i) = sink(i)
    parent(i) = child(i)

    ir = sprint(io->native_code_llvm(io, parent, Tuple{Int}))
    @test occursin(r"call .+ @julia_child_", ir)
end

@testset "sysimg" begin
    # bug: use a system image function

    function foobar(a,i)
        Base.pointerset(a, 0, mod1(i,10), 8)
    end

    ir = sprint(io->native_code_llvm(io, foobar, Tuple{Ptr{Int},Int}))
    @test !occursin("jlsys_", ir)
end

@testset "tracked pointers" begin
    function kernel(a)
        a[1] = 1
        return
    end

    # this used to throw an LLVM assertion (#223)
    native_code_llvm(devnull, kernel, Tuple{Vector{Int}}; kernel=true)
end

if VERSION >= v"1.0.2"
@testset "CUDAnative.jl#278" begin
    # codegen idempotency
    # NOTE: this isn't fixed, but surfaces here due to bad inference of checked_sub
    # NOTE: with the fix to print_to_string this doesn't error anymore,
    #       but still have a test to make sure it doesn't regress
    native_code_llvm(devnull, Base.checked_sub, Tuple{Int,Int}; optimize=false)
    native_code_llvm(devnull, Base.checked_sub, Tuple{Int,Int}; optimize=false)

    # breaking recursion in print_to_string makes it possible to compile
    # even in the presence of the above bug
    native_code_llvm(devnull, Base.print_to_string, Tuple{Int,Int}; optimize=false)
end
end

@testset "LLVM D32593" begin
    @eval struct D32593_struct
        foo::Float32
        bar::Float32
    end

    D32593(ptr) = unsafe_load(ptr).foo

    native_code_llvm(devnull, D32593, Tuple{Ptr{D32593_struct}})
end

@testset "Julia-level throw lowering" begin
    kernel(ptr) = (unsafe_store!(ptr, sqrt(unsafe_load(ptr))); nothing)

    native_code_execution(kernel, Tuple{Ptr{Float32}})
end

end

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

@testset "assembly" begin

@testset "basic reflection" begin
    valid_kernel() = return
    invalid_kernel() = 1

    @test native_code_native(devnull, valid_kernel, Tuple{}) == nothing
    @test native_code_native(devnull, invalid_kernel, Tuple{}) == nothing
    @test_throws KernelError native_code_native(devnull, invalid_kernel, Tuple{}; kernel=true)
end

@testset "idempotency" begin
    # bug: generate code twice for the same kernel (jl_to_ptx wasn't idempotent)

    kernel() = return
    native_code_native(devnull, kernel, Tuple{})
    native_code_native(devnull, kernel, Tuple{})
end

@testset "compile for host after gpu" begin
    # issue #11: re-using host functions after GPU compilation
    @noinline child(i) = sink(i+1)

    function fromhost()
        child(10)
    end

    function fromptx()
        child(10)
        return
    end

    native_code_native(devnull, fromptx, Tuple{})
    @test fromhost() == 11
end

end

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

@testset "errors" begin

# some validation happens in the emit_function hook, which is called by code_llvm

@testset "base intrinsics" begin
    foobar(i) = sin(i)

    # NOTE: we don't use test_logs in order to test all of the warning (exception, backtrace)
    logs, _ = Test.collect_test_logs(min_level=Info) do
        withenv("JULIA_DEBUG" => nothing) do
            native_code_llvm(devnull, foobar, Tuple{Int})
        end
    end
    @test length(logs) == 1
    record = logs[1]
    @test record.level == Base.CoreLogging.Warn
    @test record.message == "calls to Base intrinsics might be GPU incompatible"
    @test haskey(record.kwargs, :exception)
    err,bt = record.kwargs[:exception]
    err_msg = sprint(showerror, err)
    @test occursin(Regex("You called sin(.+) in Base.Math .+, maybe you intended to call sin(.+) in $TestRuntime .+ instead?"), err_msg)
    bt_msg = sprint(Base.show_backtrace, bt)
    @test occursin("[1] sin", bt_msg)
    @test occursin(r"\[2\] .+foobar", bt_msg)
end

# some validation happens in `compile`

@eval Main begin
struct CleverType{T}
    x::T
end
Base.unsafe_trunc(::Type{Int}, x::CleverType) = unsafe_trunc(Int, x.x)
end

@testset "non-isbits arguments" begin
    foobar(i) = (sink(unsafe_trunc(Int,i)); return)

    @test_throws_message(KernelError,
                         native_code_execution(foobar, Tuple{BigInt})) do msg
        occursin("passing and using non-bitstype argument", msg) &&
        occursin("BigInt", msg)
    end

    # test that we can handle abstract types
    @test_throws_message(KernelError,
                         native_code_execution(foobar, Tuple{Any})) do msg
        occursin("passing and using non-bitstype argument", msg) &&
        occursin("Any", msg)
    end

    @test_throws_message(KernelError,
                         native_code_execution(foobar, Tuple{Union{Int32, Int64}})) do msg
        occursin("passing and using non-bitstype argument", msg) &&
        occursin("Union{Int32, Int64}", msg)
    end

    @test_throws_message(KernelError,
                         native_code_execution(foobar, Tuple{Union{Int32, Int64}})) do msg
        occursin("passing and using non-bitstype argument", msg) &&
        occursin("Union{Int32, Int64}", msg)
    end

    # test that we get information about fields and reason why something is not isbits
    @test_throws_message(KernelError,
                         native_code_execution(foobar, Tuple{CleverType{BigInt}})) do msg
        occursin("passing and using non-bitstype argument", msg) &&
        occursin("CleverType", msg) &&
        occursin("BigInt", msg)
    end
end

@testset "invalid LLVM IR" begin
    foobar(i) = println(i)

    @test_throws_message(InvalidIRError,
                         native_code_execution(foobar, Tuple{Int})) do msg
        occursin("invalid LLVM IR", msg) &&
        occursin(GPUCompiler.RUNTIME_FUNCTION, msg) &&
        occursin("[1] println", msg) &&
        occursin(r"\[2\] .+foobar", msg)
    end
end

@testset "invalid LLVM IR (ccall)" begin
    foobar(p) = (unsafe_store!(p, ccall(:time, Cint, ())); nothing)

    @test_throws_message(InvalidIRError,
                         native_code_execution(foobar, Tuple{Ptr{Int}})) do msg
        occursin("invalid LLVM IR", msg) &&
        occursin(GPUCompiler.POINTER_FUNCTION, msg) &&
        occursin(r"\[1\] .+foobar", msg)
    end
end

@testset "delayed bindings" begin
    kernel() = (undefined; return)

    @test_throws_message(InvalidIRError,
                         native_code_execution(kernel, Tuple{})) do msg
        occursin("invalid LLVM IR", msg) &&
        occursin(GPUCompiler.DELAYED_BINDING, msg) &&
        occursin("use of 'undefined'", msg) &&
        occursin(r"\[1\] .+kernel", msg)
    end
end

@testset "dynamic call (invoke)" begin
    @eval @noinline nospecialize_child(@nospecialize(i)) = i
    kernel(a, b) = (unsafe_store!(b, nospecialize_child(a)); return)

    @test_throws_message(InvalidIRError,
                         native_code_execution(kernel, Tuple{Int,Ptr{Int}})) do msg
        occursin("invalid LLVM IR", msg) &&
        occursin(GPUCompiler.DYNAMIC_CALL, msg) &&
        occursin("call to nospecialize_child", msg) &&
        occursin(r"\[1\] .+kernel", msg)
    end
end

@testset "dynamic call (apply)" begin
    func() = println(1)

    @test_throws_message(InvalidIRError,
                         native_code_execution(func, Tuple{})) do msg
        occursin("invalid LLVM IR", msg) &&
        occursin(GPUCompiler.DYNAMIC_CALL, msg) &&
        occursin("call to println", msg) &&
        occursin("[2] func", msg)
    end
end

end

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

end


text:
module ConvertBNG

import DelimitedFiles

export convert_bng
export convert_lonlat
export convert_bng_nothread
export convert_lonlat_nothread
export convert_osgb36
export convert_etrs89
export convert_osgb36_to_ll
export convert_osgb36_to_etrs89
export convert_to_ll
export get_ostn_ref
export ostn15_shifts
export check
export convert_etrs89_to_ll
export convert_etrs89_to_osgb36

const AIRY_1830_SEMI_MAJOR = 6377563.396
const AIRY_1830_SEMI_MINOR = 6356256.909
const GRS80_SEMI_MAJOR = 6378137.000
const GRS80_SEMI_MINOR = 6356752.3141
const RAD = π / 180.0
# lon and lat of true origin
const λ₀ = RAD * -2.0
const ϕ₀ = RAD * 49.0
# Easting and Northing of origin
const E₀ = 400000.0
const N₀ = -100000.0
# convergence factor
const F₀ = 0.9996012717
# Coordinate bounds for eastings, northings
const MAX_EASTING = 700000.000
const MAX_NORTHING = 1250000.000
# Coordinate bounds for lon, lat
const MIN_LONGITUDE = -8.5790
const MAX_LONGITUDE = 1.7800
const MIN_LATITUDE = 49.9220
const MAX_LATITUDE = 60.8400
# Const container for ostn15 shift data
const OSTN15_DATA = Ref{Array{Float64,2}}()

"""
Import ostn15 shifts data file into array.
Line position in the file becomes array index and so it is critical that the
raw datafile is not midified.
"""
function __init__()
    # Import ostn15 lookup array
    ostn15_file = joinpath(@__DIR__, "..", "data", "ostn15.dat")
    global OSTN15_DATA[] = DelimitedFiles.readdlm(ostn15_file, ',', Float64)
end

"""
Bounds checking for input values
Return input value if within specified bounds, else return NaN
"""
function check(to_check::T1, bounds::Tuple{T2,T2}) where {T1<:Real, T2<:Real}
    if bounds[1] <= to_check <= bounds[2]
        return to_check
    else
        return NaN
    end
end

"""
Round a float coordinate pair to nearest mm
"""
function round_to_mm(x::T, y::T) where {T<:AbstractFloat}
    new_x = round(x * 1000.0, digits=0) / 1000.0
    new_y = round(y * 1000.0, digits=0) / 1000.0
    return new_x, new_y
end

"""
Round a float coordinate pair to eight decimal places
"""
function round_to_eight(x::T, y::T) where {T<:AbstractFloat}
    new_x = round(x * 100000000.0, digits=0) / 100000000.0
    new_y = round(y * 100000000.0, digits=0) / 100000000.0
    return new_x, new_y
end

"""
Intermediate calculation used for lon, lat to ETRS89 and reverse conversion
"""
@fastmath function compute_m(ϕ::T, b::T, n::T) where {T<:Real}
    n² = n ^ 2
    n³ = n ^ 3
    result = b * F₀ * (
        (1.0 + n + 5.0 / 4.0 * n² + 5.0 / 4.0 * n³) * (ϕ - ϕ₀) -
        (3.0 * n + 3.0 * n² + 21.0 / 8.0 * n³) * sin(ϕ - ϕ₀) * cos(ϕ + ϕ₀) +
        (15.0 / 8.0 * n² + 15.0 / 8.0 * n³) * sin(2.0 * (ϕ - ϕ₀)) * cos(2.0 * (ϕ + ϕ₀)) -
        35.0 / 24.0 * n³ * sin(3.0 * (ϕ - ϕ₀)) * cos(3.0 * (ϕ + ϕ₀))
    )
    return result
end

"""
Perform Longitude, Latitude to ETRS89 conversion
"""
@fastmath function convert_etrs89(longitude::T, latitude::T) where {T<:Real}
    # Input is restricted to the UK bounding box
    # Convert bounds-checked input to degrees, or return an Err
    λ = deg2rad(check(longitude, (MIN_LONGITUDE, MAX_LONGITUDE)))
    ϕ = deg2rad(check(latitude, (MIN_LATITUDE, MAX_LATITUDE)))
    # ellipsoid squared eccentricity constant
    e² = (GRS80_SEMI_MAJOR ^ 2 - GRS80_SEMI_MINOR ^ 2) / GRS80_SEMI_MAJOR ^ 2
    n = (GRS80_SEMI_MAJOR - GRS80_SEMI_MINOR) / (GRS80_SEMI_MAJOR + GRS80_SEMI_MINOR)

    sin²ϕ = sin(ϕ) ^ 2
    ν = GRS80_SEMI_MAJOR * F₀ * (1.0 - e² * sin²ϕ) ^ -0.5
    ρ = GRS80_SEMI_MAJOR * F₀ * (1.0 - e²) * (1.0 - e² * sin²ϕ) ^ -1.5
    η² = ν / ρ - 1.0

    m = compute_m(ϕ, GRS80_SEMI_MINOR, n)

    cosϕ = cos(ϕ)
    cos³ϕ = cosϕ ^ 3
    cos⁵ϕ = cosϕ ^ 5
    sinϕ = sin(ϕ)
    tanϕ = tan(ϕ)
    tan²ϕ = tanϕ ^ 2
    tan⁴ϕ = tanϕ ^ 4

    I = m + N₀
    II = ν / 2.0 * sinϕ * cosϕ
    III = ν / 24.0 * sinϕ * cos³ϕ * (5.0 - tan²ϕ + 9.0 * η²)
    IIIA = ν / 720.0 * sinϕ * cos⁵ϕ * (61.0 - 58.0 * tan²ϕ + tan⁴ϕ)
    IV = ν * cosϕ
    V = ν / 6.0 * cos³ϕ * (ν / ρ - tan²ϕ)
    VI = ν / 120.0 * cos⁵ϕ * (5.0 - 18.0 * tan²ϕ + tan⁴ϕ + 14.0 * η² - 58.0 * tan²ϕ * η²)
    l = λ - λ₀
    north = I + II * l ^ 2 + III * l ^ 4 + IIIA * l ^ 6
    east = E₀ + IV * l + V * l ^ 3 + VI * l ^ 5

    return east, north
end

"""
Easting and Northing to Lon, Lat conversion using a Helmert transform
Note that either GRS80 or Airy 1830 ellipsoids can be passed
"""
@fastmath function convert_to_ll(eastings::T, northings::T, ℓ_a::T, ℓ_b::T) where {T<:Real}
    # ellipsoid squared eccentricity constant
    a = ℓ_a
    b = ℓ_b
    e² = (a ^ 2 - b ^ 2) / a ^ 2
    n = (a - b) / (a + b)

    dN = check(northings, (0.0, MAX_NORTHING)) - N₀
    ϕ = ϕ₀ + dN / (a * F₀)
    m = compute_m(ϕ, b, n)
    while (dN - m) >= 0.00001
        m = compute_m(ϕ, b, n)
        ϕ += (dN - m) / (a * F₀)
    end

    sin²ϕ = sin(ϕ) ^ 2
    tanϕ = tan(ϕ)
    tan²ϕ = tanϕ ^ 2
    tan⁴ϕ = tanϕ ^ 4
    tan⁶ϕ = tan⁴ϕ * tan²ϕ
    secϕ = 1.0 / cos(ϕ)

    ν = a * F₀ * (1.0 - e² * sin²ϕ) ^ -0.5
    ρ = a * F₀ * (1.0 - e²) * (1.0 - e² * sin²ϕ) ^ -1.5
    η² = ν / ρ - 1.0

    VII = tanϕ / (2.0 * ρ * ν)
    VIII = tanϕ / (24.0 * ρ * ν ^ 3) * (5.0 + 3.0 * tan²ϕ + η² - 9.0 * tan²ϕ * η²)
    IX = tanϕ / (720.0 * ρ * ν ^ 5) * (61.0 + 90.0 * tan²ϕ + 45.0 * tan⁴ϕ)
    X = secϕ / ν
    XI = secϕ / (6.0 * ν ^ 3) * (ν / ρ + 2.0 * tan²ϕ)
    XII = secϕ / (120.0 * ν ^ 5) * (5.0 + 28.0 * tan²ϕ + 24.0 * tan⁴ϕ)
    XIIA = secϕ / (5040.0 * ν ^ 7) * (61.0 + 662.0 * tan²ϕ + 1320.0 * tan⁴ϕ + 720.0 * tan⁶ϕ)

    e = check(eastings, (0.000, MAX_EASTING)) - E₀

    ϕ = ϕ - VII * e ^ 2 + VIII * e ^ 4 - IX * e ^ 6
    λ = λ₀ + X * e - XI * e ^ 3 + XII * e ^ 5 - XIIA * e ^ 7

    ϕ = rad2deg(ϕ)
    λ = rad2deg(λ)
    λ, ϕ = round_to_eight(λ, ϕ)
    return λ, ϕ
end

function convert_etrs89_to_ll(E::T, N::T) where {T<:Real}
    # ETRS89 uses the WGS84 / GRS80 ellipsoid constants
    return convert_to_ll(E, N, GRS80_SEMI_MAJOR, GRS80_SEMI_MINOR)
end

"""
Perform Longitude, Latitude to OSGB36 conversion, using [OSTN15] data
"""
function convert_osgb36(longitude::T, latitude::T) where {T<:Real}
    # Convert input to ETRS89
    eastings, northings = convert_etrs89(longitude, latitude)
    # Obtain OSTN15 corrections, and incorporate
    return convert_etrs89_to_osgb36(eastings, northings)
end

function convert_osgb36(lonlat::Array{T,2}) where{T<:Real}
    # Create empty output array
    en_arr = zeros(T, size(lonlat))
    # Transform
    Threads.@threads for i in 1:size(lonlat)[1]
        en_arr[i,:] .= convert_osgb36(lonlat[i,1], lonlat[i,2])
    end
    return en_arr
end

"""
Non-threaded version of convert_osgb36
"""
function convert_osgb36_nothread(lonlat::Array{T,2}) where{T<:Real}
    # Create empty output array
    en_arr = zeros(T, size(lonlat))
    # Transform
    for i in 1:size(lonlat)[1]
        en_arr[i,:] .= convert_osgb36(lonlat[i,1], lonlat[i,2])
    end
    return en_arr
end

"""
Perform ETRS89 to OSGB36 conversion, using [OSTN15] data
"""
function convert_etrs89_to_osgb36(E::T, N::T) where {T<:Real}
    # Obtain OSTN15 corrections, and incorporate
    e_shift, n_shift, _ = ostn15_shifts(E, N)
    e_corr, n_corr = round_to_mm(E + e_shift, N + n_shift)
    return e_corr, n_corr
end

"""
Convert OSGB36 coordinates to Lon, Lat using OSTN15 data
"""
function convert_osgb36_to_ll(E::T, N::T) where {T<:Real}
    # Apply reverse OSTN15 adustments
    x, y = convert_osgb36_to_etrs89(E, N)
    # We've converted to ETRS89, so we need to use the WGS84/ GRS80 ellipsoid constants
    return convert_etrs89_to_ll(x, y)
end

function convert_osgb36_to_ll(en::Array{T,2}) where{T<:Real}
    # Create empty output array
    ll_arr = zeros(T, size(en))
    # Transform
    Threads.@threads for i in 1:size(en)[1]
        ll_arr[i,:] .= convert_osgb36_to_ll(en[i,1], en[i,2])
    end
    return ll_arr
end

"""
Non-threaded version of convert_osgb36_to_ll
"""
function convert_osgb36_to_ll_nothread(en::Array{T,2}) where{T<:Real}
    # Create empty output array
    ll_arr = zeros(T, size(en))
    # Transform
    for i in 1:size(en)[1]
        ll_arr[i,:] .= convert_osgb36_to_ll(en[i,1], en[i,2])
    end
    return ll_arr
end

"""
Convert OSGB36 coordinates to ETRS89 using OSTN15 data
"""
function convert_osgb36_to_etrs89(E::T, N::T) where {T<:Real}
    # Apply reverse OSTN15 adustments
    epsilon = 0.00001
    dx, dy, _ = ostn15_shifts(E, N)
    x, y = E - dx, N - dy
    last_dx, last_dy = dx, dy
    while !isnan(dx)
        dx, dy, _ = ostn15_shifts(x, y)
        x = E - dx
        y = N - dy
        if abs(dx - last_dx) < epsilon && abs(dy - last_dy) < epsilon
            break
        end
        last_dx = dx
        last_dy = dy
    end
    x = E - dx
    y = N - dy
    return x, y
end

"""
Function alias names to align with convert_bng Python library
"""
function convert_bng(lonlat) convert_osgb36(lonlat) end
function convert_bng_nothread(lonlat) convert_osgb36_nothread(lonlat) end
function convert_bng(lon, lat) convert_osgb36(lon, lat) end
function convert_lonlat(EN) convert_osgb36_to_ll(EN) end
function convert_lonlat_nothread(EN) convert_osgb36_to_ll_nothread(EN) end
function convert_lonlat(E, N) convert_osgb36_to_ll(E, N) end

"""
Try to get OSTN15 shift parameters, and calculate offsets
"""
function get_ostn_ref(x::T, y::T) where {T<:Integer}
    # Check inputs sensible
    if (0 <= x <= 700) && (0 <= y <= 1250)
        # Calculate unique index key
        index = x + (y * 701) + 1
        # Some or None, so convert to Result, which we can try!
        se = OSTN15_DATA[][index,1]
        sn = OSTN15_DATA[][index,2]
        sg = OSTN15_DATA[][index,3]
    else
        se = NaN
        sn = NaN
        sg = NaN
    end
    return se, sn, sg
end

"""
Input values must be valid ETRS89 grid references
Calculate OSTN15 shifts for a given coordinate
"""
function ostn15_shifts(x::T, y::T) where {T<:Real}
    # Ensure that we're within the boundaries
    if check(x, (0.0, MAX_EASTING)) === NaN || check(y, (0.0, MAX_NORTHING)) === NaN
        return NaN, NaN, NaN
    end
    # Calculate base ostn15 lookup indices (type must be integer)
    e_index = trunc(Int, x / 1000.0)
    n_index = trunc(Int, y / 1000.0)

    # eastings and northings of the south-west corner of the cell
    x0 = e_index * 1000
    y0 = n_index * 1000

    # The easting, northing and geoid shifts for the four corners of the cell
    # intersections
    # this is a 3 x 4 matrix (using column-major order)
    # bottom-left grid intersection
    s0_1, s0_2, s0_3 = get_ostn_ref(e_index, n_index)
    # bottom-right
    s1_1, s1_2, s1_3 = get_ostn_ref(e_index + 1, n_index)
    # top-left
    s2_1, s2_2, s2_3 = get_ostn_ref(e_index, n_index + 1)
    # top-right
    s3_1, s3_2, s3_3 = get_ostn_ref(e_index + 1, n_index + 1)

    # offset within square
    dx = x - x0
    dy = y - y0

    t = dx / 1000.0
    u = dy / 1000.0

    # Calculation of the weights for each intersection (W)
    # this is a 4 x 1 matrix
    f0 = (1.0 - t) * (1.0 - u)
    f1 = t * (1.0 - u)
    f2 = (1.0 - t) * u
    f3 = t * u

    # Bilinear interpolation, to obtain the actual shifts
    # We could also do a dot product: weights.dot(offsets)
    se = f0 * s0_1 + f1 * s1_1 + f2 * s2_1 + f3 * s3_1
    sn = f0 * s0_2 + f1 * s1_2 + f2 * s2_2 + f3 * s3_2
    # this isn't needed for this library, since it's a height offset
    sg = f0 * s0_3 + f1 * s1_3 + f2 * s2_3 + f3 * s3_3

    return se, sn, sg
end

end # module


text:
# Generate a time series of stock prices by repeatedly solving the SPDE
# from Mastromatteo et al. (2014) to model the latent order book.
# using the previously generated mid price as the new initial mid price and
# boundary mid point during each iteration to generate a sequence of T prices,
# where T is the number of measured prices stored in the loaded
# ObjectiveFunction object.

mutable struct ReactionDiffusionPricePaths
    num_paths::Int64
    T::Int64
    p₀::Float64
    M::Int64
    β::Float64
    L::Float64
    D::Float64
    σ::Float64
    nu::Float64
    source_term::SourceTerm
    x::Array{Float64, 1}
    Δx::Float64
    Δt::Float64
end
function ReactionDiffusionPricePaths(num_paths, T, p₀, M, β,
    L, D, σ, nu, source_term)
    x₀ = p₀ - 0.5*L
    xₘ = p₀ + 0.5*L
    @assert x₀ >= 0
    # if D < 1.5 * σ
    #     @warn "D=$D is less than 1.5 times σ=$σ which runs the risk of unstable LOB ghost points"
    # end
    # Δx = L/M
    # if (β*Δx/D < log(1/3)/2) | (β*Δx/D > log(3)/2)
    #     @warn "The choice of β, D, L, M might result in jump probabilities > 0.75, thus causing large shifts in the LOB density"
    # end

    x = collect(Float64, range(x₀, xₘ, length=M+1))
    Δx = L/M
    Δt = (Δx^2) / (2.0*D)
    return ReactionDiffusionPricePaths(num_paths, T, p₀, M, β,
        L, D, σ, nu, source_term, x, Δx, Δt)
end

function (rdpp::ReactionDiffusionPricePaths)(seed::Int=-1)
    if seed == -1
            seeds = Int.(rand(MersenneTwister(), UInt32, rdpp.num_paths))
        else
            seeds = Int.(rand(MersenneTwister(seed), UInt32, rdpp.num_paths))
    end

    Δx = rdpp.L/rdpp.M
    Δt = (Δx^2) / (2.0*rdpp.D)
    time_steps = floor(Int ,rdpp.T / Δt)

    raw_price_paths = ones(Float64, time_steps + 1, rdpp.num_paths)
    raw_price_paths[1, :] .= rdpp.p₀
    sample_price_paths = ones(Float64, rdpp.T + 1, rdpp.num_paths)
    sample_price_paths[1, :] .= rdpp.p₀
    lob_densities = zeros(Float64, rdpp.M+1, time_steps + 1, rdpp.num_paths)
    P⁺s = ones(Float64, time_steps, rdpp.num_paths)
    P⁻s = ones(Float64, time_steps, rdpp.num_paths)

    for path in 1:rdpp.num_paths
         Random.seed!(seeds[path])
        lob_densities[:, :, path], raw_price_paths[:, path],
            sample_price_paths[:, path], P⁺s[:, path],
            P⁻s[:, path]  = dtrw_solver(rdpp)
    end

    return lob_densities, raw_price_paths, sample_price_paths, P⁺s, P⁻s
end

ReactionDiffusionPricePaths(dict)=ReactionDiffusionPricePaths(
    dict["num_paths"], dict["T"], dict["p₀"], dict["M"],
    dict["β"], dict["L"], dict["D"],
    dict["σ"], dict["nu"], SourceTerm(dict["λ"], dict["μ"]))

ReactionDiffusionPricePaths(;num_paths=1, T::Int64=100,
    p₀::Real=100.0, M::Int64=100, β::Real=1.0,
    L::Real=50.0, D::Real=4.0, σ::Real=1.0,
    nu::Real=0.0, λ::Real=1.0, μ::Real=0.5) =
    ReactionDiffusionPricePaths(num_paths, T, p₀,
    M, β, L, D, σ, nu, SourceTerm(λ, μ))


text:
"""
    finite-difference control velocity objective
"""
struct ControlVelocityObjective <: Objective
    R
end

function control_velocity_objective(R)
    ControlVelocityObjective(R)
end

function objective(Z, obj::ControlVelocityObjective, idx, T)
    J = 0.0

    for t = 1:T-2
        u⁻ = view(Z, idx.u[t])
        u⁺ = view(Z, idx.u[t+1])
        cv = (u⁺ - u⁻)
        J += cv' * obj.R * cv
    end

    return J
end

function objective_gradient!(∇J, Z, obj::ControlVelocityObjective, idx, T)
    J = 0.0
    for t = 1:T-2
        u⁻ = view(Z, idx.u[t])
        u⁺ = view(Z, idx.u[t+1])
        cv = (u⁺ - u⁻)
        # J += cv' * obj.R * cv

        dJdv = 2.0 * obj.R * cv
        ∇J[idx.u[t]] -= dJdv
        ∇J[idx.u[t + 1]] += dJdv
    end

    return nothing
end


text:
import StatsBase:sample

mutable struct TimeBasedSampleModel <: AbstractSampleModel
    experiences::Dict{Any, Dict{Any, NamedTuple{(:reward, :isdone, :nextstate), Tuple{Float64, Bool, Any}}}}
    nactions::Int
    κ::Float64
    t::Int
    last_visit::Dict{Tuple{Any, Any}, Int}
    TimeBasedSampleModel(nactions::Int, κ::Float64=1e-4) = new(Dict{Any, Dict{Any, NamedTuple{(:reward, :isdone, :nextstate), Tuple{Float64, Bool, Any}}}}(), nactions, κ, 0, Dict{Tuple{Any, Any}, Int}())
end

update!(m::TimeBasedSampleModel, buffer, learner) = update!(m, buffer)

function update!(m::TimeBasedSampleModel, buffer::EpisodeSARDBuffer)
   s, a, r, d, s′ = buffer.state[end-1], buffer.action[end-1], buffer.reward[end], buffer.isdone[end], buffer.state[end]
   if haskey(m.experiences, s)
         m.experiences[s][a] = (reward=r, isdone=d, nextstate=s′)
   else
         m.experiences[s] = Dict{Any, NamedTuple{(:reward, :isdone, :nextstate), Tuple{Float64, Bool, Any}}}(a => (reward=r, isdone=d, nextstate=s′))
   end
   m.t += 1
   m.last_visit[(s,a)] = m.t
end

function sample(m::TimeBasedSampleModel)
   s = rand(keys(m.experiences))
   a = rand(1:m.nactions)
   r, d, s′ = get(m.experiences[s], a, (0., false, s))
   r += m.κ * sqrt(m.t - get(m.last_visit, (s, a), 0))
   s, a, r, d, s′
end


text:
"""
# Cách dùng
- `m = MarkovChain(P)`: Tạo xích Markov với ma trận chuyển `P`
- `m = MarkovProcess(G)`: Tạo quá trình Markov với ma trận sinh `G`
- `stationarydist(m)`: Phân phối dừng của `m`, với `m` là xích hoặc quá trình Markov
- `cost(m, ...)`: Chi phí lâu dài của quá trình Markov `m`, `?cost` để xem thêm.
- `imbeddedchain(m)`: Xích Markov nhúng của quá trình Markov `m`.
- `transition_matrix(m, t)`: Tạo ma trận chuyển `P(t)` tại thời điểm `t`
"""
module MarkovStuff
using LinearAlgebra

export MarkovProcess, MarkovChain, stationarydist, cost,
	imbeddedchain, steadydist, transition_matrix

"""
	MarkovChain(G::AbstractMatrix)

Xích Markov với ma trận chuyển P
"""
struct MarkovChain
	P::AbstractMatrix
	function MarkovChain(P::AbstractMatrix{Real})
		if any(1 .≉ round.(sum(P; dims=2); digits=12))
			@error "Tổng hàng khác 1"
		else
			new(P)
		end
	end
end

"""
	MarkovProcess(G::AbstractMatrix)

Quá trình Markov với ma trận sinh G
"""
struct MarkovProcess
	G::AbstractMatrix

	function MarkovProcess(G::AbstractMatrix{Real})
		if any(0 .≉ round.(sum(G; dims=2); digits=12))
			@error "Tổng hàng khác 0"
		else
			new(G)
		end
	end

	function MarkovProcess(G::AbstractMatrix)
		new(G)
	end

	function MarkovProcess(m::MarkovChain, λ::AbstractVector)
		G = copy(m.P)
		G = G .* λ
		G[diagind(G)] .= .-λ
		return new(G)
	end
end

include(joinpath(@__DIR__, "markovchain.jl"))
include(joinpath(@__DIR__, "markovprocess.jl"))

end


text:
# This file contains functions to create the actions dataframe.

"""
    function convert_to_actions(event_df)
Convert Wyscout events to SPADL actions.
Step 1: initialize a vector of WyscoutData
Step 2: create new positions (start_x, start_y and end_x, end_y)
"""
function convert_to_actions(event_df::Vector{WyscoutEvent})::Vector{RegularSPADL}

    vector_wyscout_data = create_wyscout_data_processing(event_df)
    vector_wyscout_data = make_new_positions(vector_wyscout_data)
    vector_wyscout_data = get_tags(vector_wyscout_data)
    vector_wyscout_data = fix_events(vector_wyscout_data)
    vector_spadl = create_actions(vector_wyscout_data)

    return vector_spadl
end

"""
    function create_wyscout_data_processing(event_df)
"""
function create_wyscout_data_processing(event_df::Vector{WyscoutEvent})::Vector{WyscoutData}

    vector_wyscout_data = Vector{WyscoutData}()

    for i in eachindex(event_df)
        tmp_event_fixed = WyscoutEventFixed(game_id = event_df[i].game_id,
                                event_id = event_df[i].event_id,
                                period_id = event_df[i].period_id,
                                team_id = event_df[i].team_id,
                                type_id = event_df[i].type_id,
                                player_id = event_df[i].player_id,
                                milliseconds = event_df[i].milliseconds)

        tmp = WyscoutData(event = event_df[i], event_fixed = tmp_event_fixed)
        push!(vector_wyscout_data, tmp)
    end

    return vector_wyscout_data
end

"""
    function create_actions(event_df)
Create the SciSports action format.
"""
function create_actions(vector_wyscout_data::Vector{WyscoutData})::Vector{RegularSPADL}

    configs = ActionsConfig()

    vector_spadl = Vector{RegularSPADL}()

    for i in eachindex(vector_wyscout_data)
        spadl_element = RegularSPADL(game_id = vector_wyscout_data[i].event_fixed.game_id,
                                original_event_id = vector_wyscout_data[i].event_fixed.event_id,
                                period_id = vector_wyscout_data[i].event_fixed.period_id,
                                team_id = vector_wyscout_data[i].event_fixed.team_id,
                                start_x = vector_wyscout_data[i].event_fixed.start_x,
                                start_y = vector_wyscout_data[i].event_fixed.start_y,
                                end_x = vector_wyscout_data[i].event_fixed.end_x,
                                end_y = vector_wyscout_data[i].event_fixed.end_y,
                                time_seconds = vector_wyscout_data[i].event_fixed.milliseconds # to be checked
                                )

        # determine bodypart_id
        spadl_element = determine_bodypart_id(spadl_element, vector_wyscout_data[i], configs)

        # determine type_id
        spadl_element = determine_type_id(spadl_element, vector_wyscout_data[i], configs)

        # determine result_id
        spadl_element = determine_result_id(spadl_element, vector_wyscout_data[i], configs)

        push!(vector_spadl, spadl_element)
    end

    # remove non actions
    vector_spadl = filter(x -> x.type_name != "non_action", vector_spadl)

    return vector_spadl
end

"""
    function determine_bodypart_id()
Determine the body part for each action.
"""
function determine_bodypart_id(spadl_element::RegularSPADL, wyscout_data::WyscoutData, config::ActionsConfig)::RegularSPADL

    if wyscout_data.event_fixed.subtype_id in [81, 36, 21, 90, 91]
        spadl_element.bodypart_name = "other"
    elseif wyscout_data.event_fixed.subtype_id == 82
        spadl_element.bodypart_name = "head"
    elseif  wyscout_data.event_fixed.subtype_id == 10 && wyscout_data.tags.head_body
        spadl_element.bodypart_name = "head/other"
    else
        spadl_element.bodypart_name = "foot"
    end
    spadl_element.bodypart_id = findall(x -> x == spadl_element.bodypart_name, config.bodyparts)[1] - 1
    return spadl_element
end

"""
    function determine_type_id()

This function transforms the Wyscout events, sub_events and tags
into the corresponding SciSports action type
"""
function determine_type_id(spadl_element::RegularSPADL, wyscout_data::WyscoutData, config::ActionsConfig)::RegularSPADL

    if wyscout_data.tags.own_goal
        spadl_element.type_name = "bad_touch"
    elseif wyscout_data.event_fixed.type_id == 8
        if wyscout_data.event_fixed.subtype_id == 80
            spadl_element.type_name = "cross"
        else
            spadl_element.type_name = "pass"
        end
    elseif wyscout_data.event_fixed.subtype_id == 36
        spadl_element.type_name = "throw_in"
    elseif wyscout_data.event_fixed.subtype_id == 30
        if wyscout_data.tags.high
            spadl_element.type_name = "corner_crossed"
        else
            spadl_element.type_name = "corner_short"
        end
    elseif  wyscout_data.event_fixed.subtype_id == 32
        spadl_element.type_name = "freekick_crossed"
    elseif wyscout_data.event_fixed.subtype_id == 31
        spadl_element.type_name = "freekick_short"
    elseif wyscout_data.event_fixed.subtype_id == 34
        spadl_element.type_name = "goalkick"
    elseif wyscout_data.event_fixed.type_id == 2
        spadl_element.type_name = "foul"
    elseif wyscout_data.event_fixed.type_id == 10
        spadl_element.type_name = "shot"
    elseif wyscout_data.event_fixed.subtype_id == 35
        spadl_element.type_name = "short_penalty"
    elseif wyscout_data.event_fixed.subtype_id == 33
        spadl_element.type_name = "shot_freekick"
    elseif wyscout_data.event_fixed.type_id == 9
        spadl_element.type_name = "keeper_save"
    elseif wyscout_data.event_fixed.subtype_id == 71
        spadl_element.type_name = "clearance"
    elseif wyscout_data.event_fixed.subtype_id == 72 && wyscout_data.tags.not_accurate
        spadl_element.type_name = "bad_touch"
    elseif wyscout_data.event_fixed.subtype_id == 70
        spadl_element.type_name = "dribble"
    elseif wyscout_data.tags.take_on_left || wyscout_data.tags.take_on_right
        spadl_element.type_name = "take_on"
    elseif wyscout_data.tags.sliding_tackle
        spadl_element.type_name = "tackle"
    elseif wyscout_data.tags.interception && (wyscout_data.event_fixed.subtype_id in [0, 10, 11, 12, 13, 72])
        spadl_element.type_name = "interception"
    else
        spadl_element.type_name = "non_action"
    end
    spadl_element.type_id = findall(x -> x == spadl_element.type_name, config.actionstypes)[1] - 1

    return spadl_element
end

"""
    function determine_result_id()
Determine the result of each event.
"""
function determine_result_id(spadl_element::RegularSPADL, wyscout_data::WyscoutData, config::ActionsConfig)::RegularSPADL

    if wyscout_data.tags.offside spadl_element.result_id = 2 end
    if wyscout_data.event_fixed.type_id == 2 spadl_element.result_id = 1 end
    if wyscout_data.tags.goal spadl_element.result_id = 1 end
    if wyscout_data.tags.own_goal spadl_element.result_id = 3 end
    if wyscout_data.event_fixed.subtype_id in [100,33,35] spadl_element.result_id = 0 end
    if wyscout_data.tags.accurate spadl_element.result_id = 1 end
    if wyscout_data.tags.not_accurate spadl_element.result_id = 0 end
    if wyscout_data.tags.interception || wyscout_data.tags.clearance || wyscout_data.event_fixed.subtype_id == 71
        spadl_element.result_id = 1
    end
    if wyscout_data.event_fixed.type_id == 9 spadl_element.result_id = 1 end

    # equivalent of the line 655 of the initial python code
    if isnothing(spadl_element.result_id) spadl_element.result_id = 1 end

    spadl_element.result_name = config.actionstypes[spadl_element.result_id+1]

    return spadl_element
end


text:
# Copyright (c) 2021 California Institute of Technology (“Caltech”).
# U.S. Government sponsorship acknowledged.
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# - Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# - Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# - Neither the name of Caltech nor its operating division, the Jet Propulsion
#   Laboratory, nor the names of its contributors may be used to endorse or
#   promote products derived from this software without specific prior written
#   permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
# THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
# PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
# OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
# OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
# ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

struct FormatError <: Exception end

Base.showerror(io::IO, e::FormatError) = print(io, "Error parsing file.")

const formatversion_angfab = 1
const formatversion_angradfab = 3

######## I/O AngularCryoFaB
function Base.write(io::IO, cfb::AngularCryoFaB)
    n = 0
    # assumes RingOrder
    n += write(io, Int64(formatversion_angfab))
    n += write(io, Int64(cfb.nside))
    n += write(io, Array{Float64,1}(cfb.mask))
    n += write(io, Int64(cfb.npixsurvey))
    n += write(io, Array{Int64,1}(cfb.hpix_full2survey))
    n += write(io, Vector{Float64}(cfb.ell))
    n += write(io, Matrix{Float64}(cfb.Trans))
    n += write(io, Matrix{Float64}(cfb.TransInv))
    #@show size(cfb.ell) size(cfb.Trans) size(cfb.TransInv)
    return n
end

function Base.read(io::IO, ::Type{AngularCryoFaB})
    # assumes RingOrder
    diskversion = read(io, Int64)
    if diskversion != formatversion_angfab
        @error "Format version mismatch." diskversion formatversion_angfab io
        throw(FormatError())
    end
    nside = read(io, Int64)
    npix = nside2npix(nside)
    mask = HealpixMap{Float64,Healpix.RingOrder}(read!(io, Array{Float64,1}(undef, npix)))
    npixsurvey = read(io, Int64)
    hpix_full2survey = read!(io, Array{Int64,1}(undef, npix))
    ell = read!(io, Vector{Float64}(undef, npixsurvey))
    Y = read!(io, Matrix{Float64}(undef, npixsurvey, npixsurvey))
    Yinv = read!(io, Matrix{Float64}(undef, npixsurvey, npixsurvey))
    return AngularCryoFaB(mask, nside, npixsurvey, hpix_full2survey, ell, Y, Yinv)
end

function AngularCryoFaB(fname::AbstractString)
    io = open(fname, read=true)
    cfb = read(io, AngularCryoFaB)
    close(io)
    return cfb
end

function AngRadCryoFaB(fname::AbstractString)
    io = open(fname, read=true)
    cfb = read(io, AngRadCryoFaB)
    close(io)
    return cfb
end

######## I/O BlockDiagonals
function Base.write(io::IO, bd::BlockDiagonal{T}) where {T}
    n = 0
    allblocks = blocks(bd)
    numblocks = length(allblocks)
    n += write(io, Int64(numblocks))
    for i=1:numblocks
        block = allblocks[i]
        len1, len2 = size(block)
        n += write(io, Int64(len1), Int64(len2))
        n += write(io, Array{T,2}(block))
    end
    return n
end

function Base.read(io::IO, ::Type{BlockDiagonal{T}}) where {T}
    numblocks = read(io, Int64)
    allblocks = Matrix{T}[]
    for i=1:numblocks
        len1 = read(io, Int64)
        len2 = read(io, Int64)
        block = read!(io, Array{T,2}(undef, len1, len2))
        push!(allblocks, block)
    end
    return BlockDiagonal(allblocks)
end

######## I/O AngRadCryoFaB
function Base.write(io::IO, cfb::AngRadCryoFaB)
    n = 0
    n += write(io, Int64(formatversion_angradfab))

    n += write(io, cfb.angfab)
    n += write(io, Int64(length(cfb.Angl)))
    n += write(io, Array{Float64,1}(cfb.Angl))
    n += write(io, Int64(length(cfb.reff)))  # save 'nbins'
    #n += write(io, cfb.AngTrans)    # not needed, since easily reconstructed
    #n += write(io, cfb.AngTransInv) # not needed, since easily reconstructed

    n += write(io, Int64(length(cfb.sortbyl)))  # save length
    n += write(io, Array{Int64,1}(cfb.sortbyl))

    n += write(io, Float64(cfb.rmin))
    n += write(io, Float64(cfb.rmax))
    n += write(io, Int64(length(cfb.rbounds)))  # save length
    n += write(io, Array{Float64,1}(cfb.rbounds))
    n += write(io, Int64(length(cfb.reff)))  # save length
    n += write(io, Array{Float64,1}(cfb.reff))

    n += write(io, Int64(length(cfb.RadK)))  # save length
    n += write(io, Array{Float64,1}(cfb.RadK))
    n += write(io, cfb.RadTrans)
    n += write(io, cfb.RadTransInv)
    return n
end

function Base.read(io::IO, ::Type{AngRadCryoFaB})
    # assumes RingOrder
    diskversion = read(io, Int64)
    if diskversion != formatversion_angradfab
        @error "Format version mismatch." diskversion formatversion_angradfab io
        throw(FormatError())
    end

    angfab = read(io, AngularCryoFaB)
    Angl_len = read(io, Int64)
    Angl = read!(io, Vector{Float64}(undef, Angl_len))
    nbins = read(io, Int64)
    AngTrans = BlockDiagonal(fill(angfab.Trans, nbins))
    AngTransInv = BlockDiagonal(fill(angfab.TransInv, nbins))

    len = read(io, Int64)
    sortbyl = read!(io, Array{Int64,1}(undef, len))

    rmin = read(io, Float64)
    rmax = read(io, Float64)
    len = read(io, Int64)
    rbounds = read!(io, Array{Float64,1}(undef, len))
    len = read(io, Int64)
    reff = read!(io, Array{Float64,1}(undef, len))

    len = read(io, Int64)
    RadK = read!(io, Array{Float64,1}(undef, len))
    RadTrans = read(io, BlockDiagonal{Float64})
    RadTransInv = read(io, BlockDiagonal{Float64})
    return AngRadCryoFaB(angfab, Angl, AngTrans, AngTransInv,
        sortbyl, rmin, rmax, rbounds, reff,
        RadK, RadTrans, RadTransInv)
end

# vim: set sw=4 et sts=4 :


text:
"""
    fetch_dataset(dataset_id::AbstractString; restrictions::Dict{<:AbstractString, <:AbstractArray{<:AbstractString, 1}}, vintage::AbstractString = "latest", dimensions::Symbol = :code)

Fetch dataset with dataset_id `dataset_id`.

Restrictions is an (optional) dictionary with keys (corresponding to names of dimensions)
mapping to an Array of codes. Only series that have those codes for that dimension are
returned in the dataset.

Vintage can be any of `current` (alias `latest`, default), `previous` or a date in YYYY-MM-DD
form.

`dimensions` can be one of :code or :name. This dictates whether the values for each dimension
of the dataset are returned as the codes for each dimension, or the long names that correspond
to those codes.

Returns a dataframe, including columns: `dates`, `vintage`, `series_id`, `values`. The dataframe
will also have columns for each of the dimensions of the dataset.
"""
function fetch_dataset(dataset_id::AbstractString; restrictions::Dict{<:AbstractString, <:AbstractArray{<:AbstractString, 1}} = Dict{String, Array{String, 1}}(), vintage::AbstractString = "latest", dimensions::Symbol = :code)
    check_valid_vintage(vintage)
    check_valid_code(dataset_id)

    if dimensions != :code && dimensions != :name
        throw(ErrorException("`dimensions` can only be one of `:code` or `:name`"))
    end

    restrictions_string = if length(restrictions) > 0
        """"restrictions" : $(JSON.json(restrictions)),"""
    else
        ""
    end

    # Get the data
    res = HTTP.request(
        "POST",
        "https://aeconomica.io/api/v1/dataset",
        [("Content-Type", "application/json")],
        """{
            "dataset" : "$dataset_id",
            $(restrictions_string)
            "vintage" : "$vintage",
        "apikey" : "$(apikey())"}""",
        status_exception = false)
    result = if res.status == 200
        response = JSON.parse(String(res.body), null = missing)
        df = reduce(vcat, map(x -> DataFrames.DataFrame(x), response))
        df.dates = Dates.Date.(df.dates)
        df.vintage = Dates.Date.(df.vintage)
        df.values = if any(ismissing.(df.values))
            Vector{Union{Float64, Missing}}(df.values)
        else
            Vector{Float64}(df.values)
        end

        df
    else
        error = JSON.parse(String(res.body))["error"]
        if error[1:21] == "500 Internal Error - "
            error = error[22:end]
            throw(ErrorException(error))
        else
            throw(ErrorException(error))
        end
    end

    # Convert to names if desired
    if dimensions == :name
        # Get the dataset structure
        res = HTTP.request(
            "POST",
            "https://aeconomica.io/api/v1/dataset_structure",
            [("Content-Type", "application/json")],
            """{
                "dataset" : "$dataset_id",
            "apikey" : "$(apikey())"}""",
            status_exception = false)
        structure = if res.status == 200
            JSON.parse(String(res.body), null = missing)
        else
            error = JSON.parse(String(res.body))["error"]
            if error[1:21] == "500 Internal Error - "
                error = error[22:end]
                throw(ErrorException(error))
            else
                throw(ErrorException(error))
            end
        end

        for dim in structure["dimensions"]
            #Iterate over each dimension, and replace the codes in the column of the df with their name
            result[!, dim["dimname"]] = map(code -> dim["options"][findfirst(map(x -> code == x["code"], dim["options"]))]["name"], result[:, dim["dimname"]])
        end
    end

    result
end


text:
#@testset "BigMat.jl" begin
    A        = randn(20,5)
    Attm,err = MPT_SVD(A,[5 2 2; 5 1 1],0.0)

    # test qr
    @run Qttm,R = qr(Attm,0.0)
    Q = mpo2mat(Qttm)
    norm(qr(A).R-R)
    norm(Array(qr(A).Q)-mpo2mat(Qttm))

    Q = randTTmsubspaceiter(q,A,O,acc)
    W,S,V = svd(Attm)

    @run U,S,V    = TNrSVD(Attm,30,10,10,0.0)

#end

max_stars_repo_name: clarazen/BigMat

max_issues_repo_name: clarazen/BigMat

max_forks_repo_name: clarazen/BigMat


text:
macro fill_series_1(expr)
  left, right = expr.args
  i = left.args[2]
  quote
    for $i in 1:(length($(left.args[1])))
      try $expr
        catch BoundsError
      end
    end
  end
end

function index_range(expr::AbstractString)
  min = 0
  max = 0

  s = parse(expr)
  symbol = s.args[1].args[1]

  for i in s.args[2].args
    if(typeof(i) == Expr)
      tmp = i.args[2]

      if(tmp.args[1] == symbol)
        tmp = tmp.args[2]
      end

      s = string(tmp)
      s = s[3:length(s)]
      val = parse(Int64,s,10)

      if(val < min)
        min = val
      elseif(val > max)
        max = val
      end

    end
  end

  min = abs(min) + 1
  max = N - max
  return [min, max]
end

macro fill_series_2(expr)
  minval = index_range(string(expr))[1]
  maxval = index_range(string(expr))[2]
  quote
    t = $minval
    while t <= $maxval
      $expr
      t += 1
    end
  end
end

max_stars_repo_name: wchmielarz/NobodyDoesItLikeJuliet

max_issues_repo_name: wchmielarz/NobodyDoesItLikeJuliet

max_forks_repo_name: wchmielarz/NobodyDoesItLikeJuliet


text:
using Test, Turf

@testset "Turf" begin
    include("angle.jl")
    include("bearing.jl")
    include("utils.jl")
    include("distance.jl")
    include("destination.jl")
    include("booleans.jl")
    include("centering.jl")
    include("transformations.jl")
    include("ellipse.jl")
    include("circle.jl")
    include("planes.jl")
    include("grids.jl")
    include("lines.jl")
    include("bbox.jl")
    #include("bezier.jl")
    include("square.jl")
end # begin


text:
using StaticArrays, HybridArrays, Test, LinearAlgebra

@testset "AbstractArray interface" begin
    @testset "one()" begin
        M = HybridMatrix{2, StaticArrays.Dynamic()}([1 2; 4 5])
        @test (@inferred one(M)) == @SMatrix [1 0; 0 1]
        @test isa(one(M), HybridMatrix{2,StaticArrays.Dynamic(),Int})

        Mv = view(M, :, SOneTo(2))
        @test (@inferred one(Mv)) === @SMatrix [1 0; 0 1]
    end

    @testset "zero()" begin
        M = HybridMatrix{2, StaticArrays.Dynamic()}([1 2; 4 5])
        @test (@inferred zero(M)) == @SMatrix [0 0; 0 0]
        @test isa(zero(M), HybridMatrix{2,StaticArrays.Dynamic(),Int})

        Mv = view(M, :, SOneTo(2))
        @test (@inferred zero(Mv)) === @SMatrix [0 0; 0 0]

    end

    @testset "fill!()" begin
        M = HybridMatrix{2, StaticArrays.Dynamic(), Float64}([1 2; 4 5])
        fill!(M, 3)
        @test all(M .== 3.)
    end

end