text:
Figure1.jl

######## Checking to see if the values in Figure 1, of rates varying with fertilisation, align with our calcs

#### all species 0.01g
Z = 1
M = BodyMasses(Z=Z)

#### at 20C
T = 293.15

# Check limits
fert = [1:1:20;]
allI_K = [1:0.1:20;]

######### 1a - K
K1 = map(x -> (carryingcapacity_BA(x, param, M, T)[1]), fert)
pltK1 = plot(fert, K1, c = :orange , lw = 3
    , ylabel = "K", xlabel = "fertilisation"
    , label = "1a", leg = :topleft
)

# carryingcapacity_BA(I_Kmin, param, M, T) # 0.2754228703338166
# carryingcapacity_BA(5, param, M, T)      # 1.3771143516690831
# carryingcapacity_BA(10, param, M, T)     # 2.7542287033381663
# carryingcapacity_BA(I_Kmax, param, M, T) # 5.5084574066763325
# # CORRECT

######### 1b - r
r1 = fill(growth_BA(param, M, T)[1], 20)
pltr1 = plot(fert, r1/1e-7, c = :orange , lw = 3, ylims = (3,7)
    , ylabel = "r", xlabel = "fertilisation"
    , label = "1b", leg = :topleft
)

# growth_BA(param, M, T) # 4.900749719809166e-7
# # CORRECT

######### 1c - x relative to r

xrel1 = fill(metabolism_BA(param, M, T)[2], 20) ./ r1
pltxrel1 = plot(fert, xrel1, c = :orange , lw = 3, ylims = (0.3,0.8)
    , ylabel = "x_rel", xlabel = "fertilisation"
    , label = "1c", leg = :topleft
)

# metabolism_BA(param, M, T) # 2.7338159268322533e-7
# 2.7338159268322533e-7/4.900749719809166e-7 # 0.5578362665169335
# # CORRECT

######### 1d - y relative to x

yrel1 = fill(max_ingestion_BA(param, M, Z, T)[2]/metabolism_BA(param, M, T)[2], 20)
pltyrel1 = plot(fert, yrel1, c = :orange , lw = 3, ylims = (25,55)
    , ylabel = "y_rel", xlabel = "fertilisation"
    , label = "1d", leg = :topleft
)

# max_ingestion_BA(param, M, Z, T) # 1.11074e-5
# 1.11074e-5/2.7338159268322533e-7 # 40.62965575326955
# # CORRECT

######### 1e - B0 relative to K
B0rel1 = fill(half_saturation_BA(param, M, Z , T)[2], 20) ./ K1
pltB0rel1 = plot(fert, B0rel1, c = :orange , lw = 3, ylims = (0,10)
    , ylabel = "B0_rel", xlabel = "fertilisation"
    , label = "1e", leg = :topleft
)

half_saturation_BA(param, M, Z , T) # 2.63599
2.63599/0.2754228703338166 # 9.570701215934395 # I_K = 1
2.63599/1.3771143516690831 # 1.914140243186879 # I_K = 5
2.63599/2.7542287033381663 # 0.9570701215934395 # I_K = 10
2.63599/5.5084574066763325 # 0.47853506079671976 # I_K = 20
# CORRECT

l = @layout [a ; b c; d e]
plot(pltK1, pltr1, pltxrel1, pltyrel1, pltB0rel1,layout = l)

######### 1f - biomass extremes
# fertilisation gradient
allI_K = [1:1:20;]
# time span for each simulation
tspan = (0.0, 315360000000.0)
# define extinction threshold
threshold = 10e-12
#### data frame to store outputs
df = DataFrame(
    Z = [],
    I_K = [],
    T = [],
    K = [],
    r = [],
    x_2 = [],
    y_21 = [],
    B0_21 = [],
    maxB1 = [],
    maxB2 = [],
    maxB3 = [],
    minB1 = [],
    minB2 = [],
    minB3 = [],
    sp_survival = []
)

### Loop over fertilisation gradient

for I_K in allI_K
    # generate parameters specific to this combo of Z, I_K and T
    p = ModelParameters(param, T, I_K, Z)

    # Calculate starting biomasses
    u0 = [p.K[1]/2, p.K[1]/2, p.K[1]/2]

    # Define the problem
    prob = ODEProblem(BEFW, u0, tspan, p)

    # Solve the problem
    sol = solve(prob, maxiters = 1e7)

    ###### Output metrics ######
    # turn vector output of biomasses into matrix
    matrix_u = hcat(sol.u...)' # 56×3 adjoint(::Matrix{Float64}) with eltype Float64:

    ### Biomass extremes
    maxB1 = maximum(matrix_u[2 * Int64(round(length(matrix_u[:,1])/3)):Int64(length(matrix_u[:,1])),1])
    maxB2 = maximum(matrix_u[2 * Int64(round(length(matrix_u[:,1])/3)):Int64(length(matrix_u[:,1])),2])
    maxB3 = maximum(matrix_u[2 * Int64(round(length(matrix_u[:,1])/3)):Int64(length(matrix_u[:,1])),3])

    minB1 = minimum(matrix_u[2 * Int64(round(length(matrix_u[:,1])/3)):Int64(length(matrix_u[:,1])),1])
    minB2 = minimum(matrix_u[2 * Int64(round(length(matrix_u[:,1])/3)):Int64(length(matrix_u[:,1])),2])
    minB3 = minimum(matrix_u[2 * Int64(round(length(matrix_u[:,1])/3)):Int64(length(matrix_u[:,1])),3])

    ### Number of surviving species
    # final biomasses
    final_u = matrix_u[length(sol.t),:]

    # count number of species remaining that are above extinction threshold
    survived = count(final_u .> threshold)

    ### Push outputs to DataFrame
    push!(df, [
        Z,
        I_K,
        T,
        p.K[1],
        p.r[1],
        p.x[2],
        p.y[2,1],
        p.B0[2,1],
        maxB1,
        maxB2,
        maxB3,
        minB1,
        minB2,
        minB3,
        survived
    ])

    # print some stuff - see how the simulation is progressing
    println(("Z = $Z", "I_K = $I_K", "T = $T"))
end

# plot
# initialise an empty plot
pltextr1 = Plots.plot([NaN], [NaN],
                label = "",
                xlims = (0,20),
                ylims = (0,5.5),
                foreground_colour_legend = nothing,
                xlabel = "fertilisation",
                ylabel = "biomass extremes",
                leg = :topleft)

plot!(pltextr1, allI_K, df.maxB1, ls = :dot, c = :green, lw = 3, label = "basal")
plot!(pltextr1, allI_K, df.minB1, ls = :dot, c = :green, lw = 3, label = "")
plot!(pltextr1, allI_K, df.maxB2, ls = :dot, c = :blue, lw = 3, label = "intermediate")
plot!(pltextr1, allI_K, df.minB2, ls = :dot, c = :blue, lw = 3, label = "")
plot!(pltextr1, allI_K, df.maxB3, ls = :dot, c = :red, lw = 3, label = "top")
plot!(pltextr1, allI_K, df.minB3, ls = :dot, c = :red, lw = 3, label = "")

l = @layout [a ; b c; d e; f]
plot(pltK1, pltr1, pltxrel1, pltyrel1, pltB0rel1, pltextr1, layout = l)


text:
using MosekTools
using Test
using JuMP

include("jump_sdp.jl")
include("jump_lp.jl")
include("jump_soc.jl")

@testset "JuMP tests" begin
    o = Mosek.Optimizer()
    MOI.set(o, MOI.RawOptimizerAttribute("LOG"), 0)
    test_jump_lp()
    MOI.empty!(o)
    test_jump_soc(o)
    MOI.empty!(o)
    test_jump_sdp(o)
end

ext: jl


text:
include(joinpath(@__DIR__, "..", "src/dynamics", "hopper_3D", "visuals.jl"))
T = Float64
vis = Visualizer()
open(vis)

# qet hopper model
s = get_simulation("hopper_3D", "flat_3D_lc", "flat")
model = s.model
env = s.env

nq = model.dim.q
nu = model.dim.u
nc = model.dim.c
nb = model.dim.c * friction_dim(env)
nd = nq + nc + nb
nr = nq + nu + nc + nb + nd

H = 92
# h = 0.04
# h = 0.02
h = 0.01
κ = 1.0e-8

# Design open-loop control trajectory
q0_ref = SVector{nq,T}([0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5])
q1_ref = SVector{nq,T}([0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5])
# α = 0.02937 # N=45 h=0.04
# α = 0.01468 # N=45 h=0.02
α = 0.0077 # N=46 h=0.01
# u_ref = []
# push!(u_ref,  [[0.0,  5.0*α*model.g*(model.mb+model.ml)/2] for k=1:6]...);
# push!(u_ref,  [[0.0, -0.9*α*model.g*(model.mb+model.ml)/2] for k=1:8]...);
# push!(u_ref,  [[0.0,  0.2*α*model.g*(model.mb+model.ml)/2] for k=1:14]...);
# push!(u_ref,  [[0.0,  2.1*α*model.g*(model.mb+model.ml)/2] for k=1:H-1-length(u_ref)]...);

u_ref = []
push!(u_ref,  [[0.0, 0.0,   5.0*α*model.g*(model.mb+model.ml)/2] for k=1:2*6]...);
push!(u_ref,  [[0.0, 0.0,  -0.60*α*model.g*(model.mb+model.ml)/2] for k=1:2*10]...);
push!(u_ref,  [[0.0, 0.0,   0.14*α*model.g*(model.mb+model.ml)/2] for k=1:2*15]...);
push!(u_ref,  [[0.0, 0.0,   2.19*α*model.g*(model.mb+model.ml)/2] for k=1:H-length(u_ref)]...);

# Simulate
# sim = simulator(s, q0_ref, q1_ref, h, H;
#     p = open_loop_policy(u_ref),
#     d = no_disturbances(model),
#     r! = model.res.r!,
#     rz! = model.res.rz!,
#     rθ! = model.res.rθ!,
#     rz = model.spa.rz_sp,
#     rθ = model.spa.rθ_sp,
#     ip_opts = InteriorPointOptions(r_tol=1e-12, κ_tol=2κ),
#     sim_opts = SimulatorOptions())

sim = simulator(s, q0_ref, q1_ref, h, H,
    p = open_loop_policy(u_ref),
    d = no_disturbances(model),
    ip_opts = InteriorPointOptions(
        r_tol = 1.0e-12,
        # κ_init = 1.0e-8,
        κ_tol = 1.0e-8),
    sim_opts = SimulatorOptions(warmstart = true)
    )

simulate!(sim)

sim.traj.z
sim.traj.θ
plt = plot(legend=false)
for t = 1:80
    r = residual(model, env, sim.traj.z[t], sim.traj.θ[t], [κ])
    @show norm(r)
    nz = num_var(model, env)
    # plot!(log.(10, abs.(r[1:nz])))
    # plot!(Vector(1:nz), log.(10, abs.(r[10:14)))
    plot!(log.(10, abs.(r .+ 1e-20)))
    # plot!(Vector(1:nz), log.(10, abs.(r)))
    # plot!(log.(10, abs.(r[1:nq+10])))
end
display(plt)

sim.traj.q[40][3]
sim.traj.q[4][7]
plot([q[3] for q in sim.traj.q])
plot!([q[7] for q in sim.traj.q])

plot([q[4] for q in sim.traj.q])
plot([q[5] for q in sim.traj.q])
plot([q[6] for q in sim.traj.q])

plot(hcat(Vector.(sim.traj.u)...)')
plot(hcat(Vector.(sim.traj.q)...)'[:,4:6])
anim = visualize_robot!(vis, model, sim.traj)
anim = visualize_force!(vis, model, env, sim.traj, anim=anim)

# Check ~perfect loop
@show round.(sim.traj.q[end-2] - q0_ref, digits=3)
@show round.(sim.traj.q[end-1] - q1_ref, digits=3)

# Save trajectory
traj = deepcopy(sim.traj)
gait_path = joinpath(@__DIR__, "..", "dynamics", "hopper_3D", "gaits", "gait_in_place.jld2")
@save gait_path traj

# Reload trajectory
res = JLD2.jldopen(gait_path)
loaded_traj = res["traj"]

traj = deepcopy(get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/hopper_3D/gaits/gait_in_place.jld2"),
    load_type = :joint_traj))
plot(hcat(Vector.(traj.q)...)')


text:
import Base: eachline, start, done, next, eltype

type EachLineString
    string::ASCIIString
    α::Int
    len::Int
end

function eachline(string::ASCIIString)
  len = length(string)
  if len > 0
    EachLineString(string, 0, len)
  else
    throw(ErrorException("Empty String"))
  end
end

start(itr::EachLineString) = 0

done(itr::EachLineString, state) = itr.α >= itr.len || state >= itr.len

function next(itr::EachLineString, state)
  itr.α = state + 1
  index = searchindex(itr.string, '\n', state + 1)
  state = index == 0 ? itr.len : index
  itr.string[itr.α:state], state
end

eltype(::Type{EachLineString}) = ASCIIString


text:
"""
    Atom(T <: ThreeD, ::Int64, ::Union{Real, Integer})
"""
function Atom(
    geometry::Sphere,
    N::Int64,
    kR::Union{Real,Integer};
    createFunction = ftn_AtomsOnSphere::Function,
)
    @debug "start: Shape - Sphere"

    dimensions = 3
    ρ = 3N / (4π * kR^3)
    rₘᵢₙ = get_rₘᵢₙ(ρ)
    if rₘᵢₙ ≥ kR / 10
        rₘᵢₙ = kR / 100
    end

    r = get_atoms(dimensions, N, rₘᵢₙ; createFunction, kR)

    @debug "end  : Shape - Sphere"
    return Atom(Sphere(), r, N, Float64(kR))
end

function Atom(geometry::Sphere, r::Matrix, kR::Union{Real,Integer})
    N = size(r, 2) # remember to use each collum as a atom position
    return Atom(Sphere(), r, N, Float64(kR))
end

function ftn_AtomsOnSphere(; kwargs...)
    kR = kwargs[:kR]
    new_atom = zeros(3)
    new_atom[1] = 2π * rand() # azimuth
    new_atom[2] = asin(2 * rand() - 1) # elevation
    new_atom[3] = kR * (rand()^(1 ./ 3.0)) # radii
    new_atom[:] = sph2cart(new_atom)
    return new_atom
end


text:
@defcomp impactheating begin
    regions = Index()

    heating = Variable(index=[time,regions])

    hebm = Parameter(index=[regions])
    gdp90 = Parameter(index=[regions])

    population = Parameter(index=[time,regions])
    pop90 = Parameter(index=[regions])

    income = Parameter(index=[time,regions])
    heel = Parameter(default = 0.8)

    temp = Parameter(index=[time,regions])
    henl = Parameter(default = 0.5)

    cumaeei = Parameter(index=[time,regions])

    function run_timestep(p, v, d, t)

        if !is_first(t)
            for r in d.regions
                ypc = p.income[t, r] / p.population[t, r] * 1000.0
                ypc90 = p.gdp90[r] / p.pop90[r] * 1000.0

                v.heating[t, r] = p.hebm[r] * p.cumaeei[t, r] * p.gdp90[r] * atan(p.temp[t, r]) / atan(1.0) * (ypc / ypc90)^p.heel * p.population[t, r] / p.pop90[r]
            end
        end
    end
end


text:
"""
    get_sdf(contact, x, q)

    returns the signed distance for a contact

    contact: ContactConstraint
    x: body position
    q: body orientation
"""
function get_sdf(contact::ContactConstraint{T,N,Nc,Cs}, x::AbstractVector{T},
    q::Quaternion{T}) where {T,N,Nc,Cs<:Contact{T,N}}
    model = contact.model
    return model.collision.contact_normal * (x + vector_rotate(model.collision.contact_origin, q)) - model.collision.contact_radius
end

function get_sdf(mechanism::Mechanism{T,Nn,Ne,Nb,Ni}, storage::Storage{T,N}) where {T,Nn,Ne,Nb,Ni,N}
    d = []
    for contact in mechanism.contacts
        ibody = get_body(mechanism, contact.parent_id).id - Ne
        push!(d, [get_sdf(contact, storage.x[ibody][i], storage.q[ibody][i]) for i = 1:N])
    end
    return d
end

"""
    contact_location(contact, x, q)

    location of contact point in world coordinates

    contact: ContactConstraint
    x: body position
    q: body orientation
"""
function contact_location(contact::ContactConstraint{T,N,Nc,Cs}, x::AbstractVector{T},
    q::Quaternion{T}) where {T,N,Nc,Cs<:Contact{T,N}}
    model = contact.model
    return x + vector_rotate(model.collision.contact_origin, q) - model.collision.contact_normal' * model.collision.contact_radius
end

function contact_location(contact::ContactConstraint{T,N,Nc,Cs},
    body::Body) where {T,N,Nc,Cs<:Contact{T,N}}
    x = body.state.x2
    q = body.state.q2
    return contact_location(contact, x, q)
end

function contact_location(mechanism::Mechanism, contact::ContactConstraint)
    body = mechanism.bodies[findfirst(x -> x.id == contact.parent_id, mechanism.bodies)]
    return contact_location(contact, body)
end

function contact_location(mechanism::Mechanism)
    return [contact_location(mech, contact) for contact in mechanism.contacts]
end

"""
    string_location(contact, x, q)

    location of the attach point of the string in world coordinates

    collision: StringCollision
    x: body position
    q: body orientation
"""
function string_location(collision::StringCollision, x::AbstractVector{T},
    q::Quaternion{T}; relative::Symbol=:parent) where T
    origin = (relative == :parent) ? collision.origin_parent : collision.origin_child
    return x + vector_rotate(origin, q)
end

function string_location(collision::StringCollision, body::Body; relative::Symbol=:parent)
    x = body.state.x2
    q = body.state.q2
    return string_location(collision, x, q, relative=relative)
end

function string_location(mechanism::Mechanism, collision::StringCollision)
    pbody = get_body(mechanism, contact.parent_id)
    cbody = get_body(mechanism, contact.child_id)
    x_parent = string_location(collision, pbody, relative=:parent)
    x_child = string_location(collision, cbody, relative=:child)
    return x_parent, x_child
end


text:
subs(v::PolynomialLike, s::Vararg{<:Pair{<:Variable}}) = subs(v, s)
subs(x, s::Vararg{Pair{<:Variable}}) = x

(p::Polynomial)(s::Vararg{Pair{<:Variable}}) = subs(p, s)
(p::Term)(s::Vararg{Pair{<:Variable}}) = subs(p, s)
(p::Monomial)(s::Vararg{Pair{<:Variable}}) = subs(p, s)
(p::Variable)(s::Vararg{Pair{<:Variable}}) = subs(p, s)


text:
# This file demonstrates an implementation of FLAC support for Onda.jl.

using Onda, FLAC_jll, Test, Random, Dates

#####
##### FLACFormat
#####

"""
     FLACFormat(lpcm::LPCMFormat; sample_rate, level=5)
     FLACFormat(info::SamplesInfo; level=5)

Return a `FLACFormat<:AbstractLPCMFormat` instance that represents the
FLAC format corresponding to signals whose `file_format` field is `"flac"`.

The `sample_rate` keyword argument corresponds to `flac`'s `--sample-rate` flag,
while `level` corresponds to `flac`'s `--compression-level` flag.

Note that FLAC is only applicable for `info` where `1 <= channel_count(info) <= 8`
and `sizeof(info.sample_type) in (1, 2)`.

See https://xiph.org/flac/ for details about FLAC.

See also: [`Zstd`](@ref)
"""
struct FLACFormat{S} <: Onda.AbstractLPCMFormat
    lpcm::LPCMFormat{S}
    sample_rate::Int
    level::Int
    function FLACFormat(lpcm::LPCMFormat{S}; sample_rate, level=5) where {S}
        sizeof(S) in (1, 2) || throw(ArgumentError("bit depth must be 8 or 16"))
        1 <= lpcm.channel_count <= 8 || throw(ArgumentError("channel count must be between 1 and 8"))
        return new{S}(lpcm, sample_rate, level)
    end
end

FLACFormat(info::SamplesInfo; kwargs...) = FLACFormat(LPCMFormat(info); sample_rate=info.sample_rate,
                                                      kwargs...)

Onda.register_lpcm_format!(file_format -> file_format == "flac" ? FLACFormat : nothing)

file_format_string(::FLACFormat) = "flac"

function flac_raw_specification_flags(format::FLACFormat{S}) where {S}
    return (level="--compression-level-$(format.level)",
            endian="--endian=little",
            channels="--channels=$(format.lpcm.channel_count)",
            bps="--bps=$(sizeof(S) * 8)",
            sample_rate="--sample-rate=$(format.sample_rate)",
            is_signed=string("--sign=", S <: Signed ? "signed" : "unsigned"))
end

struct FLACStream{L<:Onda.LPCMStream} <: AbstractLPCMStream
    stream::L
end

function Onda.deserializing_lpcm_stream(format::FLACFormat, io)
    flags = flac_raw_specification_flags(format)
    cmd = flac() do flac_path
        return open(`$flac_path - --totally-silent -d --force-raw-format $(flags.endian) $(flags.is_signed)`, io)
    end
    return FLACStream(Onda.LPCMStream(format.lpcm, cmd))
end

function Onda.serializing_lpcm_stream(format::FLACFormat, io)
    flags = flac_raw_specification_flags(format)
    cmd = flac() do flac_path
        return open(`$flac_path --totally-silent $(flags) -`, io; write=true)
    end
    return FLACStream(Onda.LPCMStream(format.lpcm, cmd))
end

function Onda.finalize_lpcm_stream(stream::FLACStream)
    close(stream.stream.io)
    wait(stream.stream.io)
    return true
end

Onda.deserialize_lpcm(stream::FLACStream, args...) = deserialize_lpcm(stream.stream, args...)

Onda.serialize_lpcm(stream::FLACStream, args...) = serialize_lpcm(stream.stream, args...)

function Onda.deserialize_lpcm(format::FLACFormat, bytes, args...)
    stream = deserializing_lpcm_stream(format, IOBuffer(bytes))
    results = deserialize_lpcm(stream, args...)
    finalize_lpcm_stream(stream)
    return results
end

function Onda.serialize_lpcm(format::FLACFormat, samples::AbstractMatrix)
    io = IOBuffer()
    stream = serializing_lpcm_stream(format, io)
    serialize_lpcm(stream, samples)
    finalize_lpcm_stream(stream)
    return take!(io)
end

#####
##### tests
#####

saws(info, duration) = [(j + i) % 100 * info.sample_resolution_in_unit for
                        i in 1:channel_count(info), j in 1:sample_count(info, duration)]

if VERSION >= v"1.1.0"
    @testset "FLAC example" begin
        info = SamplesInfo(; kind="test", channels=["a", "b", "c"],
                           sample_unit="unit",
                           sample_resolution_in_unit=0.25,
                           sample_offset_in_unit=-0.5,
                           sample_type=Int16,
                           sample_rate=50)
        data = saws(info, Minute(3))
        samples = encode(Samples(data, info, false))
        fmt = FLACFormat(info)

        bytes = serialize_lpcm(fmt, samples.data)
        @test deserialize_lpcm(fmt, bytes) == samples.data
        @test deserialize_lpcm(fmt, bytes, 99) == view(samples.data, :, 100:size(samples.data, 2))
        @test deserialize_lpcm(fmt, bytes, 99, 201) == view(samples.data, :, 100:300)

        io = IOBuffer()
        stream = serializing_lpcm_stream(fmt, io)
        serialize_lpcm(stream, samples.data)
        @test finalize_lpcm_stream(stream)
        seekstart(io)
        stream = deserializing_lpcm_stream(fmt, io)
        @test deserialize_lpcm(stream) == samples.data
        finalize_lpcm_stream(stream) && close(io)

        io = IOBuffer(bytes)
        stream = deserializing_lpcm_stream(fmt, io)
        @test deserialize_lpcm(stream, 49, 51) == view(samples.data, :, 50:100)
        @test deserialize_lpcm(stream, 49, 51) == view(samples.data, :, 150:200)
        @test deserialize_lpcm(stream, 9) == view(samples.data, :, 210:size(samples.data, 2))
        finalize_lpcm_stream(stream) && close(io)
    end
else
    @warn "This example may be broken on Julia versions lower than v1.1 due to https://github.com/JuliaLang/julia/issues/33117"
end


text:
using RCall
## run library in R
@rlibrary IBrokers
## run TimeSeries
using TimeSeries
using DataFrames

symbol = "GBL"
exch = "DTB"
expiry = "201809"
currency = "EUR"
duration = "1 M"
barSize = "30 mins"

function getConnectIB()

    R"""
    library(IBrokers)
    tws <- twsConnect()
    """

end

function getIB(symbol, exch, expiry, currency, barSize, duration)

    R"""
    contract <- twsFuture(symbol = $symbol, exch = $exch , expiry = $expiry, currency = $currency)
    data <- reqHistoricalData(tws, contract, barSize = $barSize, duration = $duration)['T08:00/T18:00']
    """

    data = TimeArray(rcopy(R"index(as.xts(data))"), rcopy(R"data"))
    data.colnames[:] = ["Open", "High", "Low", "Close", "Volume", "WAP", "hasGaps", "Count"]

    return data

end

function getIB2(symbol, endDateTime, exch, expiry, currency, barSize, duration)

    #RCall for IBrokers
    R"""
    contract <- twsFuture(symbol = $symbol, exch = $exch , expiry = $expiry, currency = $currency)
    data <- reqHistoricalData(tws, contract, endDateTime = $endDateTime, barSize = $barSize, duration = $duration)['T08:00/T18:00']
    """

    #jData DataFrame, with timestamps first
    jData = DataFrame(TimeStamp = rcopy(R"index(as.xts(data))"))

    #colnames for Julia DataFrame
    cNames = ["Open", "High", "Low", "Close", "Volume", "WAP", "hasGaps", "Count"]

    for i in 1:length(cNames)

        # iterate through the return array of R and combine with initial DataFrame
        jData[Symbol(cNames[i])] = rcopy(R"data")[:,i]

    end

    return jData

end

function getDailyIB(symbol, exch, expiry, currency, barSize, duration)

    R"""
    contract <- twsFuture(symbol = $symbol, exch = $exch , expiry = $expiry, currency = $currency)
    data <- reqHistoricalData(tws, contract, barSize = $barSize, duration = $duration)
    """

    data = TimeArray(Date.(rcopy(R"index(as.xts(data))")), rcopy(R"data"))
    data.colnames[:] = ["Open", "High", "Low", "Close", "Volume", "WAP", "hasGaps", "Count"]

    return data

end

# #establish R connectin to IB, chech this by R>tws1
#
# rcopy(R"index(as.xts(data))")
# rcopy(data)
#
#
# rcopy(as.data.frame(data))
#
# dates = @rget df
#
#
# data = rcopy(data)
#
# rcopy(df)
#
# BarTimes = (@rget df)
#
# DateTime(BarTimes, "yyyy-mm-dd HH:MM:SS")


text:
"""
Unit tests for quicksort.
"""

using DataStructures, Random, Test
include("../src/sorting.jl")

@testset "Quicksort Integer Functions" begin
@test pow2_ceil(1) == 1
@test pow2_ceil(2) == 2
@test pow2_ceil(3) == 4
@test pow2_ceil(5) == 8
@test pow2_ceil(8) == 8

@test Θ(0) == 0
@test Θ(1) == 1
@test Θ(2) == 1

@test ceil_div(3, 2) == 2
@test ceil_div(4, 2) == 2

@test seek_last_less_than([1, 2, 2, 3, 4, 1, 2, 2, 3, 4], 3, 0, 5, true) == 3
@test seek_last_less_than([1, 2, 2, 3, 4, 1, 2, 2, 3, 4], 3, 5, 10, true) == 8
end

function test_batch_partition(T, N, lo, hi, seed)
    my_range = lo + 1 : hi
    Random.seed!(seed)
    original = rand(T, N)
    A = CuArray(original)

    #block_N = ceil_div(hi - lo, block_dim)
    pivot = rand(original[my_range])
    sums = CuArray(zeros(Int32, ceil_div(hi - lo, 32)))
    block_N, block_dim = -1, -1

    function get_config(kernel)
        get_shmem(threads) = threads * (sizeof(Int32) + max(4, sizeof(T)))

        fun = kernel.fun
        config = launch_configuration(fun, shmem=threads->get_shmem(threads), max_threads=1024)

        threads = pow2_floor(config.threads)
        blocks = ceil(Int, (hi - lo) ./ threads)
        block_N = blocks
        block_dim = threads
        @assert block_dim >= 32 "This test assumes block size can be >= 32"
        shmem = get_shmem(threads)
        return (threads=threads, blocks=blocks, shmem=shmem)
    end
    @cuda config=get_config partition_batches_kernel(A, pivot, sums, lo, hi - lo, Val(true))
    synchronize()

    post_sort = Array(A)
    post_sums = Array(sums)

    sort_match = true
    sums_match = true

    for block in 1:block_N
        block_range = lo + 1 + (block - 1) * block_dim: min(hi, lo + block * block_dim)
        temp = original[block_range]
        #= this shows that batch partitioning is a stable sort where key for
        each value v is whether v > or <= pivot =#
        expected_sort = vcat(filter(x -> x < pivot, temp), filter(x -> x >= pivot, temp))
        sort_match &= post_sort[block_range] == expected_sort
        sums_match &= post_sums[block] == count(x -> x >= pivot, temp)
    end

    @test sort_match
    @test sums_match
end

@testset "Quicksort batch partition" begin
test_batch_partition(Int8, 10000, 2000, 6000, 0)
test_batch_partition(Int8, 10000, 2000, 6000, 1)
test_batch_partition(Int8, 10000000, 0, 10000000, 0)
test_batch_partition(Int8, 10000000, 5000, 500000, 0)
test_batch_partition(Int8, 10000, 0, 10000, 0)
test_batch_partition(Int8, 10000, 2000, 6000, 0)
test_batch_partition(Int8, 10000, 2000, 6000, 1)
test_batch_partition(Int8, 10000000, 0, 10000000, 0)
test_batch_partition(Int8, 10000000, 5000, 500000, 0)

test_batch_partition(Float32, 10000, 0, 10000, 0)
test_batch_partition(Float32, 10000, 2000, 6000, 0)
test_batch_partition(Float32, 10000, 2000, 6000, 1)
test_batch_partition(Float32, 10000000, 0, 10000000, 0)
test_batch_partition(Float32, 10000000, 5000, 500000, 0)
test_batch_partition(Float32, 10000, 0, 10000, 0)
test_batch_partition(Float32, 10000, 2000, 6000, 0)
test_batch_partition(Float32, 10000, 2000, 6000, 1)
test_batch_partition(Float32, 10000000, 0, 10000000, 0)
test_batch_partition(Float32, 10000000, 5000, 500000, 0)
end

function test_consolidate_kernel(vals, pivot, my_floor, L, N_t, N_b, b_sums, dest)
    i = threadIdx().x
    p = consolidate_batch_partition(vals, pivot, my_floor, L, N_t, N_b, i, b_sums)
    if i == 1
        dest[1] = p
    end
    return nothing
end

function test_consolidate_partition(T, N, lo, hi, seed)
    # assuming partition_batches works, we can validate consolidate by
    # checking that together they partition a large domain
    my_range = lo + 1 : hi
    Random.seed!(seed)
    original = rand(T, N)
    A = CuArray(original)
    pivot = rand(original[my_range])
    block_dim = block_N = -1
    sums = CuArray(zeros(Int32, ceil_div(hi - lo, 32)))

    function get_config(kernel)
        get_shmem(threads) = threads * (sizeof(Int32) + max(4, sizeof(T)))

        fun = kernel.fun
        config = launch_configuration(fun, shmem=threads->get_shmem(threads), max_threads=1024)

        threads = pow2_floor(config.threads)
        blocks = ceil(Int, (hi - lo) ./ threads)
        block_N = blocks
        block_dim = threads
        @assert block_dim >= 32 "This test assumes block size can be >= 32"

        shmem = get_shmem(threads)
        return (threads=threads, blocks=blocks, shmem=shmem)
    end

    @cuda config=get_config partition_batches_kernel(A, pivot, sums, lo, hi - lo, Val(true))
    synchronize()
    dest = CuArray(zeros(Int32, 1))
    @cuda threads=block_dim test_consolidate_kernel(A, pivot, lo, hi - lo, block_dim, block_N, sums, dest)
    synchronize()
    partition = Array(dest)[1]
    temp = original[my_range]
    post_sort = Array(A)
    #= consolidation is a highly unstable sort (again, by pivot comparison as
    the key) so we compare by counting each element =#
    cc(x) = x |> counter |> collect |> sort
    @test all(post_sort[lo + 1 : partition] |> cc .== filter(x -> x < pivot, temp) |> cc)
    @test all(post_sort[partition + 1 : hi] |> cc .== filter(x -> x >= pivot, temp) |> cc)
end

function test_consolidate_partition(T, N, lo, hi, seed, block_dim)
    # assuming partition_batches works, we can validate consolidate by
    # checking that together they partition a large domain
    my_range = lo + 1 : hi
    Random.seed!(seed)
    original = rand(T, N)
    A = CuArray(original)
    block_N = ceil_div(hi - lo, block_dim)
    pivot = rand(original[my_range])
    sums = CuArray(zeros(Int32, block_N))

    @cuda blocks=block_N threads=block_dim shmem=block_dim *(sizeof(Int32) + sizeof(T)) partition_batches_kernel(A, pivot, sums, lo, hi - lo, Val(true))
    synchronize()
    dest = CuArray(zeros(Int32, 1))
    @cuda threads=block_dim test_consolidate_kernel(A, pivot, lo, hi - lo, block_dim, block_N, sums, dest)
    synchronize()
    partition = Array(dest)[1]
    temp = original[my_range]
    post_sort = Array(A)
    #= consolidation is a highly unstable sort (again, by pivot comparison as
    the key) so we compare by counting each element =#
    cc(x) = x |> counter |> collect |> sort
    @test all(post_sort[lo + 1 : partition] |> cc .== filter(x -> x < pivot, temp) |> cc)
    @test all(post_sort[partition + 1 : hi] |> cc .== filter(x -> x >= pivot, temp) |> cc)
end

@testset "Quicksort consolidate partition" begin
test_consolidate_partition(Int8, 10000, 0, 10000, 0, 16)
test_consolidate_partition(Int8, 10000, 0, 10000, 0, 32)
test_consolidate_partition(Int8, 10000, 0, 10000, 0, 64)
test_consolidate_partition(Int8, 10000, 9, 6333, 0, 16)
test_consolidate_partition(Int8, 10000, 9, 6333, 0, 32)
test_consolidate_partition(Int8, 10000, 9, 6333, 0, 64)
test_consolidate_partition(Int8, 10000, 129, 9999, 0, 16)
test_consolidate_partition(Int8, 10000, 129, 9999, 0, 32)
test_consolidate_partition(Int8, 10000, 129, 9999, 0, 64)
test_consolidate_partition(Int8, 10000, 0, 10000, 1, 16)
test_consolidate_partition(Int8, 10000, 0, 10000, 2, 32)
test_consolidate_partition(Int8, 10000, 0, 10000, 3, 64)
test_consolidate_partition(Int8, 10000, 9, 6333, 4, 16)
test_consolidate_partition(Int8, 10000, 9, 6333, 5, 32)
test_consolidate_partition(Int8, 10000, 9, 6333, 6, 64)
test_consolidate_partition(Int8, 10000, 129, 9999, 7, 16)
test_consolidate_partition(Int8, 10000, 129, 9999, 8, 32)
test_consolidate_partition(Int8, 10000, 129, 9999, 9, 64)
test_consolidate_partition(Int8, 10000, 3329, 9999, 10, 16)
test_consolidate_partition(Int8, 10000, 3329, 9999, 11, 32)
test_consolidate_partition(Int8, 10000, 3329, 9999, 12, 64)

test_consolidate_partition(Int8, 10000, 0, 10000, 0)
test_consolidate_partition(Int16, 10000, 0, 10000, 0)
test_consolidate_partition(Int32, 10000, 0, 10000, 0)
test_consolidate_partition(Int64, 10000, 0, 10000, 0)

test_consolidate_partition(Float16, 10000, 0, 10000, 0)
test_consolidate_partition(Float32, 10000, 0, 10000, 1)
test_consolidate_partition(Float64, 10000, 0, 10000, 2)
end

function test_quicksort(T, f, N)
    a = map(x -> T(f(x)), 1:N)
    c = CuArray(a)
    quicksort!(c)
    @test Array(c) == sort(a)
end

@testset "Quicksort" begin
# test pre-sorted
test_quicksort(Int, x -> x, 4000000)
test_quicksort(Int32, x -> x, 4000000)
test_quicksort(Float64, x -> x, 4000000)
test_quicksort(Float32, x -> x, 4000000)

# test reverse sorted
test_quicksort(Int, x -> -x, 4000000)
test_quicksort(Int32, x -> -x, 4000000)
test_quicksort(Float64, x -> -x, 4000000)
test_quicksort(Float32, x -> -x, 4000000)

# test random arrays
Random.seed!(0)

test_quicksort(Int, x -> rand(Int), 10000)
test_quicksort(Int32, x -> rand(Int32), 10000)
test_quicksort(Int8, x -> rand(Int8), 10000)
test_quicksort(Float64, x -> rand(Float64), 10000)
test_quicksort(Float32, x -> rand(Float32), 10000)
test_quicksort(Float16, x -> rand(Float16), 10000)

test_quicksort(Int, x -> rand(Int), 4000000)
test_quicksort(Int32, x -> rand(Int32), 4000000)
test_quicksort(Int8, x -> rand(Int8), 4000000)
test_quicksort(Float64, x -> rand(Float64), 4000000)
test_quicksort(Float32, x -> rand(Float32), 4000000)
test_quicksort(Float16, x -> rand(Float16), 4000000)
end

max_stars_repo_name: StatisticalMice/experiments

max_issues_repo_name: StatisticalMice/experiments

max_forks_repo_name: StatisticalMice/experiments


text:
############################################################################################################################

# This file contains tests for EdgeIter

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

function edge_iter_sanity_test(ne::Int, eit::EdgeIter)
   @test eit == eit
   @test length(eit) == ne
   @test size(eit) == (ne,)
   @test eit == copy(eit)
   @test deepcopy(eit) == deepcopy(eit)
   @test issorted(eit) == true
   @test eltype(eit) == EdgeID
   @test EdgeIter(collect(eit)) == eit
end

function edge_iter_iteration_test(ne::Int, eit::EdgeIter)
   @test start(eit) == 1
   @test endof(eit) == ne
   @test eachindex(eit) == 1 : ne
   @test [e for e in eit] == collect(eit) == map(EdgeID, eit.us, eit.vs)
end

function edge_iter_getindex_test(ne::Int, eit::EdgeIter)
   @test eit[1] == first(eit)
   @test eit[1:ne] == eit
   @test eit[:] == eit
end

@testset "EdgeIter Construction" begin
   x = completeindxs(10)

   eit1 = EdgeIter(x)
   edge_iter_sanity_test(90, eit1)
   edge_iter_iteration_test(90, eit1)
   edge_iter_getindex_test(90, eit1)

   eit2 = EdgeIter(x, 1)
   edge_iter_sanity_test(9, eit2)
   edge_iter_iteration_test(9, eit2)
   edge_iter_getindex_test(9, eit2)
   @test eit2 == eit1[1:9]

   eit3 = EdgeIter(x, 1:5)
   edge_iter_sanity_test(45, eit3)
   edge_iter_iteration_test(45, eit3)
   edge_iter_getindex_test(45, eit3)
end


text:
module Unmarshal

# package code goes here

# Helper function
function prettyPrint(verboseLvl, str)
    tabs = ""
    for cntr = 1:verboseLvl
        tabs = tabs * "\t"
    end
    println("$(tabs)$(str)")
end

export unmarshal # returns a reconstructed variable from a JSON parsed string

using Requires

using JSON
import Missings: Missing, missing
import Nullables: Nullable

function __init__()
    @require LazyJSON="fc18253b-5e1b-504c-a4a2-9ece4944c004" include("lazyjson.jl")
end

unmarshal(::Type{Any}, x::String, verbose :: Bool = false, verboseLvl :: Int = 0) = x
unmarshal(::Type{Any}, x, verbose :: Bool = false, verboseLvl :: Int = 0) = x
function unmarshal(::Type{Any}, xs::Union{Vector{E}, AbstractArray}, verbose :: Bool = false, verboseLvl :: Int = 0) where E
    if verbose
        prettyPrint(verboseLvl, "Any{$(eltype(xs))}")
        prettyPrint(verboseLvl, "List")
        verboseLvl += 1
    end
    [(unmarshal(eltype(xs), x, verbose, verboseLvl) for x in xs)...]
end

function unmarshal(DT :: Type{String}, parsedJson :: String, verbose :: Bool = false, verboseLvl :: Int = 0)
    if verbose
        prettyPrint(verboseLvl, "$(DT) (String)")
        verboseLvl += 1
    end
    parsedJson
end

function unmarshal(::Type{Vector{E}}, parsedJson::Union{Vector, AbstractArray}, verbose :: Bool = false, verboseLvl :: Int = 0) where E
    if verbose
        prettyPrint(verboseLvl, "Vector{$E}")
        verboseLvl += 1
    end

    [(unmarshal(E, field, verbose, verboseLvl) for field in parsedJson)...]
end

unmarshal(::Type{Array{E}}, xs::Union{Vector, AbstractArray}, verbose :: Bool = false, verboseLvl :: Int = 0) where E = unmarshal(Vector{E}, xs, verbose, verboseLvl)

function unmarshal(::Type{Array{E, N}}, parsedJson::Union{Vector, AbstractArray}, verbose :: Bool = false, verboseLvl :: Int = 0) where {E, N}
    if verbose
        prettyPrint(verboseLvl, "Array{$E, $N}")
        verboseLvl += 1
    end

    cat((unmarshal(Array{E,N-1}, x, verbose, verboseLvl) for x in parsedJson)..., dims=N)
end

unmarshal(::Type{Array{E, N} where E}, xs::Union{Vector, AbstractArray}, verbose :: Bool = false, verboseLvl :: Int = 0) where N = unmarshal(Array{Any, N}, xs, verbose, verboseLvl)
unmarshal(::Type{Array}, xs::Union{Vector,AbstractArray}, verbose :: Bool = false, verboseLvl :: Int = 0) where N = unmarshal(Vector{Any}, xs, verbose, verboseLvl)

"""
    unmarshal(T, dict[, verbose[, verboselvl]])

Reconstructs an object of Type T using the dictionary output of a `JSON.parse` or now 'LazyJSON.parse`.

Set verbose `true` to get debug information about how the data hierarchy is unmarshalled. This might be useful to track down parsing errors and/or mismatches between the JSON object and the Type definition.

# Example

```jldoctest
julia> using JSON

julia> var = randn(Float64, 5);  # Should work for most other variations of types you can think of

julia> unmarshal(typeof(var), JSON.parse(JSON.json(var)) ) == var
true

or

julia> using LazyJSON

julia> var = randn(Float64, 5);  # Should work for most other variations of types you can think of

julia> unmarshal(typeof(var), LazyJSON.parse(JSON.json(var)) ) == var
true
```

"""
function unmarshal(DT :: Type, parsedJson :: AbstractDict, verbose :: Bool = false, verboseLvl :: Int = 0)
    if verbose
            prettyPrint(verboseLvl, "$(DT) AbstractDict")
        verboseLvl += 1
    end

    if !isconcretetype(DT)
        throw(ArgumentError("Cannot unmarshal a non-leaf type $(DT) without a custom specialization"))
    end

    tup = ()
    for iter in fieldnames(DT)
        DTNext = fieldtype(DT,iter)
        if verbose
           prettyPrint(verboseLvl-1, "\\--> $(iter) <: $(DTNext) ")
        end

        if !haskey(parsedJson, string(iter)) || parsedJson[string(iter)] === nothing
            # check whether DTNext is compatible with any scheme for missing values
            val = if DTNext <: Nullable
                DTNext()
            elseif Missing <: DTNext
                missing
            elseif Nothing <: DTNext
                Nothing()
            elseif !haskey(parsedJson, string(iter))
                throw(ArgumentError("Key $(string(iter)) is missing from the structure $DT, and field is neither Nullable nor Missings nor Nothing-compatible"))
            else
                throw(ArgumentError("Key $(string(iter)) is null, but the field is neither Nullable nor Missings nor Nothing-compatible"))
            end
        else
            val = unmarshal( DTNext, parsedJson[string(iter)], verbose, verboseLvl)
        end

        tup = (tup..., val)
    end

    DT(tup...)
end

function unmarshal(DT :: Type{Pair{TF, TS}}, parsedJson :: AbstractDict, verbose :: Bool = false, verboseLvl :: Int = 0) where {TF, TS}
    if verbose
          prettyPrint(verboseLvl, "Pair $(DT) AbstractDict")
          verboseLvl += 1
    end

    if (length(keys(parsedJson)) > 1)
          @warn "Expected a single pair, but found a multi-entry dictionary, just using the first key: $(collect(keys(parsedJson))[1])"
    end
    firstVal = (collect(keys(parsedJson))[1]) #, verbose, verboseLvl)
    secondVal = (parsedJson[firstVal]) #, verbose, verboseLvl)
#    @show firstVal, secondVal

    if !isa(firstVal, TF)
        try
           firstVal = TF(firstVal)
        catch ex
           firstVal = unmarshal(TF, JSON.parse(firstVal), verbose, verboseLvl)
        end
    end

    if !isa(secondVal, TS)
        try
           secondVal = TS(secondVal)
        catch ex
           throw(ArgumentError("Error trying to convert value $(secondVal) of type $(typeof(secondVal)) to a $(TS), please provide a conversion"))
        end
    end

#    @show firstVal, secondVal

    (firstVal => secondVal)
end

function unmarshal(::Type{T}, parsedJson :: Union{Vector, AbstractArray}, verbose :: Bool = true, verboseLvl :: Int = 0) where {T <: Tuple}
    if verbose
        prettyPrint(verboseLvl, "Tuple: $T")
        verboseLvl += 1
    end
    len = length(parsedJson) # fallback value
    try
        # If we have concrete lengths, rather use data type and ignore length of JSON
        len = fieldcount(T)
    catch ex
        # If we do not have concrete lengths, use length of JSON
    end
    ((unmarshal(fieldtype(T, i), parsedJson[i], verbose, verboseLvl) for i in 1:len)...,)
end

function _unmarshal(::Type{T}, key :: Symbol, parsedJson :: AbstractDict, verbose :: Bool = false, verboseLvl :: Int = 0) where T
    if verbose
        prettyPrint(verboseLvl - 1, "\\--> $(key) ")
    end
    unmarshal(T, parsedJson[string(key)], verbose, verboseLvl)
end

function unmarshal(::Type{NamedTuple{NS, TS}}, parsedJson :: AbstractDict, verbose :: Bool = false, verboseLvl :: Int = 0) where {NS, TS}
    if verbose
        prettyPrint(verboseLvl, "NamedTuple{$NS,$TS}")
        verboseLvl += 1
    end

    NamedTuple{NS}(_unmarshal(T, key, parsedJson, verbose, verboseLvl) for (i, (T, key)) in enumerate(zip(TS.parameters, NS)))
end
unmarshal(::Type{NamedTuple{NS}}, parsedJson :: AbstractDict, verbose :: Bool = false, verboseLvl :: Int = 0) where NS = unmarshal(NamedTuple{NS, NTuple{length(NS), Any}}, parsedJson, verbose, verboseLvl)
unmarshal(::Type{NamedTuple}, parsedJson :: AbstractDict, verbose :: Bool = false, verboseLvl :: Int = 0) = unmarshal(NamedTuple{(Symbol.(keys(parsedJson))...,)}, parsedJson, verbose, verboseLvl)

function unmarshal(DT :: Type{T}, parsedJson :: AbstractDict, verbose :: Bool = false, verboseLvl :: Int = 0) where T <: Dict
    if verbose
        prettyPrint(verboseLvl, "$(DT) Dict ")
        verboseLvl += 1
    end
    val = DT()
    for iter in keys(parsedJson)
        if verbose
           prettyPrint(verboseLvl - 1, "\\--> $(iter) $(valtype(DT))")
        end
        tmp = unmarshal(valtype(DT), parsedJson[iter], verbose, verboseLvl)
        if keytype(DT) <: AbstractString
            val[iter] = tmp
        else
            try
                val[unmarshal(keytype(DT),JSON.parse(iter),verbose, verboseLvl)] = tmp # Use JSON.parse and Unmarshal to cast from type of iter to ketype(DT)
            catch ex
                val[keytype(DT)(iter)] = tmp  # Try direct casting, which will hopefully generate a readable enough exception error if it fails
            end
        end
    end
    val
end

unmarshal(::Type{T}, x::Number, verbose :: Bool = false, verboseLvl :: Int = 0) where T <: Number = T(x)
unmarshal(::Type{Nullable{T}}, x, verbose :: Bool = false, verboseLvl :: Int = 0) where T = Nullable(unmarshal(T, x))
unmarshal(::Type{Nullable{T}}, x::Nothing, verbose :: Bool = false, verboseLvl :: Int = 0) where T = Nullable{T}()
unmarshal(::Type{Union{T,Missing}}, x, verbose :: Bool = false, verboseLvl :: Int = 0) where T = unmarshal(T, x, verbose, verboseLvl)
unmarshal(::Type{Union{T,Nothing}}, x::Nothing, verbose :: Bool = false, verboseLvl :: Int = 0) where T = nothing
unmarshal(::Type{Union{T,Nothing}}, x::T, verbose :: Bool = false, verboseLvl :: Int = 0) where T = unmarshal(T, x, verbose, verboseLvl)

function unmarshal(T::Type, x, verbose :: Bool = false, verboseLvl :: Int = 0)
    try
        T(x)
    catch
        throw(ArgumentError("no unmarshal function defined to convert $(typeof(x)) to $(T); consider providing a specialization"))
    end
end

end # module


text:
module AdventOfCode

end


text:
# Utils
function module_dir()
    return joinpath(@__DIR__, "..")
end

# Activate package
using Pkg; Pkg.activate(module_dir())

using LinearAlgebra
using LyceumBase
using LyceumMuJoCo
using LyceumMuJoCoViz
using MuJoCo
using Plots
using JLD2
using Random

include("methods.jl")

function mj_astronaut_simulation(; timestep=0.01, tsim=1.0, tctrl=1.0, seed::Int=0, ϵ=1e-14,
		control_amplitude=0.05, record_traj::Bool=true, record_ener::Bool=false)

	jm, jd, mjsim = mj_model("astronaut.xml", timestep=timestep)
    N = Int(floor(tsim/timestep))
	tcompute = 0.0
	ener0 = 0.0
	ener_traj = []

	traj = zeros(83,0)
	for i = 1:N
		jd.ctrl .= (i*timestep < tctrl) * control_amplitude * rand(21)
		tcompute += @elapsed mj_step(jm, jd)
		record_traj && (traj = hcat(traj, deepcopy(getstate(mjsim))))
		if !(i*timestep < tctrl) && ((i-1)*timestep < tctrl)
			# first step after the end of control inputs
			ener0 = energy(jm, jd)
		end
		if record_ener && !(i*timestep < tctrl) && (i % Int(floor(1/timestep)) == 0)
			push!(ener_traj, energy(jm, jd) - ener0)
		end
	end
	fail = norm(jd.qvel,Inf) < 1e-10 # when Mujoco fails it returns a zero qvel
	ener1 = energy(jm, jd)
	ener = ener1 - ener0
	plin = norm(linear_momentum(jm, jd))
	pang = norm(angular_momentum(jm, jd))
    return tcompute, plin, pang, ener, fail, traj, ener_traj
end

function mj_batch_astronaut_simulation(; Nsim=1, timestep=0.01, tsim=1.0, tctrl=1.0, seed::Int=0, ϵ=1e-14,
		control_amplitude=0.05, record_traj::Bool=false)
	speed = 0.0
	plin = 0.0
	pang = 0.0
	ener = 0.0
	fail = 0.0
	for i = 1:Nsim
		tcompute, pl, pa, en, fa, _ = mj_astronaut_simulation(; timestep=timestep, tsim=tsim, tctrl=tctrl, seed=seed, ϵ=ϵ,
			control_amplitude=control_amplitude, record_traj=record_traj)
		speed += tsim / tcompute
		plin += pl
		pang += pa
		ener += abs(en)
		fail += fa
	end
	speed /= Nsim
	plin /= Nsim
	pang /= Nsim
	ener /= Nsim
	fail /= Nsim
    return speed, plin, pang, ener, fail
end

function mj_benchmark_momentum(timestep; Nsim=1, tsim=1.0, tctrl=1.0, seed::Int=0, ϵ=1e-14,
		control_amplitude=0.05, record_traj::Bool=false)
	N = length(timestep)
	speed = zeros(N)
	plin = zeros(N)
	pang = zeros(N)
	ener = zeros(N)
	fail = zeros(N)
	for i = 1:N
		@show timestep[i]
		speed[i], plin[i], pang[i], ener[i], fail[i] = mj_batch_astronaut_simulation(; Nsim=Nsim, timestep=timestep[i],
			tsim=tsim, tctrl=tctrl, seed=seed, ϵ=ϵ, control_amplitude=control_amplitude, record_traj=record_traj)
	end
    return speed, plin, pang, ener, fail
end

speed, plin, pang, Δener, fail, traj, ener_traj = mj_astronaut_simulation(tsim=1, tctrl=1.0,
	timestep=0.1, control_amplitude=0.02)
speed, plin, pang, Δener, fail, traj, ener_traj = mj_astronaut_simulation(tsim=100, tctrl=1.0,
	timestep=0.1, control_amplitude=0.02)
LyceumMuJoCoViz.visualize(mjsim, trajectories=traj)

################################################################################
# momentum
################################################################################
# mj_batch_astronaut_simulation(Nsim=10)
timestep= [0.1, 0.03, 0.01, 0.003, 0.001, 0.0003, 0.0001, 0.00003]
speed, plin, pang, ener, fail = mj_benchmark_momentum(timestep, Nsim=3, control_amplitude=0.05)

# Saving results
# path = joinpath(module_dir(), "results/astronaut_momentum.jld2")
# jldsave(path, speed=speed, plin=plin, pang=pang)

plt = plot(layout=(1,2), legend=false, xlims=(1e-2,1e6), ylims=(1e-15,1e6))
plot!(plt[1,1], speed, plin, axis=:log, yaxis=:flip, linewidth=3, markersize=6)
plot!(plt[1,2], speed, pang, axis=:log, yaxis=:flip, linewidth=3, markersize=6)
scatter!(plt[1,1], speed, plin, axis=:log, yaxis=:flip, linewidth=3, markersize=6)
scatter!(plt[1,2], speed, pang, axis=:log, yaxis=:flip, linewidth=3, markersize=6)

################################################################################
# energy
################################################################################
timestep= [0.1, 0.01, 0.001]
ener_traj = []
for i = 1:3
	speed, plin, pang, Δener, fail, traj, ener_t = mj_astronaut_simulation(tsim=100, tctrl=1.0,
		timestep=timestep[i], control_amplitude=0.02, record_traj=false, record_ener=true)
	push!(ener_traj, ener_t)
end

plt = plot(layout=(3,1), legend=false)#, xlims=(0,100))
for i = 1:3
	plot!(plt[i,1],
		range(0,stop=100,length=length(ener_traj[i])),
		ener_traj[i], linewidth=3, markersize=6)
end
display(plt)

# Saving results
path = joinpath(module_dir(), "results/astronaut_energy.jld2")
jldsave(path, ener_traj=ener_traj)

################################################################################
# Demo
################################################################################

timestep=0.01
jm, jd, mjsim = mj_model("astronaut.xml", timestep=0.01)
tsim = 1.0
N = Int(floor(tsim/timestep))

# jd.qpos[3] += 10.0
# Random.seed!(100)
# jd.qvel .= 1.0 * (rand(27) .- 0.5)

traj = zeros(83, 0)
E = zeros(N)
Plin = zeros(N)
Pang = zeros(N)
tcompute = 0.0
for i = 1:N
    jd.ctrl .= 0.3 * [0; 1; zeros(19)]
    tcompute += @elapsed mj_step(jm, jd);
    traj = hcat(traj, deepcopy(getstate(mjsim)))
    E[i] = energy(jm, jd)
	Plin[i] = norm(linear_momentum(jm, jd))
	# Pang[i] = angular_momentum(jm, jd)
end
tsim / tcompute
plot(E)
plot(Plin)
plot(Pang)

LyceumMuJoCoViz.visualize(mjsim, trajectories=traj)

jd.qpos
jd.qvel


max_stars_repo_path: test/1 Algebraic functions/1.2 Trinomial products/1.2.3 General/1.2.3.5 P(x) (d x)^m (a+b x^n+c x^(2 n))^p.jl


text:
## Packages.
using LoopVectorization
using BenchmarkTools

## Functions.
function solve_tdma_kernel!(
    p, work3d, work2d,
    a, b, c,
    itot, jtot, ktot)

    @tturbo for j in 1:jtot
        for i in 1:itot
            work2d[i, j] = b[i, j, 1]
            p[i, j, 1] /= work2d[i, j]
        end
    end

    for k in 2:ktot
        @tturbo for j in 1:jtot
            for i in 1:itot
                work3d[i, j, k] = c[k-1] / work2d[i, j]
                work2d[i, j] = b[i, j, k] - a[k]*work3d[i, j, k]
                p[i, j, k] -= a[k] * p[i, j, k-1]
                p[i, j, k] /= work2d[i, j]
            end
        end
    end

    for k in ktot-1:-1:1
        @tturbo for j in 1:jtot
            for i in 1:itot
                p[i, j, k] -= work3d[i, j, k+1] * p[i, j, k+1]
            end
        end
    end
end

function solve_tdma_kernel_2!(
    p, work3d, work2d,
    a, b, c,
    itot, jtot, ktot)

    @tturbo for j in 1:jtot
        for i in 1:itot
            work2d[i, j] = b[i, j, 1]
            p[i, j, 1] /= work2d[i, j]
        end
    end

    for k in 2:ktot
        @tturbo for j in 1:jtot
            for i in 1:itot
                work3d[i, j, k] = c[k-1] / work2d[i, j]
                work2d[i, j] = b[i, j, k] - a[k]*work3d[i, j, k]
                p[i, j, k] -= a[k] * p[i, j, k-1]
                p[i, j, k] /= work2d[i, j]
            end
        end
    end

    for k in ktot-1:-1:1
        @tturbo for j in 1:jtot
            for i in 1:itot
                p[i, j, k] -= work3d[i, j, k+1] * p[i, j, k+1]
            end
        end
    end
end

## Settings.
itot = 512; jtot = 1; ktot = 512;
p_ref = rand(itot, jtot, ktot) .+ 1
work3d = zeros(itot, jtot, ktot)
work2d = zeros(itot, jtot)
a = rand(ktot) .+ 1
b = rand(itot, jtot, ktot) .+ 1
c = rand(ktot) .+ 1
p = copy(p_ref)
p2 = copy(p_ref)

## First solver.
# @btime solve_tdma_kernel!($p, $work3d, $work2d, $a, $b, $c, $itot, $jtot, $ktot)

## Second solver.
@btime solve_tdma_kernel_2!($p2, $work3d, $work2d, $a, $b, $c, $itot, $jtot, $ktot)


text:
# Create a finite state machine, using model supplied by the user.
function state_machine(model::Dict)

    fsm = StateMachine()

    # Write events and their associated states to the state machine's map
    function add(event::Dict)

        # Multiple source states are allowed
        if haskey(event, "from")
            if isa(event["from"], Array)
                from = event["from"]
            else
                from = [event["from"]]
            end
        end
        if !haskey(fsm.map, event["name"])
            fsm.map[event["name"]] = @compat Dict()
        end

        # Use "to" if it has been specified; otherwise, this is a no-op
        # transition (before and after states are the same)
        for f in from
            if haskey(event, "to")
                fsm.map[event["name"]][f] = event["to"]
            else
                fsm.map[event["name"]][f] = f
            end
        end
    end

    # Set up initial state and startup event
    if haskey(model, "initial")

        # Initial state can be specified as a string or Dict.
        # If the initial event has not been specified, default to "startup".
        if isa(model["initial"], AbstractString)
            initial = @compat Dict{AbstractString, Any}(
                "state" => model["initial"],
                "event" => "startup",
            )
        else
            initial = model["initial"]
            if !haskey(initial, "event")
                initial["event"] = "startup"
            end
        end

        # Add the startup event to the map
        add(@compat Dict{AbstractString, Any}(
            "name" => initial["event"],
            "from" => "none",
            "to" => initial["state"],
        ))
    end

    # Terminal (final) state
    if haskey(model, "terminal")
        fsm.terminal = model["terminal"]
    elseif haskey(model, "final")
        fsm.terminal = model["final"]
    end

    # Write user-specified events to the map
    if haskey(model, "events")
        for event in model["events"]
            add(event)
        end
    end

    # Set up the callable events (functions), which can be invoked by fire()
    for (name, minimap) in fsm.map
        fsm.actions[name] = () -> begin
            from = fsm.current
            to = haskey(minimap, from) ? minimap[from] : from
            if !fire(fsm, "can", name)
                error(name * " is not accessible from state " * fsm.current)
            elseif from == to
                after_event(fsm, name, from, to)
            else
                leave_state(fsm, name, from, to)
                fsm.current = to
                enter_state(fsm, name, from, to)
                change_state(fsm, name, from, to)
                after_event(fsm, name, from, to)
            end
        end
    end

    # Set up user-specified callbacks, if any were provided
    if haskey(model, "callbacks")
        for (name, callback) in model["callbacks"]
            fsm.actions[name] = callback
        end
    end

    # Fire (or defer firing of) the initial event
    if !haskey(initial, "defer") || !initial["defer"]
        fire(fsm, initial["event"])
    end

    fsm
end


text:
# Run code for previous steps with plotting turned off.
make_plots_orig_5 = isdefined(Main,:make_plots) ? make_plots : true
 make_plots = false
 include("neid_solar_5_calc_rvs_proj.jl")
 make_plots = make_plots_orig_5
 if make_plots
   using Plots
 end
using MultivariateStats

# Set parameters for plotting analysis
plt_order = 42
 plt_order_pix = 4500:5000

fm_perp = RvSpectML.compute_spectra_perp_doppler_shift(spectral_orders_matrix.flux,deriv_orders, ave_good_chunks_rvs )
  idx_good_obs = 1:length(ave_good_chunks_rvs)
  M = fit(PCA, fm_perp[:,idx_good_obs]; maxoutdim=10)
  pca_out = MultivariateStats.transform(M,fm_perp[:,idx_good_obs])
  frac_var_explained = 1.0.-cumsum(principalvars(M))./tvar(M)
  if make_plots
    plt = scatter(frac_var_explained, xlabel="Number of PCs", ylabel="Frac Variance Unexplained")
    display(plt)
  end
  frac_var_explained

if make_plots
  plt_order = 40
  # plt_order_pix = 3301:3800
  RvSpectML.plot_basis_vectors(order_grids, f_mean_orders, deriv_orders, M.proj, idx_plt = spectral_orders_matrix.chunk_map[plt_order][plt_order_pix] )
end

if make_plots
  RvSpectML.plot_basis_scores(plt_times, ave_good_chunks_rvs, pca_out )
end

if make_plots
  RvSpectML.plot_basis_scores_cor( ave_good_chunks_rvs, pca_out)
end


text:
let
    @test BayesNets.strip_arg(:symbol) == :symbol

    @test BayesNets.CPDs.paramcount(:symbol) == 1
    @test BayesNets.CPDs.paramcount(true) == 1
    @test BayesNets.CPDs.paramcount(false) == 1
    @test BayesNets.CPDs.paramcount(1.0) == 1
    @test BayesNets.CPDs.paramcount(Array{Float64}(3)) == 3
end


text:
using Base.Test

importall AES

ctx = aes_context(Array(Uint32,64),Array(Uint32,64),0)
encrypted = Array(Uint8,16)
decrypted = Array(Uint8,16)

println("AES - NIST FIPS PUB 197 test cases:")
plaintext = b"\x00\x11\x22\x33\x44\x55\x66\x77\x88\x99\xaa\xbb\xcc\xdd\xee\xff"

print(" C.1 AES-128 ")
key128 = b"\x00\x01\x02\x03\x04\x05\x06\x07\x08\x09\x0a\x0b\x0c\x0d\x0e\x0f"
aes_set_key(ctx, key128, 128)
aes_encrypt(ctx, plaintext, encrypted); @test encrypted==b"\x69\xc4\xe0\xd8\x6a\x7b\x04\x30\xd8\xcd\xb7\x80\x70\xb4\xc5\x5a"
aes_decrypt(ctx, encrypted, decrypted); @test decrypted==b"\x00\x11\x22\x33\x44\x55\x66\x77\x88\x99\xaa\xbb\xcc\xdd\xee\xff"
println("- Passed")

print(" C.2 AES-192 ")
key192 = b"\x00\x01\x02\x03\x04\x05\x06\x07\x08\x09\x0a\x0b\x0c\x0d\x0e\x0f\x10\x11\x12\x13\x14\x15\x16\x17"
aes_set_key(ctx, key192, 192)
aes_encrypt(ctx, plaintext, encrypted); @test encrypted==b"\xdd\xa9\x7c\xa4\x86\x4c\xdf\xe0\x6e\xaf\x70\xa0\xec\x0d\x71\x91"
aes_decrypt(ctx, encrypted, decrypted); @test decrypted==b"\x00\x11\x22\x33\x44\x55\x66\x77\x88\x99\xaa\xbb\xcc\xdd\xee\xff"
println("- Passed")

print(" C.3 AES-256 ")
key256 = b"\x00\x01\x02\x03\x04\x05\x06\x07\x08\x09\x0a\x0b\x0c\x0d\x0e\x0f\x10\x11\x12\x13\x14\x15\x16\x17\x18\x19\x1a\x1b\x1c\x1d\x1e\x1f"
aes_set_key(ctx, key256, 256)
aes_encrypt(ctx, plaintext, encrypted); @test encrypted==b"\x8e\xa2\xb7\xca\x51\x67\x45\xbf\xea\xfc\x49\x90\x4b\x49\x60\x89"
aes_decrypt(ctx, encrypted, decrypted); @test decrypted==b"\x00\x11\x22\x33\x44\x55\x66\x77\x88\x99\xaa\xbb\xcc\xdd\xee\xff"
println("- Passed")

#==============================================================================
For fun you could scramble your data the other way: first with aes_decrypt,
then aes_encrypt will get it back when provided with the same key.
==============================================================================#

println("NIST - Spec. Publ. 800-38A 2001 ED test cases:")
TxtBlock1 = b"\x6b\xc1\xbe\xe2\x2e\x40\x9f\x96\xe9\x3d\x7e\x11\x73\x93\x17\x2a"
TxtBlock2 = b"\xae\x2d\x8a\x57\x1e\x03\xac\x9c\x9e\xb7\x6f\xac\x45\xaf\x8e\x51"
TxtBlock3 = b"\x30\xc8\x1c\x46\xa3\x5c\xe4\x11\xe5\xfb\xc1\x19\x1a\x0a\x52\xef"
TxtBlock4 = b"\xf6\x9f\x24\x45\xdf\x4f\x9b\x17\xad\x2b\x41\x7b\xe6\x6c\x37\x10"

print(" F.1.1 ECB-AES128.Encrypt; F.1.2 ECB-AES128.Decrypt ")
Key128 = b"\x2b\x7e\x15\x16\x28\xae\xd2\xa6\xab\xf7\x15\x88\x09\xcf\x4f\x3c"
aes_set_key(ctx, Key128, 128)
aes_encrypt(ctx, TxtBlock1, encrypted); @test encrypted==b"\x3a\xd7\x7b\xb4\x0d\x7a\x36\x60\xa8\x9e\xca\xf3\x24\x66\xef\x97"
aes_decrypt(ctx, encrypted, decrypted); @test decrypted==b"\x6b\xc1\xbe\xe2\x2e\x40\x9f\x96\xe9\x3d\x7e\x11\x73\x93\x17\x2a"

aes_encrypt(ctx, TxtBlock2, encrypted); @test encrypted==b"\xf5\xd3\xd5\x85\x03\xb9\x69\x9d\xe7\x85\x89\x5a\x96\xfd\xba\xaf"
aes_decrypt(ctx, encrypted, decrypted); @test decrypted==b"\xae\x2d\x8a\x57\x1e\x03\xac\x9c\x9e\xb7\x6f\xac\x45\xaf\x8e\x51"

aes_encrypt(ctx, TxtBlock3, encrypted); @test encrypted==b"\x43\xb1\xcd\x7f\x59\x8e\xce\x23\x88\x1b\x00\xe3\xed\x03\x06\x88"
aes_decrypt(ctx, encrypted, decrypted); @test decrypted==b"\x30\xc8\x1c\x46\xa3\x5c\xe4\x11\xe5\xfb\xc1\x19\x1a\x0a\x52\xef"

aes_encrypt(ctx, TxtBlock4, encrypted); @test encrypted==b"\x7b\x0c\x78\x5e\x27\xe8\xad\x3f\x82\x23\x20\x71\x04\x72\x5d\xd4"
aes_decrypt(ctx, encrypted, decrypted); @test decrypted==b"\xf6\x9f\x24\x45\xdf\x4f\x9b\x17\xad\x2b\x41\x7b\xe6\x6c\x37\x10"
println("- Passed")

print(" F.1.3 ECB-AES192.Encrypt; F.1.4 ECB-AES192.Decrypt ")
Key192 = b"\x8e\x73\xb0\xf7\xda\x0e\x64\x52\xc8\x10\xf3\x2b\x80\x90\x79\xe5\x62\xf8\xea\xd2\x52\x2c\x6b\x7b"
aes_set_key(ctx, Key192, 192)
aes_encrypt(ctx, TxtBlock1, encrypted); @test encrypted==b"\xbd\x33\x4f\x1d\x6e\x45\xf2\x5f\xf7\x12\xa2\x14\x57\x1f\xa5\xcc"
aes_decrypt(ctx, encrypted, decrypted); @test decrypted==b"\x6b\xc1\xbe\xe2\x2e\x40\x9f\x96\xe9\x3d\x7e\x11\x73\x93\x17\x2a"

aes_encrypt(ctx, TxtBlock2, encrypted); @test encrypted==b"\x97\x41\x04\x84\x6d\x0a\xd3\xad\x77\x34\xec\xb3\xec\xee\x4e\xef"
aes_decrypt(ctx, encrypted, decrypted); @test decrypted==b"\xae\x2d\x8a\x57\x1e\x03\xac\x9c\x9e\xb7\x6f\xac\x45\xaf\x8e\x51"

aes_encrypt(ctx, TxtBlock3, encrypted); @test encrypted==b"\xef\x7a\xfd\x22\x70\xe2\xe6\x0a\xdc\xe0\xba\x2f\xac\xe6\x44\x4e"
aes_decrypt(ctx, encrypted, decrypted); @test decrypted==b"\x30\xc8\x1c\x46\xa3\x5c\xe4\x11\xe5\xfb\xc1\x19\x1a\x0a\x52\xef"

aes_encrypt(ctx, TxtBlock4, encrypted); @test encrypted==b"\x9a\x4b\x41\xba\x73\x8d\x6c\x72\xfb\x16\x69\x16\x03\xc1\x8e\x0e"
aes_decrypt(ctx, encrypted, decrypted); @test decrypted==b"\xf6\x9f\x24\x45\xdf\x4f\x9b\x17\xad\x2b\x41\x7b\xe6\x6c\x37\x10"
println("- Passed")

print(" F.1.5 ECB-AES256.Encrypt; F.1.6 ECB-AES256.Decrypt ")
Key256 = b"\x60\x3d\xeb\x10\x15\xca\x71\xbe\x2b\x73\xae\xf0\x85\x7d\x77\x81\x1f\x35\x2c\x07\x3b\x61\x08\xd7\x2d\x98\x10\xa3\x09\x14\xdf\xf4"
aes_set_key(ctx, Key256, 256)
aes_encrypt(ctx, TxtBlock1, encrypted); @test encrypted==b"\xf3\xee\xd1\xbd\xb5\xd2\xa0\x3c\x06\x4b\x5a\x7e\x3d\xb1\x81\xf8"
aes_decrypt(ctx, encrypted, decrypted); @test decrypted==b"\x6b\xc1\xbe\xe2\x2e\x40\x9f\x96\xe9\x3d\x7e\x11\x73\x93\x17\x2a"

aes_encrypt(ctx, TxtBlock2, encrypted); @test encrypted==b"\x59\x1c\xcb\x10\xd4\x10\xed\x26\xdc\x5b\xa7\x4a\x31\x36\x28\x70"
aes_decrypt(ctx, encrypted, decrypted); @test decrypted==b"\xae\x2d\x8a\x57\x1e\x03\xac\x9c\x9e\xb7\x6f\xac\x45\xaf\x8e\x51"

aes_encrypt(ctx, TxtBlock3, encrypted); @test encrypted==b"\xb6\xed\x21\xb9\x9c\xa6\xf4\xf9\xf1\x53\xe7\xb1\xbe\xaf\xed\x1d"
aes_decrypt(ctx, encrypted, decrypted); @test decrypted==b"\x30\xc8\x1c\x46\xa3\x5c\xe4\x11\xe5\xfb\xc1\x19\x1a\x0a\x52\xef"

aes_encrypt(ctx, TxtBlock4, encrypted); @test encrypted==b"\x23\x30\x4b\x7a\x39\xf9\xf3\xff\x06\x7d\x8d\x8f\x9e\x24\xec\xc7"
aes_decrypt(ctx, encrypted, decrypted); @test decrypted==b"\xf6\x9f\x24\x45\xdf\x4f\x9b\x17\xad\x2b\x41\x7b\xe6\x6c\x37\x10"
println("- Passed")

#==============================================================================
File encryption/decryption utility examples:

file_encrypt("original_file.txt", "encrypted_128.aes", 128, "up_to_16char_pwd")
file_decrypt("encrypted_128.aes", "original_file.txt", 128, "up_to_16char_pwd")

file_encrypt("original_file.txt", "encrypted_192.aes", 192, "upto24character_password")
file_decrypt("encrypted_192.aes", "original_file.txt", 192, "upto24character_password")

file_encrypt("original_file.txt", "encrypted_256.aes", 256, "longest_upto32character_password")
file_decrypt("encrypted_256.aes", "original_file.txt", 256, "longest_upto32character_password")

==============================================================================#


text:
function addkey!(membernames, nam)
    if !haskey(membernames, nam)
        membernames[nam] = gensym()
    end
    membernames[nam]
end

onearg(e::Expr, f) = e.head == :call && length(e.args) == 2 && e.args[1] == f
onearg(e, f) = false

mapexpr(f, e) = Expr(e.head, map(f, e.args)...)

replace_syms!(x, membernames) = x
replace_syms!(q::QuoteNode, membernames) =
    replace_syms!(Meta.quot(q.value), membernames)
replace_syms!(e::Expr, membernames) =
    if onearg(e, :^)
        e.args[2]
    elseif onearg(e, :_I_)
        @warn "_I_() for escaping variables is deprecated, use cols() instead"
        addkey!(membernames, :($(e.args[2])))
    elseif onearg(e, :cols)
        addkey!(membernames, :($(e.args[2])))
    elseif e.head == :quote
        addkey!(membernames, Meta.quot(e.args[1]) )
    elseif e.head == :.
        replace_dotted!(e, membernames)
    else
        mapexpr(x -> replace_syms!(x, membernames), e)
    end

is_simple_non_broadcast_call(x) = false
function is_simple_non_broadcast_call(expr::Expr)
    expr.head == :call &&
        length(expr.args) >= 2 &&
        expr.args[1] isa Symbol &&
        all(x -> x isa QuoteNode || onearg(x, :cols), expr.args[2:end])
end

is_simple_broadcast_call(x) = false
function is_simple_broadcast_call(expr::Expr)
    expr.head == :. &&
        length(expr.args) == 2 &&
        expr.args[1] isa Symbol &&
        expr.args[2] isa Expr &&
        expr.args[2].head == :tuple &&
        all(x -> x isa QuoteNode || onearg(x, :cols), expr.args[2].args)
end

function args_to_selectors(v)
    t = map(v) do arg
        if arg isa QuoteNode
            arg
        elseif onearg(arg, :cols)
            arg.args[2]
        else
            throw(ArgumentError("This path should not be reached, arg: $(arg)"))
        end
    end

    :(DataFramesMeta.make_source_concrete($(Expr(:vect, t...))))
end

is_macro_head(ex, name) = false
is_macro_head(ex::Expr, name) = ex.head == :macrocall && ex.args[1] == Symbol(name)

extract_macro_flags(ex, exprflags = (;Symbol("@byrow") => Ref(false),)) = (ex, exprflags)
function extract_macro_flags(ex::Expr, exprflags = (;Symbol("@byrow") => Ref(false),))
    if ex.head == :macrocall
        macroname = ex.args[1]
        if macroname in keys(exprflags)
            exprflag = exprflags[macroname]
            if exprflag[] == true
                throw(ArgumentError("Redundant flag $macroname used."))
            end
            exprflag[] = true
            return extract_macro_flags(MacroTools.unblock(ex.args[3]), exprflags)
        else
            return (ex, exprflags)
        end
    end

    return (ex, exprflags)
end
"""
    get_source_fun(function_expr; wrap_byrow::Bool=false)

Given an expression that may contain `QuoteNode`s (`:x`)
and items wrapped in `cols`, return a function
that is equivalent to that expression where the
`QuoteNode`s and `cols` items are the inputs
to the function.

For fast compilation `get_source_fun` returns
the name of a called function where possible.

* `f(:x, :y)` will return `f`
* `f.(:x, :y)` will return `ByRow(f)`
* `:x .+ :y` will return `.+`

`get_source_fun` also returns an expression
representing the vector of inputs that will be
used as the `src` in the `src => fun => dest`
call later on.

If `wrap_byrow=true` then the function gets wrapped
in `ByRow`. If the expression begins with `@byrow`,
then `get_source_fun` is recurively called on the
expression that `@byrow` acts on, with `wrap_byrow=true`.

### Examples

julia> using MacroTools

julia> ex = :(:x + :y);

julia> DataFramesMeta.get_source_fun(ex)
(:(DataFramesMeta.make_source_concrete([:x, :y])), :+)

julia> ex = quote
           :x .+ 1 .* :y
       end |> MacroTools.prettify

julia> src, fun = DataFramesMeta.get_source_fun(ex);

julia> MacroTools.prettify(fun)
:((mammoth, goat)->mammoth .+ 1 .* goat)

julia> ex = :(@byrow :x * :y);

julia> src, fun = DataFramesMeta.get_source_fun(ex);

julia> MacroTools.prettify(fun)
:(ByRow(*))
```

"""
function get_source_fun(function_expr; wrap_byrow::Bool=false)
    function_expr = MacroTools.unblock(function_expr)

    if is_simple_non_broadcast_call(function_expr)
        source = args_to_selectors(function_expr.args[2:end])
        fun_t = function_expr.args[1]

        # .+ to +
        if startswith(string(fun_t), '.')
            f_sym_without_dot = Symbol(chop(string(fun_t), head = 1, tail = 0))
            fun = :(DataFrames.ByRow($f_sym_without_dot))
        else
            fun = fun_t
        end
    elseif is_simple_broadcast_call(function_expr)
        # extract source symbols from quotenodes
        source = args_to_selectors(function_expr.args[2].args)
        fun_t = function_expr.args[1]
        fun = :(DataFrames.ByRow($fun_t))
    else
        membernames = Dict{Any, Symbol}()

        body = replace_syms!(function_expr, membernames)

        source = :(DataFramesMeta.make_source_concrete($(Expr(:vect, keys(membernames)...))))
        inputargs = Expr(:tuple, values(membernames)...)

        fun = quote
            $inputargs -> begin
                $body
            end
        end
    end

    if wrap_byrow
        fun = :(ByRow($fun))
    end

    return source, fun
end

# `nolhs` needs to be `true` when we have syntax of the form
# `@combine(gd, fun(:x, :y))` where `fun` returns a `table` object.
# We don't create the "new name" pair because new names are
# given by the table.
# We need wrap_byrow as a keyword argument here in case someone
# uses `@transform df @byrow begin ... end`, which we
# deal with outside of this function.
function fun_to_vec(ex::Expr;
                    gensym_names::Bool=false,
                    outer_flags::Union{NamedTuple, Nothing}=nothing,
                    no_dest::Bool=false)
    # classify the type of expression
    # :x # handled via dispatch
    # cols(:x) # handled as though above
    # f(:x) # requires pass_as_is, for `@with` and `@subset` in future
    # y = :x # :x is a QuoteNode
    # y = cols(:x) # use cols on RHS
    # cols(:y) = :x # RHS in :block
    # cols(:y) = cols(:x)
    # y = f(:x) # re-write as simple call
    # y = f(cols(:x)) # re-write as simple call, use cols
    # y = :x + 1 # re-write as complicated call
    # y = cols(:x) + 1 # re-write as complicated call, with cols
    # cols(:y) = f(:x) # re-write as simple call, but RHS is :block
    # cols(:y) = f(cols(:x)) # re-write as simple call, RHS is block, use cols
    # cols(y) = :x + 1 # re-write as complicated col, but RHS is :block
    # cols(:y) = cols(:x) + 1 # re-write as complicated call, RHS is block, use cols
    # `@byrow` before any of the above
    ex, inner_flags = extract_macro_flags(MacroTools.unblock(ex))

    # Use tuple syntax in future when we add more flags
    inner_wrap_byrow = inner_flags[Symbol("@byrow")][]
    outer_wrap_byrow = outer_flags === nothing ? false : outer_flags[Symbol("@byrow")][]

    if inner_wrap_byrow && outer_wrap_byrow
        throw(ArgumentError("Redundant @byrow calls."))
    else
        wrap_byrow = inner_wrap_byrow || outer_wrap_byrow
    end

    if gensym_names
        ex = Expr(:kw, gensym(), ex)
    end

    # :x
    # handled below via dispatch on ::QuoteNode

    # cols(:x)
    if onearg(ex, :cols)
        return ex.args[2]
    end

    if no_dest
        source, fun = get_source_fun(ex, wrap_byrow = wrap_byrow)
        return quote
            $source => $fun
        end
    end

    @assert ex.head == :kw || ex.head == :(=)
    lhs = ex.args[1]
    rhs = MacroTools.unblock(ex.args[2])

    if is_macro_head(rhs, "@byrow")
        s = "In keyword argument inputs, `@byrow` must be on the left hand side. " *
        "Did you write `y = @byrow f(:x)` instead of `@byrow y = f(:x)`?"
        throw(ArgumentError(s))
    end

    # y = ...
    if lhs isa Symbol
        msg = "Using an un-quoted Symbol on the LHS is deprecated. " *
              "Write $(QuoteNode(lhs)) = ... instead."
        Base.depwarn(msg, "")
        lhs = QuoteNode(lhs)
    end

    # :y = :x
    if lhs isa QuoteNode && rhs isa QuoteNode
        source = rhs
        dest = lhs

        return quote
            $source => $dest
        end
    end

    # :y = cols(:x)
    if lhs isa QuoteNode && onearg(rhs, :cols)
        source = rhs.args[2]
        dest = lhs

        return quote
            $source => $dest
        end
    end

    # cols(:y) = :x
    if onearg(lhs, :cols) && rhs isa QuoteNode
        source = rhs
        dest = lhs.args[2]

        return quote
            $source => $dest
        end
    end

    # cols(:y) = cols(:x)
    if onearg(lhs, :cols) && onearg(rhs, :cols)
        source = rhs.args[2]
        dest = lhs.args[2]
        return quote
            $source => $dest
        end
    end

    # :y = f(:x)
    # :y = f(cols(:x))
    # :y = :x + 1
    # :y = cols(:x) + 1
    source, fun = get_source_fun(rhs; wrap_byrow = wrap_byrow)
    if lhs isa QuoteNode
        dest = lhs
        return quote
            $source => $fun => $dest
        end
    end

    # cols(:y) = f(:x)
    if onearg(lhs, :cols)
        dest = lhs.args[2]

        return quote
            $source => $fun => $dest
        end
    end

    throw(ArgumentError("This path should not be reached"))
end
fun_to_vec(ex::QuoteNode;
           no_dest::Bool=false,
           gensym_names::Bool=false,
           outer_flags::Union{NamedTuple, Nothing}=nothing) = ex

function make_source_concrete(x::AbstractVector)
    if isempty(x) || isconcretetype(eltype(x))
        return x
    elseif all(t -> t isa Union{AbstractString, Symbol}, x)
        return Symbol.(x)
    else
        throw(ArgumentError("Column references must be either all the same " *
                            "type or a a combination of `Symbol`s and strings"))
    end
end

protect_replace_syms!(e, membernames) = e
function protect_replace_syms!(e::Expr, membernames)
    if e.head == :quote
        e
    else
        replace_syms!(e, membernames)
    end
end

function replace_dotted!(e, membernames)
    x_new = replace_syms!(e.args[1], membernames)
    y_new = protect_replace_syms!(e.args[2], membernames)
    Expr(:., x_new, y_new)
end

function create_args_vector(args...)
    create_args_vector(Expr(:block, args...))
end

"""
   create_args_vector(arg) -> vec, outer_flags

Given an expression return a vector of operations
and a `NamedTuple` of the macro-flags that appear
in the expression.

If a `:block` expression, return the `args` of
the block as an array. If a simple expression,
wrap the expression in a one-element vector.
"""
function create_args_vector(arg)
    arg, outer_flags = extract_macro_flags(MacroTools.unblock(arg))

    if arg isa Expr && arg.head == :block
        x = MacroTools.rmlines(arg).args
    else
        x = Any[arg]
    end

    return x, outer_flags
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    SineHoleVariogram(range=r, sill=s, nugget=n)
    SineHoleVariogram(ball; sill=s, nugget=n)

A sine hole variogram with range `r`, sill `s` and nugget `n`.
Optionally, use a custom metric `ball`.
"""
struct SineHoleVariogram{T,B} <: Variogram
  sill::T
  nugget::T
  ball::B
end

SineHoleVariogram(ball; sill=1.0, nugget=0.0) =
  SineHoleVariogram(sill, nugget, ball)

SineHoleVariogram(; range=1.0, sill=1.0, nugget=0.0) =
  SineHoleVariogram(sill, nugget, MetricBall(range))

function (γ::SineHoleVariogram)(h::T) where {T}
  r = radius(γ.ball)
  s = γ.sill
  n = γ.nugget

  # shift lag by machine precision to
  # avoid explosion at the origin
  h′ = h + eps(T)
  c  = T(π)

  (s - n) * (1 - sin(c*h′/r)/(c*h′/r)) + (h′ > 0) * n
end

isstationary(::Type{<:SineHoleVariogram}) = true


text:
using DocExample
using Test

@testset "DocExample.jl" begin
    # Write your tests here.
end


text:
import Libdl

# the following functions are adapted from the julia-config.jl script and are used to
# determine parameters for the Mex file build process

isDebug() = ccall(:jl_is_debugbuild, Cint, ()) != 0

function shell_escape(str)
    str = replace(str, "'" => "'\''")
    return "\"$str\""
end

function matlab_escape(str)
    str = replace(str, "'" => "''")
    return "'$str'"
end

function libDir()
    return if isDebug()
        dirname(abspath(Libdl.dlpath("libjulia-debug")))
    else
        dirname(abspath(Libdl.dlpath("libjulia")))
    end
end

private_libDir() = abspath(Sys.BINDIR, Base.PRIVATE_LIBDIR)

function includeDir()
    return abspath(Sys.BINDIR, Base.INCLUDEDIR, "julia")
end

function cflags()
    fl = "-I $(shell_escape(includeDir()))"
    if Sys.isunix()
        fl = fl * " -fPIC"
    end
    return fl
end

function ldflags()
    fl = "-L$(shell_escape(libDir()))"
    if Sys.isunix()
        fl = fl * " -Wl,-rpath $(shell_escape(libDir()))"
    end
    return fl
end

function ldlibs()
    libname = if isDebug()
        "julia-debug"
    else
        "julia"
    end
    if Sys.isunix()
        return "-l$libname -ldl"
    else
        return "\'$(normpath(joinpath(libDir(), "..", "lib", "lib$libname.dll.a")))\'"
    end
end

# the following two functions are taken from the MATLAB.jl build script and are used to
# locate MATLAB and its relevant libraries and commands

function find_matlab_root()
    # Determine MATLAB library path and provide facilities to load libraries with this path
    matlab_root = get(ENV, "MATLAB_ROOT",
                        get(ENV, "MATLAB_HOME", nothing))
    if isnothing(matlab_root)
        matlab_exe = Sys.which("matlab")
        if !isnothing(matlab_exe)
            matlab_root = dirname(dirname(matlab_exe))
        else
            if Sys.isapple()
                default_dir = "/Applications"
                if isdir(default_dir)
                    dirs = readdir(default_dir)
                    filter!(app -> occursin(r"^MATLAB_R[0-9]+[ab]\.app$", app), dirs)
                    if !isempty(dirs)
                        matlab_root = joinpath(default_dir, maximum(dirs))
                    end
                end
            elseif Sys.iswindows()
                default_dir = Sys.WORD_SIZE == 32 ? "C:\\Program Files (x86)\\MATLAB" : "C:\\Program Files\\MATLAB"
                if isdir(default_dir)
                    dirs = readdir(default_dir)
                    filter!(dir -> occursin(r"^R[0-9]+[ab]$", dir), dirs)
                    if !isempty(dirs)
                        matlab_root = joinpath(default_dir, maximum(dirs))
                    end
                end
            end
        end
    end
    !isnothing(matlab_root) && isdir(matlab_root) && @info("Detected MATLAB root folder at \"$matlab_root\"")
    return matlab_root
end

function find_matlab_cmd(matlab_root)
    if Sys.iswindows()
        matlab_cmd = joinpath(matlab_root, "bin", (Sys.WORD_SIZE == 32 ? "win32" : "win64"), "matlab.exe")
        isfile(matlab_cmd) && @info("Detected MATLAB executable at \"$matlab_cmd\"")
    else
        matlab_exe = joinpath(matlab_root, "bin", "matlab")
        isfile(matlab_exe) && @info("Detected MATLAB executable at \"$matlab_exe\"")
        matlab_cmd = "$(Base.shell_escape(matlab_exe))"
    end
    return matlab_cmd
end

function mex_extension()
    if Sys.islinux()
        ext = ".mexa64"
    elseif Sys.isapple()
        ext = ".mexmaci64"
    elseif Sys.iswindows()
        ext = ".mexw"*string(Sys.WORD_SIZE)
    end
    return ext
end

is_ci() = lowercase(get(ENV, "CI", "false")) == "true"

# find matlab root
matlab_root = find_matlab_root()

if !isnothing(matlab_root)

    # get matlab command
    matlab_cmd = find_matlab_cmd(matlab_root)

    # get build parameters
    is_debug = isDebug()
    julia_bin = is_debug ?
        joinpath(unsafe_string(Base.JLOptions().julia_bindir), "julia-debug") :
        joinpath(unsafe_string(Base.JLOptions().julia_bindir), "julia")
    julia_home = unsafe_string(Base.JLOptions().julia_bindir)
    sys_image = unsafe_string(Base.JLOptions().image_file)
    lib_base = is_debug ? "julia-debug" : "julia"
    lib_path = Libdl.dlpath("lib$lib_base")
    lib_dir = Sys.iswindows() ? joinpath(dirname(julia_home), "lib") : dirname(lib_path)
    inc_dir = includeDir()
    build_cflags = cflags()
    build_ldflags = ldflags()
    build_ldlibs = ldlibs()
    build_src = abspath("mexjulia.cpp")
    outdir = normpath(joinpath(pwd(),"..","mexjulia"))
    mex_cmd = "mex -v -largeArrayDims -outdir \'$outdir\' LDFLAGS=\'$(build_ldflags) \$LDFLAGS\' CFLAGS=\'$(build_cflags) \$CFLAGS\' \'$(build_src)\' $(build_ldlibs)"

    # generate Mex-file build script
    dict_file = joinpath(outdir, "jldict.mat")
    build_file = joinpath(pwd(), "build.m")

    open(build_file, "w") do io
        println(io,
            """
            % This file is automatically generated, do not edit.

            % Save build parameters to a .mat file
            is_debug = $is_debug;
            julia_bin = $(matlab_escape(julia_bin));
            julia_home = $(matlab_escape(julia_home));
            sys_image = $(matlab_escape(sys_image));
            lib_base = $(matlab_escape(lib_base));
            lib_path = $(matlab_escape(lib_path));
            lib_dir = $(matlab_escape(lib_dir));
            inc_dir = $(matlab_escape(inc_dir));
            build_cflags = $(matlab_escape(build_cflags));
            build_ldflags = $(matlab_escape(build_ldflags));
            build_ldlibs = $(matlab_escape(build_ldlibs));
            build_src = $(matlab_escape(build_src));
            mex_cmd = $(matlab_escape(mex_cmd));
            save($(matlab_escape(dict_file)), "is_debug", "julia_bin", "julia_home",...
                "sys_image", "lib_base", "lib_path", "lib_dir", "inc_dir", "build_cflags",...
                "build_ldflags", "build_ldlibs", "build_src", "mex_cmd");

            % Run Mex Command
            $mex_cmd;

            """
        )

        if !is_ci()
            println(io,
                """
                % Check if the `mexjulia` directory is already on the path
                path_dirs = regexp(path, pathsep, 'split');
                if ispc
                    on_path = any(strcmpi($(matlab_escape(outdir)), path_dirs));
                else
                    on_path = any(strcmp($(matlab_escape(outdir)), path_dirs));
                end

                % Add the `mexjulia` directory to the path and save the path (if necessary)
                if ~on_path
                    fprintf('%s is not on the MATLAB path. Adding it and saving...\\n\', $(matlab_escape(outdir)));
                    path($(matlab_escape(outdir)), path);
                    savepath;
                end
                """
            )
        end
    end

    # We have to run the build script separately for CI due to licensing issues
    if !is_ci()

        # run the Mex-file build script
        run(`$matlab_cmd -nodesktop -nosplash -r "run('$build_file');exit"`)

        # check that the build information has been saved in the mexjulia directory
        @assert isfile(joinpath(outdir, "jldict.mat"))

        # check that the compiled mexjulia file has been saved in the mexjulia directory
        @assert isfile(joinpath(outdir, "mexjulia" * mex_extension()))

    end

elseif get(ENV, "JULIA_REGISTRYCI_AUTOMERGE", nothing) == "true"
    # We need to be able to install and load this package without error for
    # Julia's registry AutoMerge to work, so we just skip the mex file build process.
else
    error("MATLAB cannot be found. Set the \"MATLAB_ROOT\" environment variable to the MATLAB root directory and re-run Pkg.build(\"Mex\").")
end

max_stars_repo_name: taylormcd/mexjulia

max_issues_repo_name: taylormcd/mexjulia

max_forks_repo_name: taylormcd/mexjulia


text:
function replace_low!(A::Array{Int32,3}, trs::Int32; nodata = missing)
    inds = CartesianIndices(A)
    n = length(inds)
    @inbounds for ind = inds
        if !ismissing(A[ind]) && (A[ind] < trs)
            A[ind] = nodata;
        end
    end
end

function mapvalues(A::AbstractArray{T,2}, levs_old, levs_new) where T <: Real
    A2 = deepcopy(A);
    inds = CartesianIndices(A)
    n = length(inds)
    nlev = length(levs_old)
    @inbounds for ind = inds
        if !ismissing(A[ind])
            for i in 1:nlev
                if A[ind] == levs_old[i]
                    A2[ind] = levs_new[i]
                    break
                end
            end
        end
    end
    A2
end

# change ArcGIS direction into taudem direction
function flowdir_gis2tau(A::AbstractArray{T,2}) where T <: Real
    dir_gis = Int32.([1, 2, 4, 8, 16, 32, 64, 128]);
    dir_taudem = Int32.([1, 8, 7, 6, 5, 4, 3, 2]);
    mapvalues(A, dir_gis, dir_taudem)
end

function flowdir_tau2gis(A::AbstractArray{T,2}) where T <: Real
    dir_gis = Int32.([1, 2, 4, 8, 16, 32, 64, 128]);
    dir_taudem = Int32.([1, 8, 7, 6, 5, 4, 3, 2]);
    mapvalues(A, dir_taudem, dir_gis)
end

export flowdir_gis2tau, flowdir_tau2gis, mapvalues, replace_low!


text:
module CTD

import DataFrames: DataFrame, DataFrameRow
import TimeZones: ZonedDateTime
import Dates: Second
import ..ARLToolkit: dts

function read(filename; mindepth=1.0, minmaxdepth=5.0)
  df = DataFrame(
    time = ZonedDateTime[],
    maxdepth = Float64[],
    lno = UnitRange{Int}[],
    filename = String[])
  t0 = nothing
  Δt = 2
  start = 0
  n0 = 0
  maxp = 0.0
  for (lno, line) ∈ enumerate(eachline(filename))
    m = match(r"^Date (.*) Time (.*) GMT Interval: (\d+)sec", line)
    if m !== nothing
      t0 = dts("20" * m[1] * " " * m[2]; format="yyyy-mm-dd HH-MM-SS")
      Δt = parse(Int, m[3])
      n0 = 0
      continue
    end
    m = match(r"^N(\d+) C(.+) T(.+) P(.+) c(.+)$", line)
    if m !== nothing
      n = parse(Int, m[1])
      n0 == 0 && (n0 = n)
      C = parse(Float64, m[2])
      p = parse(Float64, m[4])
      maxp = max(p, maxp)
      if p > mindepth && C > 1.0
        if start == 0
          start = lno
          maxp = 0.0
        end
      else
        if start > 0 && maxp > minmaxdepth
          push!(df, (t0 + Second(Δt * (n-n0)), maxp, start:(lno-1), filename))
          start = 0
        end
      end
      continue
    end
  end
  df
end

read(filenames::AbstractVector) = sort!(vcat(read.(f)...), :time)

function read(row::DataFrameRow)
  df = DataFrame(
    n = Int[],
    conductivity = Float64[],
    temperature = Float64[],
    depth = Float64[],
    soundspeed = Float64[])
  for (lno, line) ∈ enumerate(eachline(row.filename))
    lno ∈ row.lno || continue
    m = match(r"^N(\d+) C(.+) T(.+) P(.+) c(.+)$", line)
    if m !== nothing
      n = parse(Int, m[1])
      C = parse(Float64, m[2])
      T = parse(Float64, m[3])
      D = parse(Float64, m[4])
      c = parse(Float64, m[5])
      push!(df, (n, C, T, D, c))
    end
  end
  df
end

read(df::DataFrame, n) = read(df[n,:])

end # module


text:
#
# SoapySDR Modules API
#
# https://github.com/pothosware/SoapySDR/blob/1cf5a539a21414ff509ff7d0eedfc5fa8edb90c6/include/SoapySDR/Modules.h
#
# NOTE: This is an API for manually loading Soapy Modules. It is recommended that developers
# package via BinaryBuilder+Yggdrasil for Julia.

""" Query the root installation path"""
function SoapySDR_getRootPath()
    #SOAPY_SDR_API const char *SoapySDR_getRootPath(void);
    @check_error ccall((:SoapySDR_getRootPath, lib), Cstring, ())
end

"""
The list of paths automatically searched by loadModules().

param [out] length the number of elements in the result.
return a list of automatically searched file paths
"""
function SoapySDR_listSearchPaths()
    #SOAPY_SDR_API char **SoapySDR_listSearchPaths(size_t *length);
    len = Ref{Csize_t}()
    ptr = @check_error ccall((:SoapySDR_listSearchPaths, lib), Ptr{Cstring}, (Ref{Csize_t},), len)
    (ptr, len[])
end

"""
List all modules found in default path.
The result is an array of strings owned by the caller.

param [out] length the number of elements in the result.
return a list of file paths to loadable modules
"""
function SoapySDR_listModules()
    #SOAPY_SDR_API char **SoapySDR_listModules(size_t *length);
    len = Ref{Csize_t}()
    ptr = @check_error ccall((:SoapySDR_listModules, lib), Ptr{Cstring}, (Ref{Csize_t},), len)
    (ptr, len[])
end

"""
List all modules found in the given path.
The result is an array of strings owned by the caller.

param path a directory on the system
param [out] length the number of elements in the result.
return a list of file paths to loadable modules
"""
function SoapySDR_listModulesPath(path)
    #SOAPY_SDR_API char **SoapySDR_listModulesPath(const char *path, size_t *length);
    len = Ref{Csize_t}()
    ptr = @check_error ccall((:SoapySDR_listModulesPath, lib), Ptr{Cstring}, (Cstring, Ref{Csize_t}), path, len)
    (ptr, len[])
end

"""
Load a single module given its file system path.
The caller must free the result error string.

param path the path to a specific module file]
return an error message, empty on success
"""
function SoapySDR_loadModule(path)
    #SOAPY_SDR_API char *SoapySDR_loadModule(const char *path);
    err = @check_error ccall((:SoapySDR_loadModule, lib), SoapySDRKwargs, (Cstring,), path)
    err
end

"""
List all registration loader errors for a given module path.
The resulting dictionary contains all registry entry names
provided by the specified module. The value of each entry
is an error message string or empty on successful load.

param path the path to a specific module file
return a dictionary of registry names to error messages
"""
function SoapySDR_listLoaderResult(path)
    #SOAPY_SDR_API SoapySDRKwargs SoapySDR_getLoaderResult(const char *path);
    kwargs = @check_error ccall((:SoapySDR_listLoaderResult, lib), SoapySDRKwargs, (Cstring,), path)
    kwargs
end

"""
Get a version string for the specified module.
Modules may optionally provide version strings.

param path the path to a specific module file
return a version string or empty if no version provided
"""
function SoapySDR_getModuleVersion(path)
    #SOAPY_SDR_API char *SoapySDR_getModuleVersion(const char *path);
    ver = @check_error ccall((:SoapySDR_getModuleVersion, lib), Cstring, (Cstring,), path)
    ver
end

"""
Unload a module that was loaded with loadModule().
The caller must free the result error string.

param path the path to a specific module file
return an error message, empty on success
"""
function SoapySDR_unloadModule(path)
    #SOAPY_SDR_API char *SoapySDR_unloadModule(const char *path);
    err = @check_error ccall((:SoapySDR_loadModules, lib), Cstring, (Cstring,), path)
    err
end

"""
Load the support modules installed on this system.
This call will only actually perform the load once.
Subsequent calls are a NOP.
"""
function SoapySDR_loadModules()
    #SOAPY_SDR_API void SoapySDR_loadModules(void);
    @check_error ccall((:SoapySDR_loadModules, lib), Cvoid, ())
end

"""
Unload all currently loaded support modules.
"""
function SoapySDR_unloadModules()
    #SOAPY_SDR_API void SoapySDR_unloadModules(void);
    @check_error ccall((:SoapySDR_unloadModule, lib), Cvoid, ())
end


text:
# This file was generated, do not modify it. # hide
r = report(mtm)

res = r.plotting

ss = res.parameter_values[:,1]
cbt = res.parameter_values[:,2]

figure(figsize=(8,6))
tricontourf(ss, cbt, res.measurements)

xlabel("Sub sample", fontsize=14)
ylabel("Col sample by tree", fontsize=14)
xticks(fontsize=12)
yticks(fontsize=12)

savefig(joinpath(@OUTPUT, "EX-crabs-xgb-heatmap2.svg")) # hide


text:
"""
```
orientations = corner_orientations(img)
orientations = corner_orientations(img, corners)
orientations = corner_orientations(img, corners, kernel)
```

Returns the orientations of corner patches in an image. The orientation of a corner patch
is denoted by the orientation of the vector between intensity centroid and the corner.
The intensity centroid can be calculated as `C = (m01/m00, m10/m00)` where mpq is defined as -

    `mpq = (x^p)(y^q)I(y, x) for each p, q in the corner patch`

The kernel used for the patch can be given through the `kernel` argument. The default kernel used is a gaussian kernel of size `5x5`.
"""
function corner_orientations(img::AbstractArray{T, 2}, corners::Keypoints, kernel::Array{K, 2}) where {T<:Gray, K<:Real}
    h, w = size(kernel)
    pre_y = ceil(Int, (h - 1) / 2)
    pre_x = ceil(Int, (w - 1) / 2)
    post_y = floor(Int, (h - 1) / 2)
       post_x = floor(Int, (w - 1) / 2)
       img_padded = parent(padarray(img, Fill(0, (pre_y, pre_x), (post_y, post_x))))
    orientations = Float64[]
    for c in corners
        m10 = zero(T)
        m01 = zero(T)
        for i in 1:w
            col_sum = zero(T)
            for j in 1:h
                pixel = img_padded[c + CartesianIndex(j - 1, i - 1)] * kernel[j, i]
                m01 += pixel * (j - pre_y - 1)
                col_sum += pixel
            end
            m10 += col_sum * (i - pre_x - 1)
        end
        push!(orientations, atan(Float64(m01), Float64(m10)))
    end
    orientations
end

corner_orientations(img::AbstractArray{T, 2}, corners::Keypoints, kernel::Array{K, 2}) where {T, K<:Real} = corner_orientations(Gray.(img), corners, kernel)

function corner_orientations(img::AbstractArray, kernel::Array{K, 2}) where K<:Real
    corners = imcorner(img)
    corner_indexes = Keypoints(corners)
    corner_orientations(img, corner_indexes, kernel)
end

function corner_orientations(img::AbstractArray)
    corners = imcorner(img)
    corner_indexes = Keypoints(corners)
    kernel = Kernel.gaussian((2,2), (5, 5))
    kernel /= maxfinite(kernel)
    corner_orientations(img, corner_indexes, parent(kernel))
end

function corner_orientations(img::AbstractArray, corners::Keypoints)
    kernel = Kernel.gaussian((2,2), (5, 5))
    kernel /= maxfinite(kernel)
    corner_orientations(img, corners, parent(kernel))
end


text:
"""
    fit!(model, obs::AbstractArray)

Fit a model to data with least squares regression, using `curve_fit` from
LsqFit.jl. The passed in model should be initialised with sensible defaults,
these will be used as the initial parameters for the optimization.

Any (nested) `Real` fields on the struct are flattened to a parameter vector using
[Flatten.jl](http://github.com/rafaqz/Flatten.jl). Fields can be marked to ignore
using the `@flattenable` macro from [FieldMetadata.jl](http://github.com/rafaqz/FieldMetadata.jl).

## Arguments

- `model`: Any constructed [`RateModel`](@ref) or a `Tuple` of `RateModel`.
- `obs`: A `Vector` of `(val, rate)` tuples where `val` is the value of the
  x-axis variable (such as temperature), and `rate` is the growth rate observed.

## Returns

An updated `Model` containing the fitted parameter values.
"""
function fit!(model::Model, obs::AbstractArray)
    fit = curve_fit(first.(obs), last.(obs), collect(model[:val])) do xs, vals
        model[:val] = vals
        GrowthMaps.conditionalrate.(Ref(stripparams(model)), xs)
    end
    model[:val] = fit.param
    return model
end

"""
    manualfit!(model::Model, data::NamedTuple; obs=[], kw...) =

Returns the fitted model.

# Arguments

- `obs`: A `Vector` of `(val, rate)` tuples where `val` is the value of the
  x-axis variable (such as temperature), and `rate` is the growth rate observed.
- `data`: A `NamedTuple` of `AbstractVector`. The `NamedTuple` key/s must match
  the key required by the `Layer`/s.
- `obs`: Optional observations to scatter-plot over the curve. A `Vector` of `(val, rate)`
  tuples where `val` is the value of the x-axis variable (such as temperature), and `rate`
  is the growth rate observed.
- `kwargs`: passed to `plot`

# Example

```julia
p = 3e-01
ΔH_A = 3e4cal/mol
ΔH_L = -1e5cal/mol
ΔH_H = 3e5cal/mol
Thalf_L = 2e2K
Thalf_H = 3e2K
T_ref = K(25.0°C)
growthmodel = SchoolfieldIntrinsicGrowth(p, ΔH_A, ΔH_L, Thalf_L, ΔH_H, Thalf_H, T_ref)
model = Model(Layer(:surface_temp, K, growthmodel))
obs = []

tempdata=(surface_temp=(270.0:0.1:310.0)K,)
manualfit!(model, tempdata; obs=obs)

To use the interface in a desktop app, use Blink.jl:

```julia; eval=false
using Blink
w = Blink.Window()
body!(w, interface)
```
"""
function manualfit!(
    model::AbstractModel, data::NamedTuple{<:Any,<:Tuple{Vararg{<:AbstractVector}}};
    observations=[],
    throttle=0.1,
    kwargs...
)
    InteractModel(model; throttle=throttle, submodel=RateModel) do updated_model
        @show params(updated_model)
        ModelParameters.setparent!(model, updated_model)
        manualfit(stripparams(updated_model), (observations, data); kwargs...)
    end
end

function manualfit(model, (observations, data); kwargs...)
    predictions = combinelayers(model, data)
    p = plot(first(data), predictions; label="predicted", legend=false, kwargs...)
    scatter!(p, observations; label="observed")
    return p
end

"""
    mapfit!(model::Model, modelkwargs;
            occurrence=[],
            precomputed=nothing,
            throttle=0.1,
            window=(Band(1),),
            kwargs...
    )

Fit a model to the map.

# Arguments

- `occurence`: a `Vector` of occurence locations, as `(lon, lat)` tuples.
- `modelkwargs`: are passed to the `mapgrowth` with the model.
- `mapkwargs`: are passed to the `plot` function the plots the `GeoArray`
- `throttle`: the response time of Interact.jl sliders.
- `window`: selects a window of the output to plot. By default this is `(Band(1),)`,
  which just removes the `Band` dimension from a `GeoArray`, if it exists.
- `kwargs`: passed to `Plots.scatter!`

# Example

```julia
model = Model(wiltstress, coldstress, heatstress)
interface = mapfit!(model, modelkwargs;
    occurrence=occurrence,
    precomputed=precomputed,
    throttle=0.2,
    window=(Band(1),),
    markershape=:cross,
    markercolor=:lightblue,
    markeropacity=0.4
)
display(interface)
```

To use the interface in a desktop app, use Blink.jl:

```julia; eval=false
using Blink
w = Blink.Window()
body!(w, interface)
```
"""
function mapfit!(model::AbstractModel, modelkwargs;
    mapkwargs=(),
    occurrence=[],
    precomputed=nothing,
    throttle=0.1,
    scatterkwargs...
)
    title = "GrowthMaps: mapfit interface"
    InteractModel(model; throttle=throttle, title=title, submodel=RateModel) do updated_model
        ModelParameters.setparent!(model, updated_model)
        mapfit(
            stripparams(updated_model), (modelkwargs, occurrence, precomputed);
            scatterkwargs...
        )
    end
end

function mapfit(model, (modelkwargs, occurrence, precomputed);
    window=(Band(1),),
    levels=10,
    mapkwargs=(),
    markercolor=:white,
    markersize=2.0,
    clims=(0.0, 0.25),
    scatterkwargs...
)
    output = mapgrowth(model; modelkwargs...)
    output = isnothing(precomputed) ? output : output .+ precomputed
    windowed = output[window...]
    p = plot(windowed; legend=:none, levels=levels, clims=clims, mapkwargs...)
    scatter(; t...) = scatter!(
        p, occurrence;
        t..., markercolor=markercolor, markersize=markersize, scatterkwargs...
    )
    if hasdim(windowed, Ti())
        for t in 1:length(dims(windowed, Ti()))
            scatter(; subplot=t)
        end
    else
        scatter()
    end
    return p
end


text:
module PEPS

using ITensors, ITensorsGPU
using CUDA, TimerOutputs
import ITensors: data, store
import ITensors.NDTensors: data
using Random, Logging, LinearAlgebra, DelimitedFiles

const q_dict = Dict{Pair{Int, Symbol}, MPO}()

include("peps_util.jl")

export fPEPS, randomfPEPS, checkerboardfPEPS, cufPEPS, randomCufPEPS, buildEdgeEnvironment, buildNextEnvironment, buildLs, buildRs

end # module


text:
function choice_gradients(trace::VectorTrace{UnfoldType,T,U}, selection::Selection,
                        retval_grad) where {T,U}
    error("Not implemented")
end

function accumulate_param_gradients!(trace::VectorTrace{UnfoldType,T,U}, retval_grad) where {T,U}
    args = get_args(trace)
    (len, init_state, params) = unpack_args(args)

    kernel_has_grads = has_argument_grads(trace.gen_fn.kernel)
    if kernel_has_grads[1]
        error("Cannot differentiate with respect to index in unfold")
    end
    if kernel_has_grads[2]
        # backpropagation through the state not yet implemented
        error("Not implemented")
    end
    params_has_grad = kernel_has_grads[3:end]

    params_grad = Vector(undef, length(params_has_grad))
    for (i, has_grad) in enumerate(params_has_grad)
        if has_grad
            params_grad[i] = Vector(undef, len)
        else
            params_grad[i] = nothing
        end
    end

    for key=1:len
        subtrace = trace.subtraces[key]
        kernel_retval_grad = (retval_grad == nothing) ? nothing : retval_grad[key]
        kernel_arg_grad::Tuple = accumulate_param_gradients!(subtrace, kernel_retval_grad)
        @assert kernel_arg_grad[1] == nothing
        @assert kernel_arg_grad[2] == nothing
        for (i, (grad, has_grad)) in enumerate(zip(kernel_arg_grad[3:end], params_has_grad))
            if has_grad
                params_grad[i][key] = grad
            end
        end
    end
    (nothing, nothing, params_grad...)
end