text:
push!(LOAD_PATH,"../src/") # load module from local directory
using PolaronMobility

# Ivan Biaggio emailed 7th Nov 2018, suggesting there may be an error in these
# codes, due to a typographic error in the equations (particularly 62c) in
# Hellwarth1999
#
# As part of investigating this, these tests reproduce Table III in
# Hellwarth1999, to play with different ways of specifying the electron-phonon
# action.

"Print out F(alpha,beta) for a specific v,w; as a test"
function test_fns()
    @printf("\t\t")
    for α in 1:5
        @printf("α=%d\t\t",α)
    end
    @printf("\n")

    for β in 1:0.25:3.0
        v=w=4
        print("β: $β\t|\t")
        for α in 1:5
            @printf("%f\t",F(v,w,β,α))
        end
        println()
    end
end

# Our original definition of B, directly following the equation as written in Hellwarth1999
# 62c
println("Original B (62c, electron-phonon free energy):")
PolaronMobility.B(v,w,β,α) = α*v/(sqrt(π)*(exp(β)-1)) * quadgk(x->PolaronMobility.f(x,v,w,β),0,β/2)[1]
test_fns()  # OK - very primitive!
# Following private communication with Ivan Biaggio, and crosschecking against Osaka1959

println("Corrected B (62c, electron-phonon free energy, but corrected as Osaka.")
PolaronMobility.B(v,w,β,α) = exp(β)* α*v/(sqrt(π)*(exp(β)-1)) * quadgk(x->PolaronMobility.f(x,v,w,β),0,β/2)[1]
test_fns()

# Reproduce HellwarthTableIII, v and w solutions with finite Beta
function HellwarthTableIII()
    @printf("\t\t")
    for α in 1:5
        @printf("α=%d\t\t",α)
    end
    @printf("\n")

    for β in 1:0.25:3.0
        print("β: $β\t|\t")
        for α in 1:5
            v,w = feynmanvw(α, β)
            @printf("%.2f %.2f\t",v,w)
        end
        println()
    end

    print("Athermal\n β=Inf\t|\t")
    for α in 1:5
        v,w = feynmanvw(α)
        @printf("%.2f %.2f\t",v,w)
    end
    println()
end

# Our original definition of B, directly following the equation as written in Hellwarth1999
# 62c
println("Original B (62c, electron-phonon free energy):")
PolaronMobility.B(v,w,β,α) = α*v/(sqrt(π)*(exp(β)-1)) * quadgk(x->PolaronMobility.f(x,v,w,β),0,β/2)[1]
HellwarthTableIII()  # OK - very primitive!
# Following private communication with Ivan Biaggio, and crosschecking against Osaka1959

println("Corrected B (62c, electron-phonon free energy, but corrected as Osaka.")
PolaronMobility.B(v,w,β,α) = exp(β)* α*v/(sqrt(π)*(exp(β)-1)) * quadgk(x->PolaronMobility.f(x,v,w,β),0,β/2)[1]
HellwarthTableIII()

end

# Taking it back to where it all started; directly doing the full integral in Osaka1959
OsakaIntegrand(t,v,w,β) = exp(-t) / sqrt( w^2*t*(1-t/β) + (v^2-w^2)/v  *(1-exp(-v*t)+(1-coth(v/2*β)*(cosh(v*t)-1))) )
PolaronMobility.B(v,w,β,α) = 1/sqrt(pi) * α*v* exp(β)/(exp(β)-1) * quadgk(x->OsakaIntegrand(x,v,w,β),0,β)
# CURRENTLY BROKEN (!)
# sqrt(x) in the integrand is taken negative; so Julia throws a hissy fit - not expecting to find Imaginary values.
HellwarthTableIII()


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

module GeoEstimation

using Meshes
using GeoStatsBase
using Variography
using KrigingEstimators

using NearestNeighbors
using Distances
using LinearAlgebra

import GeoStatsBase: solve

include("idw.jl")
include("lwr.jl")
include("kriging.jl")

export IDW, LWR, Kriging

end

ext: jl


text:
module control

using Ipopt
using InfiniteOpt

import Parameters: @with_kw
using IOCapture: IOCapture
import Term: Panel, RenderableText
import Term.style: apply_style
import MyterialColors: blue_light

import jcontrol: Track
import ..bicycle: Bicycle, State

export ModelVariables, ControlOptions, create_and_solve_control, Bounds, State
export DynamicsProblem, KinematicsProblem, realistict_control_options

abstract type MTMproblem end

⋅ = *

# ---------------------------------------------------------------------------- #
#                                    OPTIONS                                   #
# ---------------------------------------------------------------------------- #
struct Bounds{T<:Number}
    lower::T
    upper::T
end

"""
Constructor with an additional parametr to convert angles in degrees to radians
"""
Bounds(lower, upper, angle) = Bounds(deg2rad(lower), deg2rad(upper))

function Base.show(io::IO, b::Bounds)
    return print(
        io,
        apply_style(
            "[green]$(round(b.lower; digits=2))[/green] [red]≤[/red] var [red]≤[/red][green] $(round(b.upper; digits=2))[/green]",
        ),
    )
end

"""
Options for solving the control problem.
Include options for the optimizer itself as well
other parameters such as bounds on allowed errors.
"""
@with_kw struct ControlOptions
    # errors bounds
    track_safety::Float64 = 1
    ψ_bounds::Bounds = Bounds(-30, 30, :angle)

    # control bounds
    u̇_bounds::Bounds = Bounds(-200, 200)
    δ̇_bounds::Bounds = Bounds(-1.5, 1.5)

    # varibles bounds
    u_bounds::Bounds = Bounds(5, 100)
    δ_bounds::Bounds = Bounds(-45, 45, :angle)
    ω_bounds::Bounds = Bounds(-500, 500, :angle)

    # dynamic system variables
    Fy_bounds::Bounds = Bounds(-100, 100)  # lateral forces
    Fu_bounds::Bounds = Bounds(-200, 200)  # driving force
    v_bounds::Bounds = Bounds(-125, 125)
end

realistict_control_options = ControlOptions(;
    u_bounds=Bounds(5, 80),
    u̇_bounds=Bounds(-180, 200),
    δ_bounds=Bounds(-50, 50, :angle),
    δ̇_bounds=Bounds(-4, 4),
    ω_bounds=Bounds(-400, 400, :angle),
    Fy_bounds=Bounds(-100, 100),
    v_bounds=Bounds(-125, 125),
)

"""
These are th best controls for the Dynamic model
as of 04/04/2022, they're also the very close
to the realistic values ranges.
"""
default_control_options = ControlOptions(;
    u_bounds=Bounds(10, 80),
    δ_bounds=Bounds(-50, 50, :angle),
    δ̇_bounds=Bounds(-6, 6),
    ω_bounds=Bounds(-400, 400, :angle),
    v_bounds=Bounds(-15, 15),
    Fu_bounds=Bounds(-2000, 3000),
)

# ---------------------------------------------------------------------------- #
#                                KINEMATIC MODEL                               #
# ---------------------------------------------------------------------------- #

"""
Model description.

Bicyle model with kinematics in the CoM reference frame.
    The CoM has position (x,y) and is at a distance l from
    the rear wheel, the total length is L.
    The front wheel has a steering angle δ and there's a slip
    angle β between the velocity vector at the CoM and the
    bike's orientation angle θ.

    Variables:
        L: total length
        l: rear wheel to CoM length

        x, y: position of CoM
        θ: orientation
        δ: steering angle
        β: slip angle of velocity vector
        u: velocity (speed)

    The model then has equations:
        ẋ = u * cos(θ + β)
        ẏ = u * sin(θ + β)
        θ̇ = ω = u * (tan(δ) * cos(β)) / L

    with:
        β = tan^{-1}(l * tan(δ)/L)

    There are two (bounded) controls:
        v̇
        δ̇

    The model is solved in the track's curvilinear coordinates system
    in which `n` represents the lateral distance from the center line
    and `ψ` the angular error (θ - track angle). The kinamtics need
    to be expressed as a function of distance along the track `s` instead
    of time `t`, that's done with the scalig factor

        SF = (1 - n * κ(s))/(u * cos(ψ + β) + eps())

    where κ(s) represents the track's curvature.

    Excluding x, y and θ's equation (because we are not doing things in the
    allocentric reference frame), we can write the model as:

        ∂(n, s) == SF * u * sin(ψ + β)
        ∂(ψ, s) == SF * ω - κ(s)

        ∂(u, s) == SF * v̇
        ∂(δ, s) == SF * δ̇
"""
struct KinematicsProblem <: MTMproblem end

"""
Create a `InfiniteOpt` model given a set of parameters and solve it to
get the controls for the MTM problem.
"""
function create_and_solve_control(
    problem_type::KinematicsProblem,
    num_supports::Int,
    track::Track,
    bike::Bicycle,
    options::ControlOptions,
    initial_conditions::State,
    final_conditions::State;
    quiet::Bool=false,
    n_iter::Int=1000,
    tollerance::Float64=1e-10,
    verbose::Int=0,
)
    # initialize optimizer
    model = InfiniteModel(Ipopt.Optimizer)
    set_optimizer_attribute(model, "max_iter", n_iter)
    set_optimizer_attribute(model, "acceptable_tol", tollerance)
    set_optimizer_attribute(model, "print_level", verbose)
    set_optimizer_attribute(model, "max_wall_time", 40.0)

    # register curvature function
    κ(s) = track.κ(s)
    @register(model, κ(s))

    # ----------------------------- define variables ----------------------------- #
    @infinite_parameter(model, s ∈ [0, track.S_f], num_supports = num_supports)

    @variables(
        model,
        begin
            # CONTROLS
            options.u̇_bounds.lower ≤ u̇ ≤ options.u̇_bounds.upper, Infinite(s)    # wheel acceleration
            options.δ̇_bounds.lower ≤ δ̇ ≤ options.δ̇_bounds.upper, Infinite(s)    # steering acceleration

            # track errors
            n, Infinite(s)  # constraints defined separtely
            options.ψ_bounds.lower ≤ ψ ≤ options.ψ_bounds.upper, Infinite(s)

            # other variables
            options.u_bounds.lower ≤ u ≤ options.u_bounds.upper, Infinite(s)
            options.δ_bounds.lower ≤ δ ≤ options.δ_bounds.upper, Infinite(s)

            β, Infinite(s)
            options.ω_bounds.lower ≤ ω ≤ options.ω_bounds.upper, Infinite(s)
            SF, Infinite(s)

            # time
            0 ≤ t ≤ 60, Infinite(s), (start = 10)
        end
    )

    # -------------------------- track width constraints ------------------------- #
    @parameter_function(model, allowed_track_width == track.width(s))
    @constraint(model, -allowed_track_width + bike.width ≤ n)
    @constraint(model, allowed_track_width - bike.width ≥ n)

    # ----------------------------- define kinematics ---------------------------- #
    l = bike.l_r
    L = bike.L

    @constraints(
        model,
        begin
            β == atan(l * tan(δ) / L)
            ω == u * (tan(δ) * cos(β)) / L
            SF == (1 - n * κ(s)) / (u * cos(ψ + β) + eps())  # time -> space domain conversion factor

            ∂(n, s) == SF * u * sin(ψ + β)
            ∂(ψ, s) == SF * ω - κ(s)

            ∂(u, s) == SF * u̇
            ∂(δ, s) == SF * δ̇

            ∂(t, s) == SF
        end
    )

    # ----------------------- set initial/final conditions ----------------------- #
    @constraints(
        model,
        begin
            # initial conditions
            n(0) == initial_conditions.n
            ψ(0) == initial_conditions.ψ
            u(0) == initial_conditions.u
            δ(0) == initial_conditions.δ
            t(0) == 0
            β(0) == initial_conditions.β
            ω(0) == initial_conditions.ω

            # final conditions
            u(track.S_f) == final_conditions.u
            ω(track.S_f) == final_conditions.ω
        end
    )

    # --------------------------------- optimize --------------------------------- #
    # solve
    @objective(model, Min, ∫(SF, s))
    optimize!(model)

    # print info
    if !quiet
        c = IOCapture.capture() do
            println(solution_summary(optimizer_model(model)))
        end
        print(
            "\n" *
            Panel(
                RenderableText(c.output, "$blue_light italic");
                style="yellow1",
                title="IPoPT output",
                title_style="red bold",
                justify=:center,
            ) *
            "\n\n",
        )
    end

    return model
end

# ---------------------------------------------------------------------------- #
#                                DYNAMICS MODEL                                #
# ---------------------------------------------------------------------------- #
"""
The dynamic bicycle model is taken from:
    Dynamics and optimal control of road Vehicles (Massaro 2018)

The model has the following variables
    x, y:   position
    θ:      orientation
    δ:      steering angle
    u:      longitudinal velocity component
    v:      lateral velocity component
    ω:      angular velocity (around the CoM)
    β:      the slip angle of the velocity vector V

and these forces
    Fu:     logitudinal force applied to back tire
    Fr:     lateral force at the rear wheel
    Ff:     lateral force at the front wheel

The bike has the following constants
    l_f:     distance from CoM to front  wheel
    l_r:     distance from rear wheel to CoM
    m:       mass
    Iz:      angular momentum
    c:  tyre cornering stiffness

The two controls are
    δ̇ and Fu

# --------------------------------- dynamics --------------------------------- #
The bike's dynamics are described by four equations of motion:
    δ̇    = δ̇   # the control
    mu̇   = mωv - Ff⋅sinδ + Fu
    mv̇   = -mωu + Ff⋅cosδ + Fr
    Iz⋅ω = l_f⋅Ff⋅cosδ - b⋅Fr

we can then compute the magnitude and slip angle of the velocity vector
    V̂ = √(u² + v²)
    β = arctan(v/u)

Finally the kinematics are given by:
    ẋ = V̂⋅cos(θ+β)
    ẏ = V̂⋅sin(θ+β)
    θ̇ = ω

# --------------------------------- solution --------------------------------- #
The strategy for the solution is the same as for the kinematics probelm.

Recast everything to the track's space using SF and using the track errors n and ψ,
like for the kinematics problem:
    ∂(n, s) == SF * u * sin(ψ + β)
    ∂(ψ, s) == SF * ω - κ(s)

With `SF` being the scaling factor:
    SF = (1 - n * κ(s))/(u * cos(ψ + β) + eps())
"""
struct DynamicsProblem <: MTMproblem end

"""
Create a `InfiniteOpt` model given a set of parameters and solve it to
get the controls for the MTM problem.
"""
function create_and_solve_control(
    problem_type::DynamicsProblem,
    num_supports::Int,
    track::Track,
    bike::Bicycle,
    options::ControlOptions,
    initial_conditions::State,
    final_conditions::Union{Nothing, State, Symbol};
    quiet::Bool=false,
    n_iter::Int=1000,
    tollerance::Float64=1e-4,
    verbose::Int=0,
    α::Float64=0.0, # cost of Fu
    γ::Float64=0.0,  # cost of δ̇
    waypoint=nothing,  # pass (s, State) to set the state at a track.s value
)

    # initialize optimizer
    model = InfiniteModel(Ipopt.Optimizer)
    set_optimizer_attribute(model, "max_iter", n_iter)
    set_optimizer_attribute(model, "acceptable_tol", tollerance)
    set_optimizer_attribute(model, "print_level", verbose)
    # set_optimizer_attribute(model, "max_wall_time", 180.0)

    # register curvature function
    κ(s) = track.κ(s)
    @register(model, κ(s))

    # ----------------------------- define variables ----------------------------- #
    @infinite_parameter(model, s ∈ [track.S[1], track.S[end]], num_supports = max(5, num_supports))

    @variables(
        model,
        begin
            # track errors
            n, Infinite(s)  # constraints defined separtely
            options.ψ_bounds.lower ≤ ψ ≤ options.ψ_bounds.upper, Infinite(s)

            # steering
            options.δ_bounds.lower ≤ δ ≤ options.δ_bounds.upper, Infinite(s)
            options.δ̇_bounds.lower ≤ δ̇ ≤ options.δ̇_bounds.upper, Infinite(s)    # control

            # long/lat/angular velocities
            options.u_bounds.lower ≤ u ≤ options.u_bounds.upper, Infinite(s)
            options.v_bounds.lower ≤ v ≤ options.v_bounds.upper, Infinite(s)
            options.ω_bounds.lower ≤ ω ≤ options.ω_bounds.upper, Infinite(s)

            # driving force
            options.Fu_bounds.lower ≤ Fu ≤ options.Fu_bounds.upper, Infinite(s)  # control

            # time
            0 ≤ t ≤ 60, Infinite(s), (start = 3)
       end
   )

    # -------------------------- track width constraints ------------------------- #
    @parameter_function(model, allowed_track_width == track.width(s))
    @constraint(model, -allowed_track_width + bike.width ≤ n)
    @constraint(model, allowed_track_width - bike.width ≥ n)

    # ----------------------------- define EOM       ---------------------------- #
    l_r, l_f = bike.l_r, bike.l_f
    m, Iz, c = bike.m, bike.Iz, bike.c

    β = atan(v / (u))  # slip angle
    V = √(u^2 + v^2)
    SF = (1 - n * κ(s)) / (V ⋅ cos(ψ + β))  # time -> space domain conversion factor

    # lateral forces
    Ff = c⋅(δ - (l_f⋅ω + v)/u)
    Fr = c⋅(l_r⋅ω - v)/u

    @constraints(
        model,
        begin
            # errors
            ∂(n, s) == SF * u ⋅ sin(ψ + β)
            ∂(ψ, s) == SF * ω - κ(s)

            # EOM
            ∂(δ, s) == SF * δ̇
            ∂(u, s) == SF / m * (m⋅ω⋅v - Ff⋅sin(δ) + Fu)
            ∂(v, s) == SF / m * (-m⋅ω⋅u + Ff⋅cos(δ) + Fr)
            ∂(ω, s) == SF / Iz * (l_f⋅Ff⋅cos(δ) - l_r⋅Fr)

            # time
            ∂(t, s) == SF
        end
    )

    # ----------------------- set initial/final conditions ----------------------- #
    # initial conditions
    @constraints(
        model,
        begin
            n(track.S[1]) == initial_conditions.n
            ψ(track.S[1]) == initial_conditions.ψ

            δ(track.S[1]) == initial_conditions.δ
            δ̇(track.S[1]) == initial_conditions.δ̇

            u(track.S[1]) == initial_conditions.u
            v(track.S[1]) == initial_conditions.v
            ω(track.S[1]) == initial_conditions.ω
            Fu(track.S[1]) == initial_conditions.Fu

            t(track.S[1]) == 0
        end
    )

    # final conditions
    if !isnothing(final_conditions)
        if final_conditions == :minimal
            @constraints(
                model,
                begin
                n(track.S[end]) == 0
                ψ(track.S[end]) == 0

                end
            )
        else
            @constraints(
                model,
                begin
                n(track.S[end]) == final_conditions.n
                ψ(track.S[end]) == final_conditions.ψ

                δ(track.S[end]) == final_conditions.δ
                δ̇(track.S[end]) == final_conditions.δ̇

                u(track.S[end]) == final_conditions.u
                v(track.S[end]) == final_conditions.v
                ω(track.S[end]) == final_conditions.ω
                Fu(track.S[end]) == final_conditions.Fu
                end
            )
            end
        end

    # additional waypoint
    if !isnothing(waypoint)
        sval, constraint = waypoint
        @constraints(
            model,
            begin
            n(sval) == constraint.n
            ψ(sval) == constraint.ψ

            δ(sval) == constraint.δ
            δ̇(sval) == constraint.δ̇

            u(sval) == constraint.u
            v(sval) == constraint.v
            ω(sval) == constraint.ω
            Fu(sval) == constraint.Fu
            end
        )
    end

    # --------------------------------- optimize --------------------------------- #
    set_all_derivative_methods(model, OrthogonalCollocation(2))
    @objective(model, Min, ∫(SF + α*Fu + γ*δ̇, s))
    optimize!(model)

    # print info
    if !quiet
        c = IOCapture.capture() do
            println(solution_summary(optimizer_model(model)))
        end
        print(
            "\n" *
            Panel(
                c.output;
                style="yellow1",
                title="IPoPT output",
                title_style="red bold",
                justify=:center,
            ) *
            "\n\n",
        )
    end
    return model
end

end

max_stars_repo_name: FedeClaudi/LocomotionControl

max_issues_repo_name: FedeClaudi/LocomotionControl

max_forks_repo_name: FedeClaudi/LocomotionControl


text:
#TODO: clean up this gawdawful struct
mutable struct ROSPF{C, P, D, POPUB, PAPUB, ODPUB, M, LP, LP2}
    conf::C
    pf::P
    dtype::D

    pose_pub::POPUB
    particle_pub::PAPUB
    odom_pub::ODPUB

    occmap::M

    # scratch
    laserpose_pub::LP

    mr_W::Float64
    maxviz::Int
    nrays::Int

    last_laser::Float64
    last_acker::Float64
    last_vesc::Float64
    last_servo::Float64

    speed2erpm_offset::Float64
    speed2erpm_gain::Float64
    steering2servo_offset::Float64
    steering2servo_gain::Float64

    initialized::Bool

    rand_part_idxs::Vector{Int}
    angles::Vector{Float64}
    laseridxs::Vector{Int}
    laserposes::LP2
end

function ROSPF(conf, pf, pose_pub, particle_pub, odom_pub, occmap, laserposepub)
    mr_W = conf[:pfconf][:reweightconf][:max_range_meters]
    maxviz = conf[:rosconf][:max_viz_particles]
    nrays = conf[:pfconf][:reweightconf][:nrays]
    ridxs = Vector{Int}(undef, maxviz)
    angles = Vector{Float64}(undef, 1)
    laseridxs = Vector{Int}(undef, 1)
    spo = conf

    spo = conf[:rosconf][:vescconf][:speed2erpm_offset]
    spg = conf[:rosconf][:vescconf][:speed2erpm_gain]
    sto = conf[:rosconf][:vescconf][:steering2servo_offset]
    stg = conf[:rosconf][:vescconf][:steering2servo_gain]

    laserposes = Vector{Pose2D{conf[:pfconf][:dtype]}}(undef, nrays)
    ROSPF(conf, pf, conf[:pfconf][:dtype], pose_pub, particle_pub, odom_pub, occmap, laserposepub, convert(Float64, mr_W), maxviz, nrays, -1.0, -1.0, -1.0,-1.0, spo, spg, sto, stg, false, ridxs,angles,laseridxs, laserposes)
end

macro debugtask(ex)
  quote
    try
      $(esc(ex))
    catch e
      bt = stacktrace(catch_backtrace())
      io = IOBuffer()
      showerror(io, e, bt)
      errstr = String(take!(io))
      RobotOS.logfatal("Error: $errstr")
      exit()
    end
  end
end

function laserposes!(rospf, angles, ranges, meanpose)
    @inbounds for (i, (a, r)) in enumerate(zip(angles, ranges))
        x, y, theta = meanpose.statev
        theta += a
        xhit, yhit = rangebearing2point(x, y, theta, r)
        rospf.laserposes[i] = Pose2D{rospf.dtype}(xhit, yhit, theta)
    end
    rospf
end

function rosheader(frame_id, stamp=nothing)
    header = Header()
    header.frame_id = frame_id
    header.stamp = stamp===nothing ? get_rostime() : stamp
    header
end

@inline rosheader!(msg, ::Nothing, ::Nothing) = msg
@inline function rosheader!(msg, frame_id=nothing, stamp=nothing)
    msg.header = rosheader(frame_id, stamp)
    msg
end

@inline function rosodometrystamped(pose::Pose2D, frame_id, stamp=nothing)
    odom = Odometry()
    odom.pose = roscovarpose(pose)
    rosheader!(odom, frame_id, stamp)
end

@inline function roscovarposestamped(pose::Pose2D, frame_id, stamp=nothing)
    rcps = PoseWithCovarianceStamped()
    rcps.pose = roscovarpose(pose)
    rosheader!(rcps, frame_id, stamp)
end

@inline function roscovarpose(pose::Pose2D)
    covarpose = PoseWithCovariance()
    covarpose.pose = rospose(pose)
    covarpose
end

@inline function rosposestamped(pose::Pose2D, frame_id, stamp=nothing)
    rps = PoseStamped()
    rps.pose = rospose(pose)
    rosheader!(rps, frame_id, stamp)
end

@inline function rospose(pose::Pose2D)
    po = Pose()
    po.position.x = pose.x
    po.position.y = pose.y
    q=Quat(RotXYZ(0, 0, pose.theta))
    po.orientation = Quaternion(q.x, q.y, q.z, q.w)
    po
end

function validindices(ranges::Vector, rospf::ROSPF)
    j = 0
    @inbounds for i in eachindex(ranges)
        r = ranges[i]
        if !isnan(r) && 0<r<rospf.mr_W
            j += 1
            rospf.laseridxs[j] = i
        end
    end
    @assert 0< j <= length(ranges)
    view(rospf.laseridxs, 1:j)
end

function rosposearraystamped(ps::AbstractVector, frame_id, stamp=nothing)
    arr = PoseArray()
    arr.poses = Vector{Pose}(undef, length(ps))
    @inbounds @simd for i in eachindex(arr.poses)
        arr.poses[i] = rospose(ps[i])
    end
    rosheader!(arr, frame_id, stamp)
end


text:
# Wrapper for level-1 BLIS vector routines.
#

# Level1v API of common forms are put together.
#
macro blis_group_level1v_form1(funcname)
    return quote
        @blis_ccall_group($funcname,
                          Cvoid,
                          conjx, BliConj,
                          n,     BliDim,
                          x,     Ptr{xType}, incx, BliInc,
                          y,     Ptr{xType}, incy, BliInc)
    end
end

@blis_group_level1v_form1 addv
@blis_group_level1v_form1 copyv
@blis_group_level1v_form1 subv

macro blis_group_level1v_form2(funcname)
    return quote
        @blis_ccall_group($funcname,
                          Cvoid,
                          conjx, BliConj,
                          n,     BliDim,
                          α,     Ptr{xType},
                          x,     Ptr{xType}, incx, BliInc,
                          y,     Ptr{xType}, incy, BliInc)
    end
end

@blis_group_level1v_form2 axpyv
@blis_group_level1v_form2 scal2v

macro blis_group_level1v_form3(funcname)
    return quote
        @blis_ccall_group($funcname,
                          Cvoid,
                          conjα, BliConj,
                          n,     BliDim,
                          α,     Ptr{xType},
                          x,     Ptr{xType}, incx, BliInc)
    end
end

@blis_group_level1v_form3 scalv
@blis_group_level1v_form3 setv

# Forms with only 1 member are defined directly.
#
@blis_ccall_group(amaxv,
                  Cvoid,
                  n,     BliDim,
                  x,     Ptr{xType}, incx, BliInc,
                  index, Ptr{Clong})
@blis_ccall_group(axpbyv,
                  Cvoid,
                  conjx, BliConj,
                  n,     BliDim,
                  α,     Ptr{xType},
                  x,     Ptr{xType}, incx, BliInc,
                  beta,  Ptr{xType},
                  y,     Ptr{xType}, incy, BliInc)
@blis_ccall_group(dotv,
                  Cvoid,
                  conjx, BliConj,
                  conjy, BliConj,
                  n,     BliDim,
                  x,     Ptr{xType}, incx, BliInc,
                  y,     Ptr{xType}, incy, BliInc,
                  ρ,     Ptr{xType})
@blis_ccall_group(dotxv,
                  Cvoid,
                  conjx, BliConj,
                  conjy, BliConj,
                  n,     BliDim,
                  α,     Ptr{xType},
                  x,     Ptr{xType}, incx, BliInc,
                  y,     Ptr{xType}, incy, BliInc,
                  β,     Ptr{xType},
                  ρ,     Ptr{xType})
@blis_ccall_group(invertv,
                  Cvoid,
                  n,     BliDim,
                  x,     Ptr{xType}, incx, BliInc)
@blis_ccall_group(swapv,
                  Cvoid,
                  n,     BliDim,
                  x,     Ptr{xType}, incx, BliInc,
                  y,     Ptr{xType}, incy, BliInc)
@blis_ccall_group(xpbyv,
                  Cvoid,
                  conjx, BliConj,
                  n,     BliDim,
                  x,     Ptr{xType}, incx, BliInc,
                  β,     Ptr{xType},
                  y,     Ptr{xType}, incy, BliInc)


text:
module problem129

description = """
Repunit divisibility

A number consisting entirely of ones is called a repunit. We shall define R(k) to be a repunit of length k; for example, R(6) = 111111.
Given that n is a positive integer and GCD(n, 10) = 1, it can be shown that there always exists a value, k, for which R(k) is divisible by n,
and let A(n) be the least such value of k; for example, A(7) = 6 and A(41) = 5.
The least value of n for which A(n) first exceeds ten is 17.

Find the least value of n for which A(n) first exceeds one-million.
"""

# Calculate repunit length k = (10^k - 1) / 9
function R(k :: Int)
  return div(BigInt(10)^k - 1, 9)
end

# find the least value of k such that R(k) is a multiple of n
# Insight from the Project Euler forums:
# A(n) is the minimum k, such that n divides R(k)
# => (10^k - 1)/9 == 0 mod n
# => 10^k - 1     == 0 mod 9n
# => 10^k         == 1 mod 9n
# Iteratively generate powers of 10 mod 9n, and terminate when the result is 1.
function A(n)
  tenK = 10
  for k = 1:Inf
    if tenK == 1
      return k
    end
    tenK = (tenK * 10) % (9*n)
  end
end

# Search for the least value of n, such that A(n) > x
# The main insight here is that search(x) >= x, so search 1,000,000 doesn't require
#  searching 1..1000000.
function search(x = 1000000)
  for n = x:Inf
    if gcd(n, 10) == 1
      k = A(n)
      if k > x
        return n
      end
    end
  end
end

using Base.Test

@test R(6) == 111111
@test A(7) == 6
@test A(41) == 5
@test search(10) == 17

end

max_stars_repo_name: mbuhot/mbuhot-euler-solutions

max_issues_repo_name: mbuhot/mbuhot-euler-solutions

max_forks_repo_name: mbuhot/mbuhot-euler-solutions


text:
module NIFTy

import Base: size, length, eltype, getindex, setindex!, IndexStyle, randn, similar
import Base: adjoint
import Base: exp, log, sin, cos, +, -, *, /

import FFTW: fft, ifft
using FFTW

using LoopVectorization

include("Domain.jl")
include("Field.jl")
include("Operator.jl")
include("Functions.jl")
include("Distributor.jl")
include("OperatorLibrary.jl")

export AbstractDomain, Domain, UnstructuredDomain, StructuredDomain,
       RGDomain, PowerDomain, DomainTuple,
       rgdomain,
       lengths, isharmonic, distances, dims, shape,
       getcodomain,
       domains,
       #Field.jl
       Field, MultiField,
       similar, getdomain,
       #Functions.jl
       hartley, ihartley,
       #Operator.jl
       AbstractOperator, LinearOperator, PointwiseOperator, OperatorChain, Operator,
       domain, target, adjoint,
       combine_operators, apply,
       #Distributor.jl
       DistributionOperator, PowerDistributor,
       #OperatorLibrary.jl
       add, scale,
       AdditionOperator, SubtractionOperator, MultiplicationOperator, DivisionOperator,
       HartleyOperator, FourierOperator
end


text:
using PICDataStructures, Test
using PICDataStructures: dimensionality
using Unitful
using CairoMakie

@testset "Scalar fields" begin
    grids =
        SparseAxisGrid(0.01:0.1:1),
        SparseAxisGrid(-1:0.012:1, -1:0.012:1),
        SparseAxisGrid(-1:0.0052:1, -1:0.012:1, -1:0.012:1),
        SparseAxisGrid(0.01u"m":0.1u"m":1.0u"m"),
        SparseAxisGrid(-1u"m":0.012u"m":1u"m", -1u"m":0.012u"m":1u"m"),
        SparseAxisGrid(-1u"m":0.0052u"m":1u"m", -1u"m":0.012u"m":1u"m", -1u"m":0.012u"m":1u"m")

    fields = [
        scalarfield(x->inv(x...), grids[1]),
        scalarfield(grids[2]) do (x,y)
            1 / √(x^2 + y^2)
        end,
        scalarfield(grids[3]) do (x,y,z)
            1 / √(x^2 + y^2 + z^2)
        end,
        scalarfield(grids[4]) do (x,)
            1u"V" / x
        end,
        scalarfield(grids[5]) do (x,y)
            1u"V" / √(x^2 + y^2)
        end,
        scalarfield(grids[6]) do (x,y,z)
            1u"V" / √(x^2 + y^2 + z^2)
        end
    ]

    @testset "Unitless" begin
        f = fields[1]
        fig = plotdata(f)
        @test fig.current_axis[].xlabel[] == "x"

        f = fields[2]
        fig = plotdata(f)
        @test fig.current_axis[].xlabel[] == "x"
        @test fig.current_axis[].ylabel[] == "y"

        f = fields[3]
        fig = plotdata(f)
        @test fig.current_axis[].scene[OldAxis][:names, :axisnames][] == ("x", "y", "z")
    end

    @testset "Unitful" begin
        f = fields[4]
        fig = plotdata(f)
        u = unitname(getdomain(f))
        @test u == "m"
        @test fig.current_axis[].xlabel[] == "x ($u)"

        f = fields[5]
        fig = plotdata(f)
        u = unitname(getdomain(f))
        @test fig.current_axis[].xlabel[] == "x ($u)"
        @test fig.current_axis[].ylabel[] == "y ($u)"

        f = fields[6]
        fig = plotdata(f)
        u = unitname(getdomain(f))
        # VolumeLike doesn't work on CairoMakie
        @test_broken fig.current_axis[].scene[OldAxis][:names, :axisnames][] == ("x ($u)", "y ($u)", "z ($u)")
    end
end

@testset "Scalar variables" begin
    grids = [
        ParticlePositions(sin.(0:0.1:2π)),
        ParticlePositions(sin.(0:0.1:2π),cos.(0:0.1:2π)),
        ParticlePositions(sin.(0:0.1:2π),cos.(0:0.1:2π),sin.(0:0.1:2π)),
        ParticlePositions(sin.(0:0.1:2π).*u"m"),
        ParticlePositions(sin.(0:0.1:2π).*u"m",cos.(0:0.1:2π).*u"m"),
        ParticlePositions(sin.(0:0.1:2π).*u"m",cos.(0:0.1:2π).*u"m",sin.(0:0.1:2π).*u"m"),
    ]
    vars = [
        scalarvariable(grids[1]) do (x,)
            x^2
        end,
        scalarvariable(grids[2]) do (x,y)
            x^2 + y^2
        end,
        scalarvariable(grids[3]) do (x,y,z)
            (x^2 + y^2 + z^2)
        end,
        scalarvariable(grids[4]) do (x,)
            x^2
        end,
        scalarvariable(grids[5]) do (x,y)
            (x^2 + y^2)
        end,
        scalarvariable(grids[6]) do (x,y,z)
            (x^2 + y^2 + z^2)
        end,
    ]

    @testset "Unitless" begin
        v = vars[1]
        fig = plotdata(v)
        @test fig.current_axis[].xlabel[] == "x"

        v = vars[2]
        fig = plotdata(v)
        @test fig.current_axis[].xlabel[] == "x"
        @test fig.current_axis[].ylabel[] == "y"
        @test all(fig.content[2].limits[] .≈ (-1,1))

        v = vars[3]
        fig = plotdata(v)
        @test fig.current_axis[].scene[OldAxis][:names, :axisnames][] == ("x", "y", "z")
    end

    @testset "Unitful" begin
        v = vars[4]
        fig = plotdata(v)
        u = unitname(getdomain(v))
        @test u == "m"
        @test fig.current_axis[].xlabel[] == "x ($u)"

        v = vars[5]
        fig = plotdata(v)
        u = unitname(getdomain(v))
        @test u == "m"
        @test fig.current_axis[].xlabel[] == "x ($u)"
        @test fig.current_axis[].ylabel[] == "y ($u)"
        @test all(fig.content[2].limits[] .≈ (-1,1))

        v = vars[6]
        fig = plotdata(v)
        @test_broken fig.current_axis[].scene[OldAxis][:names, :axisnames][] == ("x ($u)", "y ($u)", "z ($u)")
    end
end

@testset "Vector fields" begin
    grids =
        SparseAxisGrid(0.01:0.1:1),
        SparseAxisGrid(-1:0.012:1, -1:0.012:1),
        SparseAxisGrid(-1:0.0052:1, -1:0.012:1, -1:0.012:1),
        SparseAxisGrid(0.01u"m":0.1u"m":1.0u"m"),
        SparseAxisGrid(-1u"m":0.012u"m":1u"m", -1u"m":0.012u"m":1u"m"),
        SparseAxisGrid(-1u"m":0.0052u"m":1u"m", -1u"m":0.012u"m":1u"m", -1u"m":0.012u"m":1u"m")

    scalar_fields = [
        scalarfield(x->inv(x...), grids[1]),
        scalarfield(grids[2]) do (x,y)
            1 / √(x^2 + y^2)
        end,
        scalarfield(grids[3]) do (x,y,z)
            1 / √(x^2 + y^2 + z^2)
        end,
        scalarfield(grids[4]) do (x,)
            1u"V" / x
        end,
        scalarfield(grids[5]) do (x,y)
            1u"V" / √(x^2 + y^2)
        end,
        scalarfield(grids[6]) do (x,y,z)
            1u"V" / √(x^2 + y^2 + z^2)
        end
    ]
    fields = [
        build_vector((scalar_fields[1], ), (:x,)),
        build_vector((scalar_fields[2], scalar_fields[2]), (:x, :y)),
        build_vector((scalar_fields[3], scalar_fields[3], scalar_fields[3]), (:x, :y, :z)),
        build_vector((scalar_fields[4], ), (:x,)),
        build_vector((scalar_fields[5], scalar_fields[5]), (:x, :y)),
        build_vector((scalar_fields[6], scalar_fields[6], scalar_fields[6]), (:x, :y, :z))
    ]

    @testset "Unitless" begin
        # Does 1D make sense?
        # f = fields[1]
        # fig = plotdata(f)
        # @test fig.current_axis[].xlabel[] == "x"

        f = fields[2]
        fig = plotdata(f)
        @test fig.current_axis[].xlabel[] == "x"
        @test fig.current_axis[].ylabel[] == "y"

        f = fields[3]
        fig = plotdata(f)
        @test fig.current_axis[].scene[OldAxis][:names, :axisnames][] == ("x", "y", "z")
    end

    @testset "Unitful" begin
        # f = fields[4]
        # fig = plotdata(f)
        # u = unitname(getdomain(f))
        # @test u == "m"
        # @test fig.current_axis[].xlabel[] == "x ($u)"

        f = fields[5]
        fig = plotdata(f)
        u = unitname(getdomain(f))
        @test fig.current_axis[].xlabel[] == "x ($u)"
        @test fig.current_axis[].ylabel[] == "y ($u)"

        f = fields[6]
        fig = plotdata(f)
        u = unitname(getdomain(f))
        # VolumeLike doesn't work on CairoMakie
        @test_broken fig.current_axis[].scene[OldAxis][:names, :axisnames][] == ("x ($u)", "y ($u)", "z ($u)")
    end
end

@testset "Vector variables" begin
    grids = [
        ParticlePositions(sin.(0:0.1:2π)),
        ParticlePositions(sin.(0:0.1:2π),cos.(0:0.1:2π)),
        ParticlePositions(sin.(0:0.1:2π),cos.(0:0.1:2π),sin.(0:0.1:2π)),
        ParticlePositions(sin.(0:0.1:2π).*u"m"),
        ParticlePositions(sin.(0:0.1:2π).*u"m",cos.(0:0.1:2π).*u"m"),
        ParticlePositions(sin.(0:0.1:2π).*u"m",cos.(0:0.1:2π).*u"m",sin.(0:0.1:2π).*u"m"),
    ]
    scalar_vars = [
        scalarvariable(grids[1]) do (x,)
            x^2
        end,
        scalarvariable(grids[2]) do (x,y)
            (x)/10
        end,
        scalarvariable(grids[3]) do (x,y,z)
            (x^2 + y^2 + z^2)/10
        end,
        scalarvariable(grids[4]) do (x,)
            x^2
        end,
        scalarvariable(grids[5]) do (x,y)
            (x)/10
        end,
        scalarvariable(grids[6]) do (x,y,z)
            (x^2 + y^2 + z^2)/10
        end,
    ]
    vars = [
        build_vector((scalar_vars[1], ), (:x,)),
        build_vector((scalar_vars[2], scalar_vars[2]), (:x, :y)),
        build_vector((scalar_vars[3], scalar_vars[3], scalar_vars[3]), (:x, :y, :z)),
        build_vector((scalar_vars[4], ), (:x,)),
        build_vector((scalar_vars[5], scalar_vars[5]), (:x, :y)),
        build_vector((scalar_vars[6], scalar_vars[6], scalar_vars[6]), (:x, :y, :z))
    ]

    @testset "Unitless" begin
        # v = vars[1]
        # fig = plotdata(v)
        # @test fig.current_axis[].xlabel[] == "x"

        v = vars[2]
        fig = plotdata(v)
        @test fig.current_axis[].xlabel[] == "x"
        @test fig.current_axis[].ylabel[] == "y"
        @test all(fig.content[2].limits[] .≈ (0,√2/10))

        v = vars[3]
        fig = plotdata(v)
        @test fig.current_axis[].scene[OldAxis][:names, :axisnames][] == ("x", "y", "z")
    end

    @testset "Unitful" begin
        # v = vars[4]
        # fig = plotdata(v)
        # u = unitname(getdomain(v))
        # @test u == "m"
        # @test fig.current_axis[].xlabel[] == "x ($u)"

        v = vars[5]
        fig = plotdata(v)
        u = unitname(getdomain(v))
        @test u == "m"
        @test fig.current_axis[].xlabel[] == "x ($u)"
        @test fig.current_axis[].ylabel[] == "y ($u)"
        @test all(fig.content[2].limits[] .≈ (0,√2/10))

        v = vars[6]
        fig = plotdata(v)
        @test_broken fig.current_axis[].scene[OldAxis][:names, :axisnames][] == ("x ($u)", "y ($u)", "z ($u)")
    end
end


text:
# Autogenerated wrapper script for libwebp_jll for i686-linux-gnu
export cwebp, dwebp, gif2webp, img2webp, libwebp, libwebpdecoder, libwebpdemux, libwebpmux, webpinfo, webpmux

using Giflib_jll
using JpegTurbo_jll
using libpng_jll
using Libtiff_jll
using Libglvnd_jll
JLLWrappers.@generate_wrapper_header("libwebp")
JLLWrappers.@declare_executable_product(cwebp)
JLLWrappers.@declare_executable_product(dwebp)
JLLWrappers.@declare_executable_product(gif2webp)
JLLWrappers.@declare_executable_product(img2webp)
JLLWrappers.@declare_library_product(libwebp, "libwebp.so.7")
JLLWrappers.@declare_library_product(libwebpdecoder, "libwebpdecoder.so.3")
JLLWrappers.@declare_library_product(libwebpdemux, "libwebpdemux.so.2")
JLLWrappers.@declare_library_product(libwebpmux, "libwebpmux.so.3")
JLLWrappers.@declare_executable_product(webpinfo)
JLLWrappers.@declare_executable_product(webpmux)
function __init__()
    JLLWrappers.@generate_init_header(Giflib_jll, JpegTurbo_jll, libpng_jll, Libtiff_jll, Libglvnd_jll)
    JLLWrappers.@init_executable_product(
        cwebp,
        "bin/cwebp",
    )

    JLLWrappers.@init_executable_product(
        dwebp,
        "bin/dwebp",
    )

    JLLWrappers.@init_executable_product(
        gif2webp,
        "bin/gif2webp",
    )

    JLLWrappers.@init_executable_product(
        img2webp,
        "bin/img2webp",
    )

    JLLWrappers.@init_library_product(
        libwebp,
        "lib/libwebp.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libwebpdecoder,
        "lib/libwebpdecoder.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libwebpdemux,
        "lib/libwebpdemux.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libwebpmux,
        "lib/libwebpmux.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        webpinfo,
        "bin/webpinfo",
    )

    JLLWrappers.@init_executable_product(
        webpmux,
        "bin/webpmux",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using EngEconomics, Roots, Plots, ColorSchemes
plotly()

uwEngColorScheme = ColorScheme([parse(Colorant, "#D0B4E7"), parse(Colorant, "#57058B"), parse(Colorant, "#8100B4"), parse(Colorant, "#000000")])

# Given
n = 5
timeDelay = 3
A = 10
G = 1
MARR = 0.15
sensitivity = [0.85, 0.925, 1.0, 1.075, 1.15]

# Find
# Construct a sensitivity graph for the interest rate:
NPW_interest = (A .+ G * sinkingFundPaymentArithmetic.(MARR * sensitivity, n + 1)) .* uniformSeriesCompoundAmoundFactor.(MARR * sensitivity, (n + 1)) .* presentWorthFactor.(MARR * sensitivity, n + timeDelay)
plot(sensitivity, NPW_interest, label="Interest", legend=:outerbottomright,
	xlabel="Sensitivity (+/- 7.5%, 15%)",
	ylabel="Net Present Worth (\$1000)",
	title="Sensitivity of Interest Rate, Base Savings,<br> and Saving Gradient",
	palette=cgrad(uwEngColorScheme))

# Construct a sensitivity graph for the base savings
NPW_baseSavings = (A * sensitivity .+ G * sinkingFundPaymentArithmetic.(MARR, n + 1)) .* uniformSeriesCompoundAmoundFactor.(MARR, n + 1) .* presentWorthFactor.(MARR, n + timeDelay)
plot!(sensitivity, NPW_baseSavings, label="Base Savings")

# Construct a sensitivity graph for the savings gradient
NPW_gradientSavings = (A .+ G * sensitivity * sinkingFundPaymentArithmetic.(MARR, n + 1)) .* uniformSeriesCompoundAmoundFactor.(MARR, n + 1) .* presentWorthFactor.(MARR, n + timeDelay)
plot!(sensitivity, NPW_gradientSavings, label="Gradient Savings")


text:
#load the ff data
function loadff(ffname::String;
  ffpath::String = PARAM[:ffpath],
  ffdateformat::DateFormat = PARAM[:wrdsdateformat],
  fieldprefix::String = "F_$(lowercase(ffname))",
  ffscale::Bool = PARAM[:ffscale][ffname])

  ff::DataFrame = CSV.File("$ffpath\\$ffname.csv") |> DataFrame

  #rename the columns
  ffnames::Vector{String} = names(ff)
  for (i,str) ∈ enumerate(ffnames)
    if i==1 #date column is a special case
      ffnames[1] = "date"
      continue
    end

    str = replace(str, "-"=>"")
    str = replace(str, " "=>"")
    str = replace(str, "_"=>"")
    str = lowercase(str)
    ffnames[i] = length(fieldprefix) ≥ 1 ? string("$(fieldprefix)_", str) : str
  end
  rename!(ff, ffnames)
  ff.date = (i::Int->Date("$i", ffdateformat)).(ff.date)

  #check for missing data and scale
  Frets::Vector{String} = setdiff(ffnames, ["date"])
  Threads.@threads for Fret ∈ Frets #iterate the focal columns
    for r ∈ eachrow(ff) #now iterate down the rows
      if (r[Fret] == -99.99) || (r[Fret] < -100.) #missing codes for ff
        r[Fret] = missing
      elseif ffscale
        r[Fret] /= 100.
      end
    end
  end

  return ff
end

#loads the ff returns and converts the returns to an index
function loadffindexes(ffname::String)
  rawff::DataFrame = loadff(ffname)
  sort!(rawff, :date)

  #want the return names which include the file prefix
  retnames::Vector{Symbol} = (s->Symbol("F_$(ffname)_", s)).(PARAM[:ffretfieldidx][ffname])
  Findexes, indexdf::DataFrame = returns2index(rawff,Frets=retnames)

  #writes out the entire df and index
  if PARAM[:testff]
    testdf = leftjoin(rawff, indexdf, on=:date)
    sort!(testdf, :date)
    testdf |> CSV.write("$(PARAM[:testpath])\\$(ffname)_idxs.csv")
  end

  return Findexes, indexdf
end

#loads ff data and merges it into the existing dataframe
function mergeff(df::DataFrame;
  datefrequency::Symbol = throw("date frequency is required (:day, :month, :week)"),
  testlabel::String = "FF",
  Fgroup::NSymbol = nothing, #only used for testing purposes if provided
  validateexactfrequency = datefrequency != :day)

  local ffnames::Vector{String}=PARAM[:ffnames]
  local ffretfieldidx::Dict = PARAM[:ffretfieldidx]

  #load the idnex values, perform a basic validation and merge them into the main df
  for (i,ffname) ∈ enumerate(ffnames)
    Findexes, indexdf::DataFrame = loadffindexes(ffname)

    #make sure we have the right fields
    Frets = (s->Symbol("F_$(ffname)_", s)).(ffretfieldidx[ffname])
    @assert issetequal(Findexes, (Fret->Symbol(Fret,:_idx)).(Frets))
    @assert setdiff(Findexes, propertynames(indexdf)) |> isempty

    df = index2returns!(df, indexdf;
      Frets, Findexes, datefrequency, testlabel, Fgroup, validateexactfrequency)

    if datefrequency === :month
      select!(df, Not([:indexdate]))
    end
  end

  return df
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
import Base: read

@enum MirrorModes begin
  Horizontal = 0
  Vertical = 1
  Single0 = 2
  Single1 = 3
  Four = 4
end

mutable struct Cartridge
  # PRG-ROM banks
  PRG::Vector{UInt8}
  # CHR-ROM banks
  CHR::Vector{UInt8}
  # Save RAM
  SRAM::Vector{UInt8}
  # mapper type
  mapper::UInt8
  # mirroring mode
  mirror::UInt8
  ## battery present
  battery::UInt8

  function Cartridge(PRG::Vector{UInt8}, CHR::Vector{UInt8}, mapper::UInt8, mirror::UInt8, battery::UInt8)
    new(PRG, CHR, zeros(UInt8, 0x2000), mapper, mirror, battery)
  end
end

const iNESFileMagic = 0x1a53454e

function fromfile(filename::String)::Cartridge
  open(filename) do f
    # read file header
    magic = read(f, UInt32)
    numPRG = read(f, UInt8)
    numCHR = read(f, UInt8)
    control1 = read(f, UInt8)
    control2 = read(f, UInt8)
    numRAM = read(f, UInt8)
    padding = read(f, 7)
    if magic != iNESFileMagic
      throw(ErrorException("File does not start with magic number."))
    end

    # mapper type
    mapper1 = control1 >> 0x04
    mapper2 = control2 >> 0x04
    mapper = mapper1 | (mapper2 << 0x04)

    # mirroring type
    mirror1 = control1 & 0x01
    mirror2 = (control1 >> 0x03) & 0x01
    mirror = mirror1 | (mirror2 << 0x01)

    # battery-backed RAM
    battery = (control1 >> 0x01) & 0x01

    # read trainer if present (unused)
    if (control1 & 0x04) == 0x04
      throw(ErrorException("Shouldn't be here."))
      # todo
    end

    # read prg-rom bank(s)
    prgLength = numPRG * 16384
    prg = read(f, prgLength)

    # read chr-rom bank(s)
    chrLength = 8192
    if numCHR != 0
      chrLength = numCHR * 8192
      chr = read(f, chrLength)
    else
      chr = zeros(UInt8, chrLength)
    end

    # success
    Cartridge(prg, chr, mapper, mirror, battery)
  end
end


text:
using Coverage, Pkg, Glob
rm.(glob("src/*.cov"))
Pkg.test(;coverage=true)
coverage = process_folder()
covered_lines, total_lines = get_summary(coverage)
println("Coverage: ", round(100.0*covered_lines/total_lines), " %")


text:
using Dictionaries
using Test

# The following tests are run without the --bounds-check option.

# Issue #63 - Broadcasting only omitted bounds check due to inbounds propagation
@static if VERSION < v"1.4" # Duplicate of Base implementation
    Base.@propagate_inbounds function only(x)
        i = iterate(x)
        @boundscheck if i === nothing
            throw(ArgumentError("Collection is empty, must contain exactly 1 element"))
        end
        (ret, state) = i
        @boundscheck if iterate(x, state) !== nothing
            throw(ArgumentError("Collection has multiple elements, must contain exactly 1 element"))
        end
        return ret
    end
end
d3 = Dictionary([:a, :b], [[1,2], [3,4]])
@test_throws ArgumentError only.(d3)


text:
using Cairo
using Rsvg

abstract type AbstractGraphics <: AbstractResource end

mutable struct Texture{T} <: AbstractGraphics
    ptr::Ptr{T}
    width::Int
    height::Int
    center_x::Int
    center_y::Int
end

function Texture(render_ptr::Ptr{SDL.Renderer}, sdl_surface::Ptr{SDL.Surface}; halign = 0.5, valign = 0.5)
    texture_ptr = SDL.CreateTextureFromSurface(render_ptr, sdl_surface)
    SDL.FreeSurface(sdl_surface)
    width, height = Int[0], Int[0]
    SDL.QueryTexture(texture_ptr, C_NULL, C_NULL, width, height)
    self = Texture(texture_ptr, width[], height[], round(Int, width[]*halign), round(Int, height[]*valign))
    finalizer(destroy!, self)
    return self
end

function Texture(resources::Resources, filename::String; width::Int = -1, height::Int = -1, scale::Real = 1.0)
    if filename in keys(resources)
        @debug("resource already loaded: '$filename'")
        return resources[filename]::Texture
    end

    ptr, width, height = load(query(filename), resources, width, height, scale)

    self = Texture(ptr, width, height, round(Int, width/2), round(Int, height/2))
    finalizer(destroy!, self)
    resources[filename] = self
    return self
end

function load(f::File{format"SVG"}, resources::Resources, width::Int, height::Int, scale::Real)
    rsvg_handle = Rsvg.handle_new_from_file(f.filename)
    Int(rsvg_handle.ptr) == 0 && error("'$(f.filename)' is not a valid SVG file")

    rsvg_dim = Rsvg.handle_get_dimensions(rsvg_handle)
    width = ceil(Int, scale*(width >= 0 ? width : rsvg_dim.width))
    height = ceil(Int, scale*(height >= 0 ? height : rsvg_dim.height))

    cairo_surface = Cairo.CairoImageSurface(fill(UInt32(0), (height, width)), Cairo.FORMAT_ARGB32)
    cairo_context = Cairo.CairoContext(cairo_surface)
    scale_x, scale_y = width / rsvg_dim.width, height / rsvg_dim.height
    Cairo.scale(cairo_context, scale_x, scale_y)
    Rsvg.handle_render_cairo(cairo_context, rsvg_handle)
    Cairo.destroy(cairo_context)

    width, height = ceil(Int, cairo_surface.width), ceil(Int, cairo_surface.height)
    sdl_surface = SDL.CreateRGBSurfaceFrom(pointer(cairo_surface.data),
                                           width, height, 32, Int32(width*4),
                                           0x00_ff_00_00, 0x00_00_ff_00, 0x00_00_00_ff, 0xff_00_00_00)
    texture_ptr = SDL.CreateTextureFromSurface(resources.render_ptr, sdl_surface)
    SDL.FreeSurface(sdl_surface)
    Cairo.destroy(cairo_surface)

    if texture_ptr == C_NULL
        error("Failed to load texture from: '$(f.filename)'")
    end

    return texture_ptr, width, height
end

function destroy!(texture::Texture{SDL.Texture})
    SDL.DestroyTexture(texture.ptr)
    texture.ptr = C_NULL
    return nothing
end


text:
using Bingomatic
using Test

@testset "Bingomatic.jl" begin
    @test_throws ErrorException sample_words(["a","b"], 2)
    @test_throws ErrorException sample_words(["a", "b", "c"], 4)
    n_rows = 3
    word_list = ["a","b","c","d","e","f","g","h","i"]
    words = sample_words(word_list, n_rows, n_rows)
    @test size(words) == (n_rows,n_rows)
    @test words[5] == "Free Space"
    card = make_card(words)
end


text:
using LinearAlgebra
using LoopVectorization

function RCCSDpT(rhf::RHF, Alg::abc)
    aoints = IntegralHelper()
    moints = IntegralHelper(orbitals=rhf.orbitals)
    ccsd = RCCSD(moints, aoints)
    return RCCSDpT(ccsd, moints, Alg)
end

function RCCSDpT(Alg::abc)
    aoints = IntegralHelper{Float64}()
    rhf = RHF(aoints)
    moints = IntegralHelper(orbitals=rhf.orbitals)
    ccsd = RCCSD(moints, aoints)
    return RCCSDpT(ccsd, moints, Alg)
end

function RCCSDpT(ccsd::RCCSD, moints::IntegralHelper{T,E,O}, Alg::abc) where {T<:AbstractFloat, E<:AbstractERI, O<:AbstractRestrictedOrbitals}

    output("\n   • Perturbative Triples Started\n")
    output("   - Using abc algorithm")

    T1 = ccsd.T1.data
    T2 = ccsd.T2.data

    # Use physicists notation
    ovvv = permutedims(moints["OVVV"].data, (1,3,2,4))
    ooov = permutedims(moints["OOOV"].data, (1,3,2,4))
    oovv = permutedims(moints["OVOV"].data, (1,3,2,4))

    o,v = size(T1)

    fo = moints["Fii"].data
    fv = moints["Faa"].data

    # Pre-allocate Intermediate arrays
    Ws  = [Array{T}(undef, o,o,o) for _ = 1:Threads.nthreads()]
    Vs  = [Array{T}(undef, o,o,o) for _ = 1:Threads.nthreads()]
    Evals = zeros(T, Threads.nthreads())

    output("Computing energy contribution from occupied orbitals:")
    BLAS_THREADS = BLAS.get_num_threads()
    BLAS.set_num_threads(1)
    TBLIS_THREADS = TBLIS.set_num_threads(3)

    t = @elapsed begin
    Threads.@threads for a in 1:v
    @inbounds begin
        id = Threads.threadid()
        W = Ws[id]
        V = Vs[id]

        # T1 views → V
        @views T1_a   = T1[:,a]

        # T2 views → W
        @views T2_a   = T2[:,:,:,a]

        # ooov views
        @views ooov_a = ooov[:,:,:,a]

        for b in 1:a
            Dab = fv[a] + fv[b]

            # T1 views → V
            @views T1_b  = T1[:,b]

            # T2 views → W
            @views T2_b  = T2[:,:,:,b]
            @views T2_ab = T2[:,:,a,b]
            @views T2_ba = T2[:,:,b,a]

            # ooov views
            @views ooov_b = ooov[:,:,:,b]

            # oovv views
            @views oovv_ab = oovv[:,:,a,b]
            @views oovv_ba = oovv[:,:,b,a]

            # ovvv views
            @views ovvv_ab = ovvv[:,:,a,b]
            @views ovvv_ba = ovvv[:,:,b,a]

            δab = a == b
            for c in 1:b
                Dabc = Dab + fv[c]

                # T1 views → V
                @views T1_c = T1[:,c]

                # T2 views → W
                @views T2_c = T2[:,:,:,c]
                @views T2_ac = T2[:,:,a,c]
                @views T2_ca = T2[:,:,c,a]
                @views T2_bc = T2[:,:,b,c]
                @views T2_cb = T2[:,:,c,b]

                # ooov views
                @views ooov_c = ooov[:,:,:,c]

                # oovv views
                @views oovv_ac = oovv[:,:,a,c]
                @views oovv_ca = oovv[:,:,c,a]
                @views oovv_bc = oovv[:,:,b,c]
                @views oovv_cb = oovv[:,:,c,b]

                # ovvv views
                @views ovvv_ac = ovvv[:,:,a,c]
                @views ovvv_ca = ovvv[:,:,c,a]
                @views ovvv_bc = ovvv[:,:,b,c]
                @views ovvv_cb = ovvv[:,:,c,b]

                # Build intermediates W and V for the set abc
                @tensor begin
                    W[i,j,k]  = ovvv_ab[i,d]*T2_c[j,k,d] - ooov_c[j,k,l]*T2_ab[i,l] #(ijk)(abc)
                    W[i,j,k] += ovvv_ac[i,d]*T2_b[k,j,d] - ooov_b[k,j,l]*T2_ac[i,l] #(ikj)(acb)

                    W[i,j,k] += ovvv_ba[j,d]*T2_c[i,k,d] - ooov_c[i,k,l]*T2_ba[j,l] #(jik)(bac)
                    W[i,j,k] += ovvv_bc[j,d]*T2_c[k,i,d] - ooov_a[k,i,l]*T2_bc[j,l] #(jki)(bca)

                    W[i,j,k] += ovvv_ca[k,d]*T2_b[i,j,d] - ooov_b[i,j,l]*T2_ca[k,l] #(kij)(cab)
                    W[i,j,k] += ovvv_cb[k,d]*T2_c[j,i,d] - ooov_a[j,i,l]*T2_cb[k,l] #(kji)(cba)

                    V[i,j,k] = W[i,j,k] + oovv_bc[j,k]*T1_a[i] + oovv_ac[i,k]*T1_b[j] + oovv_ab[i,j]*T1_c[k]
                end

                # Compute Energy contribution
                δbc = b === c
                for i in 1:o
                    Diabc = fo[i] - Dabc
                    for j in 1:i
                        Dijabc = fo[j] + Diabc
                        δij = i === j
                        @fastmath for k in 1:j
                            Dd = Dijabc + fo[k]
                            δjk = j === k
                            X = (W[i,j,k]*V[i,j,k] + W[i,k,j]*V[i,k,j] + W[j,i,k]*V[j,i,k] + W[j,k,i]*V[j,k,i] + W[k,i,j]*V[k,i,j] + W[k,j,i]*V[k,j,i])
                            Y = (V[i,j,k] + V[j,k,i] + V[k,i,j])
                            Z = (V[i,k,j] + V[j,i,k] + V[k,j,i])
                            Ef  = (Y - 2*Z)*(W[i,j,k] + W[k,i,j] + W[j,k,i])
                            Ef += (Z - 2*Y)*(W[i,k,j] + W[j,i,k] + W[k,j,i]) + 3*X
                            Evals[id] += 2*Ef / (Dd * (1+ δij + δjk)*(1 + δab + δbc))
                        end
                    end
                end
            end
        end
        #output("  Orbital {} ✔", a)
    end # Threads (i loop)
    end # inbounds
    end # time

    BLAS.set_num_threads(BLAS_THREADS)
    TBLIS.set_num_threads(TBLIS_THREADS)
    Et = sum(Evals)
    output("Finished in {:5.5f} s", t)
    output("Final (T) contribution: {:15.10f}", Et)
    output("CCSD(T) energy:         {:15.10f}", Et+ccsd.energy)
    return RCCSDpT{T}(ccsd, Et+ccsd.energy, Et)
end


text:
using NPZ
using CSV
using Plots
pyplot()

simdata = CSV.read(joinpath(@__DIR__, "benchmark_data.csv"))
julia_times = CSV.read(joinpath(@__DIR__, "julia_times.csv"))
python_times = CSV.read(joinpath(@__DIR__, "python_times.csv"))

# Global times
jtime = julia_times.time
ptime = python_times.time
ratios = ptime ./ jtime

# @layout [ a b ;
#            c   ]

histogram(jtime, alpha=0.3)

histogram(ptime, alpha=0.3)

using Statistics

scatter(ratios, yscale=:log10, markerstrokealpha=0, label="")
m = mean(ratios)
hline!([m], c=1, label="μ=$(round(m, digits=1))")
hline!([1], label="", c=:black, ls=:dash)
ylabel!("Time Ratio (Python/Julia)")
xlabel!("Run Number")

# histogram(ratios, alpha=0.3, xscale=:log10, xlabel="Time ratio in seconds (Python/Julia)")

max_stars_repo_name: zclaytor/butterpy

max_issues_repo_name: zclaytor/butterpy

max_forks_repo_name: zclaytor/butterpy


text:
using WinRPM
using Test

@testset "WinRPM.jl" begin
    @testset "WinRPM.update" begin
        for (ignorecache, allow_remote) in [(false, true), (true, true)]
            WinRPM.update(ignorecache, allow_remote)
        end
    end

    if Sys.iswindows() && lowercase(get(ENV, "CI", "")) == "true"
        @testset "WinRPM.install" begin
            WinRPM.update()
            pkg = "bzip2"
            @info("Attempting to install $(pkg)")
            WinRPM.install(pkg; yes = true)
            @test isfile(WinRPM.installedlist)
            installedlist = read(WinRPM.installedlist, String)
            @test occursin(pkg, installedlist)
            @info("Successfully installed $(pkg)")
        end
    end
end


text:
abstract type MeasureOp{n} <: Measure end
abstract type UnaryOp{A} <: MeasureOp{1} end
abstract type ScalarOp{A} <: MeasureOp{2} end
abstract type BinaryOp{A, B} <: MeasureOp{2} end

struct Neg{A <: Measure} <: UnaryOp{A}
    a::A
end

struct Add{A <: Measure, B <: Measure} <: BinaryOp{A, B}
    a::A
    b::B
end

struct Min{A <: Measure, B <: Measure} <: BinaryOp{A, B}
    a::A
    b::B
end

struct Max{A <: Measure, B <: Measure} <: BinaryOp{A, B}
    a::A
    b::B
end

struct Div{A <: Measure} <: ScalarOp{A}
    a::A
    b::Number
end

struct Mul{A <: Measure} <: ScalarOp{A}
    a::A
    b::Number
end

# Easy simplifications
# TODO: Add more simplifications
Neg(x::Neg) = x.a
iszero(x::Measure) = false

Base.:+(a::Measure, b::Measure) = Add(a, b)
Base.:-(a::Measure) = Neg(a)
Base.:-(a::Neg) = a.value
Base.:-(a::Measure, b::Measure) = Add(a, -b)
Base.:/(a::Measure, b::Number) = Div(a, b)
Base.:/(a::T, b::T) where T <: Measure = Div(a, b)
Base.:*(a::Measure, b::Number) = Mul(a, b)
Base.:*(a::Number, b::Measure) = Mul(b, a)
Base.min(a::Measure, b::Measure) = Min(a, b)
Base.max(a::Measure, b::Measure) = Max(a, b)


text:
using Documenter
using GRIB

makedocs(
    sitename = "GRIB.jl Documentation",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    modules = [GRIB],
    authors = "Alex Weech and contributors",
    pages = [
        "Home" => "index.md",
        "Manual" => ["gribfile.md", "message.md", "indexer.md", "nearest.md"]
    ]
)

deploydocs(
    repo = "github.com/weech/GRIB.jl.git"
)


text:
"""
    get_xsep()

Return x-position of the separation point closest to the trailing edge on the
upper and lower airfoil surfaces. This location is defined as the point where
the friction coefficient becomes negative.
"""
get_xsep

"""
    get_xsep_cs()

`get_xsep` for the complex step enabled version of XFOIL.
"""
get_xsep_cs

# definition of get_xsep
for (T, name, globals, bldump) in
    ((:Float64, :get_xsep, :xfoilglobals, :bldump),
    (:ComplexF64, :get_xsep_cs, :xfoilglobals_cs, :bldump_cs))

    @eval begin

        function $(name)(; xtol=1e-6)

            s, x, y, ue, dstar, theta, cf = $(bldump)()

            # Find indices of LE and TE edge
            xile = findfirst(x -> abs(x - $(globals).xle[1]) <= xtol, x)
            xite1 = findfirst(x -> abs(x - $(globals).xte[1]) <= xtol, x)
            xite2 = findnext(x -> abs(x - $(globals).xte[1]) <= xtol, x, xite1+1)

            # Get both faces going from TE to LE
            iup = xite1:(-1)^(xite1 > xile):xile
            ilo = xite2:(-1)^(xite2 > xile):xile

            # xup, yup = x[iup], y[iup]
            # xlo, ylo = x[ilo], y[ilo]

            # Find first separation point from the TE to LE on each side
            iupsep = findfirst( val -> val >= 0, cf[iup])
            ilosep = findfirst( val -> val >= 0, cf[ilo])

            iupsep == nothing ? iupsep = length(iup) : nothing
            ilosep == nothing ? ilosep = length(ilo) : nothing

            xsepup = x[iup[iupsep]]
            xseplo = x[ilo[ilosep]]

            return xsepup, xseplo
        end

    end

end

"""
    alpha_sweep(x, y, alpha, re; kwargs...)

Perform angle of attack sweep using XFOIL.  Return cl, cd, cdp, cm, converged.

# Arguments
- `x`: Airfoil coordinates start from trailing edge looping counterclockwise
- `y`:
- `alpha`: Array of angle of attacks in degrees
- `re`: Reynolds number
- `ma`: Mach number
- `iter=50`: Maximum number of iterations for each XFOIL run
- `npan=140`: Number of panels
- `percussive_maintenance=true`: Flag to call `do_percussive_maintenance` upon convergence failure
- `printdata=false`: Flag to indicate whether to print data obtained from XFOIL during the solution
- `zeroinit=true`: Flag to indicate whether to start angle of attack sweeps from zero angle of attack
- `clmaxstop=false`: Stop if lift coefficient decreases twice consecutively going up
- `clminstop=false`: Stop if lift coefficient increases twice consecutively going down
"""
alpha_sweep

"""
    alpha_sweep_cs(x, y, alpha, re; kwargs...)

`alpha_sweep` for the complex step enabled version of XFOIL.
"""
alpha_sweep_cs

# definition of alpha_sweep
for (T, name, set_coordinates, pane, solve_alpha, do_percussive_maintenance) in
    ((:Float64, :alpha_sweep, :set_coordinates, :pane, :solve_alpha, :do_percussive_maintenance),
    (:ComplexF64, :alpha_sweep_cs, :set_coordinates_cs, :pane_cs, :solve_alpha_cs, :do_percussive_maintenance_cs))

    @eval begin

        function $(name)(x, y, alpha, re; mach=0.0, iter=50, npan=140,
            percussive_maintenance=true, printdata=false, zeroinit=true,
            clmaxstop=false, clminstop=false, ncrit=9)

            @assert length(x) == length(y) "x and y arrays must have the same length"

            naoa = length(alpha)

            if zeroinit
                # Set up angle of attack range going up and going down from zero if specified
                # This helps XFOIL to converge more consistently

                # separate out negative angles of attack
                aoaneg = sort(alpha[findall(real.(alpha) .< 0.0)], by=real, rev = true)

                # add zero angle of attack to negative angles of attack
                pushfirst!(aoaneg, 0.0)

                # perform angle of attack sweep for negative angles of attack
                clneg, cdneg, cdpneg, cmneg, convneg = $(name)(x,
                    y, aoaneg, re, mach, iter, npan, percussive_maintenance, printdata,
                    false, clminstop, ncrit)

                # separate out positive angles of attack
                aoapos = sort(alpha[findall(real.(alpha) .>= 0.0)], by=real)

                # add zero angle of attack to positive angles of attack
                pushfirst!(aoapos,0.0)

                # perform angle of attack sweep for positive angles of attack
                clpos, cdpos, cdppos, cmpos, convpos = $(name)(x,
                    y, aoapos, re, mach, iter, npan, percussive_maintenance, printdata,
                    clmaxstop, false, ncrit)

                # combine results from negative and positive runs (excluding zero angle of attack runs)
                cl = vcat(clneg[end:-1:2], clpos[2:end])
                cd = vcat(cdneg[end:-1:2], cdpos[2:end])
                cdp = vcat(cdpneg[end:-1:2], cdppos[2:end])
                cm = vcat(cmneg[end:-1:2], cmpos[2:end])
                conv = vcat(convneg[end:-1:2], convpos[2:end])
            else
                cl, cd, cdp, cm, conv = $(name)(x, y, alpha, re,
                    mach, iter, npan, percussive_maintenance, printdata, clminstop,
                    clmaxstop, ncrit)
            end

            return cl, cd, cdp, cm, conv
        end

    end

    @eval begin

        function $(name)(x, y, alpha, re, mach, iter, npan, percussive_maintenance,
            printdata, clmaxstop, clminstop, ncrit)

            # Set up storage arrays
            naoa = length(alpha)
            cl = zeros($T, naoa)
            cd = zeros($T, naoa)
            cdp = zeros($T, naoa)
            cm = zeros($T, naoa)
            converged = zeros(Bool, naoa)

            # print header for data #TODO: add option to print to file rather than terminal
            if printdata == true
                println("\nAngle\t\tCl\t\tCd\t\tCm\t\tConverged")
            end

            # start unconverged
            for i in eachindex(alpha)

                # set coordinates if its the first iteration, also set the coordinates
                # if convergence has failed since it resets XFOIL's initial guess
                # TODO: reset XFOIL's initial guess without setting coordinates again
                if i == 1 || converged[i-1] == false
                    $(set_coordinates)(x, y)
                    $(pane)(npan = npan)
                    # TODO: remove hard-coded repaneling, add it as a default option
                end

                # run XFOIL
                cl[i], cd[i], cdp[i], cm[i], converged[i] = $(solve_alpha)(alpha[i], re, mach=mach, iter=iter, ncrit=ncrit)

                # check convergence, do percussive maintenance if desired and necessary
                if !converged[i] && percussive_maintenance
                    cl[i], cd[i], cdp[i], cm[i], converged[i] = $(do_percussive_maintenance)(x, y, alpha[i], re, mach, iter, npan, ncrit)
                end

                # print data from the run #TODO: add option to print to file rather than terminal
                if printdata == true
                    @printf("%8f\t%8f\t%8f\t%8f\t%d\n", real(alpha[i]), real(cl[i]), real(cd[i]), real(cm[i]), converged[i])
                end

                aoaconv = alpha[findall(converged)]
                clconv = cl[findall(converged)]
                if (real(alpha[i]) > 0.0) && (length(clconv) >= 4) # TODO: add minimum number of angle of attacks (currently 4) as an input option
                    if clmaxstop #TODO: allow clmaxstop to trigger after user specified number of decreasing lift runs
                        # break if maximum found #TODO: allow user specified function definition for what clmax is
                        if (real(clconv[end]) < real(clconv[end-1])) && (real(clconv[end-1]) < real(clconv[end-2]))
                            break
                        end
                    end
                end
                if (real(alpha[i]) < 0.0) && (length(clconv) >= 4) # TODO: add minimum number of angles of attacks (currently 4) as an input option
                    if clminstop # TODO: allow clminstop to trigger after user specified number of increaseing lift runs
                        # break if minimum found #TODO: allow user specified function definition for what clmin is
                        if (real(clconv[end]) > real(clconv[end-1])) && (real(clconv[end-1]) > real(clconv[end-2]))
                            break
                        end
                    end
                end
            end

            return cl, cd, cdp, cm, converged
        end

    end

end

"""
    do_percussive_maintenance(x, y, alpha, re, iter, npan, ncrit)

Attempts to converge previously unconverged XFOIL solutions by perturbing the
starting point for the XFOIL solve. Returns cl, cd, cdp, cm, converged
"""
do_percussive_maintenance

"""
    do_percussive_maintenance_cs(x, y, alpha, re, iter, npan, ncrit)

`do_percussive_maintenance` for the complex step enabled version of XFOIL.
"""
do_percussive_maintenance_cs

# definition of do_percussive_maintenance
for (T, name, set_coordinates, pane, solve_alpha) in
    ((:Float64, :do_percussive_maintenance, :set_coordinates, :pane, :solve_alpha),
    (:ComplexF64, :do_percussive_maintenance_cs, :set_coordinates_cs, :pane_cs, :solve_alpha_cs))

    @eval begin

        function $(name)(x, y, alpha, re, mach, iter, npan, ncrit)

            # set new parameters to original parameters
            remod = re
            alphamod = alpha

            f = (alphamod, remod) -> $(solve_alpha)(alphamod, remod, mach=mach, iter=iter, ncrit=ncrit)

            # TODO: add number of iterations of percussive maintenance as an input parameter
            # perhaps by adding a PercussiveMaintenance type with customizable parameters
            for j = 1:25
                # perturb new parameters
                remod = remod+1000
                alphamod = 0.95*alphamod

                # reinput coordinates in order to re-initialize the XFOIL solution
                # TODO: allow reinitializing XFOIL without re-inputting coordinates
                # TODO: change hard-coded pane command to an option
                $(set_coordinates)(x, y)
                $(pane)(npan=npan)

                # try modified parameters
                cl, cd, cdp, cm, converged = f(alphamod, remod)

                # exit loop if converged
                if converged
                    break
                end
            end

            # Gradually bring the parameters back to those of the original problem
            if converged
                # get difference between actual parameters and modified parameters
                rediff = remod - re
                aoadiff = alphamod-alpha
                # over ten iterations, gradually return to the original problem
                # TODO: add number of iterations taken to return to the original problem
                # as a customizable parameter
                for j = 1:9
                    remod = remod - rediff/10
                    alphamod = alphamod - aoadiff/10
                    cl, cd, cdp, cm, converged = f(alphamod, remod)
                end
                # now try solving the original problem again
                cl, cd, cdp, cm, converged = f(alpha, re)
            end

            return cl, cd, cdp, cm, converged
        end

    end

end


text:
"""
    cache(results::Vector{T}) where {T <: Union{MangalReferenceTaxon,MangalNode,MangalNetwork}}

Internally, the `Mangal` package uses a cache to store some objects that are
likely to be queried more than once. These are `MangalNode` and
`MangalReferenceTaxon`, which are called in a nested way during the querying of
*e.g.* `Interactions`. This is **not** a fancy mechanism, and it only works when
calling the nodes or backbones by their `id` (which is what the resources-hungry
functions do internally anyways).
"""
function cache(results::Vector{T}) where {T <: Union{MangalReferenceTaxon,MangalNode,MangalNetwork,MangalAttribute}}
    for result in results
        if !haskey(_MANGAL_CACHES[T], result.id)
            global _MANGAL_CACHES[T][result.id] = result
        end
    end
end

"""
    generate_base_header()

If a bearer token is present, this function will add it to the header.
"""
function generate_base_header()
    headers = ["Content-Type" => "application/json"]
    if haskey(ENV, "MANGAL_BEARER_TOKEN")
        push!(headers, "Authorization" => "bearer $(ENV["MANGAL_BEARER_TOKEN"])")
    end
    return headers
end

"""
    generate_request_query(parameters::Pair...)

Takes a series of `Pairs`, and returns an URL-ready query string.
"""
function generate_request_query(parameters::Pair...)
    query = ""
    for (i, pair) in enumerate(promote(parameters...))
        i == 1 || (query *= "&")
        i == 1 && (query *= "?")
        query *= pair.first * "=" * string(pair.second)
    end
    return replace(query, " " => "+")
end

"""
    search_objects_by_query(endpoint::AbstractString, formatter::Function, query::Pair...)

In all cases, it is assumed that the functions will be wrapepd
in calls to query objects until no further objects are found.
"""
function search_objects_by_query(ReturnType::Type, query::Pair...)
    # Headers
    headers = Mangal.generate_base_header()

    # Full endpoint
    endpoint = Mangal.api_root * Mangal._MANGAL_ENDPOINTS[ReturnType]

    # Convert query parameters
    request_url = length(query) == 0 ? endpoint : endpoint*Mangal.generate_request_query(query...)

    # Perform the request
    this_request = HTTP.get(request_url, headers)
    request_body = String(this_request.body)
    request_body = replace(request_body, "" => "-") # This shouldn't be necessary anymore but...

    # Returns the collection
    parsed_json = JSON.parse.(request_body)

    # Return the formatted object(s)
    formatter = (x) -> format_mangal_response(ReturnType, x)
    return convert(Vector{ReturnType}, formatter.(parsed_json))
end

"""
    number_of_objects(endpoint::AbstractString, query::Pair...)

This function will return the number of objects already queried, and remaining
to query, for any given query.
"""
function number_of_objects(endpoint::AbstractString, query::Pair...)
    # Headers
    headers = Mangal.generate_base_header()

    # Full endpoint
    endpoint = Mangal.api_root * endpoint

    # Convert query parameters
    request_url = query == nothing ? endpoint : endpoint*Mangal.generate_request_query("count" => 1, query...)

    # Perform the request
    this_request = HTTP.get(request_url, headers)

    content_range = first(filter(head -> head.first == "Content-Range", this_request.headers))
    positions = parse.(Int64, split(content_range.second, [' ', '-', '/'])[2:end])
    range_begin = positions[1]+1
    range_stop = positions[2]+1
    range_ends = positions[3]

    return range_ends
end


text:
maxdims["dmrg_2d_conserve_ky"] = 1_000:1_000:10_000
descriptions["dmrg_2d_conserve_ky"] = "DMRG, 2D Hubbard model (U/t = 4)\n(Nx, Ny) = (8, 4), 10 sweeps\nHybrid real and momentum space"

function runbenchmark(::Val{:dmrg_2d_conserve_ky};
                      maxdim::Int, nsweeps::Int = 10,
                      outputlevel::Int = 0, conserve_ky::Bool = true,
                      cutoff::Float64 = 0.0, Nx::Int = 8, Ny::Int = 4,
                      U::Float64 = 4.0, t::Float64 = 1.0,
                      splitblocks::Bool = false)
  N = Nx * Ny
  sweeps = Sweeps(nsweeps)
  maxdims = min.(maxdim, [100, 200, 400, 800, 2000, 10_000, maxdim])
  maxdim!(sweeps, maxdims...)
  cutoff!(sweeps, cutoff)
  noise!(sweeps, 1e-6, 1e-7, 1e-8, 0.0)
  sites = siteinds("ElecK", N;
                   conserve_qns = true,
                   conserve_ky = conserve_ky,
                   modulus_ky = Ny)
  ampo = hubbard(Nx = Nx, Ny = Ny, t = t, U = U, ky = true)
  H = MPO(ampo, sites)
  if splitblocks
    H = ITensors.splitblocks(linkinds, H)
  end
  # Create start state
  state = Vector{String}(undef, N)
  for i in 1:N
    x = (i - 1) ÷ Ny
    y = (i - 1) % Ny
    if x % 2 == 0
      if y % 2 == 0
        state[i] = "Up"
      else
        state[i] = "Dn"
      end
    else
      if y % 2 == 0
        state[i] = "Dn"
      else
        state[i] = "Up"
      end
    end
  end
  ψ0 = productMPS(sites, state)
  energy, ψ = dmrg(H, ψ0, sweeps; outputlevel = outputlevel)
  if outputlevel > 0
    @show nsweeps
    @show Nx, Ny
    @show U, t
    @show cutoff
    @show energy
    @show flux(ψ)
  end
  return maxlinkdim(ψ)
end


text:
"""
    pauli_are_commuting(operators) -> Bool

Return true if the Pauli operators mutually commute, or false otherwise, where `operators`
is an iterable of `AbstractVector{Int}`.

# Examples
```jldoctest
julia> stabilizers = [[1, 3, 3, 1, 0], [0, 1, 3, 3, 1], [1, 0, 1, 3, 3], [3, 1, 0, 1, 3]];

julia> pauli_are_commuting(stabilizers)  # stabilizers mutually commute
true

julia> logical_x, logical_z = [1, 1, 1, 1, 1], [3, 3, 3, 3, 3];  # XXXXX, ZZZZZ

julia> pauli_are_commuting([logical_x, logical_z])  # logicals do not commute
false
```
"""
function pauli_are_commuting(operators)
    return all(==(0), pauli_commutation(a, b) for (a, b) in combinations(operators, 2))
end

"""
    pauli_are_independent(operators) -> Bool

Return true if the Pauli operators are independent, or false otherwise, where
`operators` is an iterable of `AbstractVector{Int}`.  Here independence means
that no operator can be expressed as a product of the others.

# Examples
```jldoctest
julia> pauli_are_independent([[3, 3, 0], [0, 3, 3]])  # ZZI, IZZ
true

julia> pauli_are_independent([[3, 3, 0], [0, 3, 3], [3, 0, 3]])  # ZZI, IZZ, ZIZ
false
```
"""
function pauli_are_independent(operators)
    operator_list = collect(operators)  # collect iterable to use length and index access

    num_operators = length(operator_list)
    if num_operators == 0  # empty operator_list is independent
        return true
    end
    num_qubits = length(operator_list[1])
    remaining = collect(1:num_operators)

    for qubit in 1:num_qubits
        paulis = [1,2,3]
        indices = Int[]
        for α in remaining
            if operator_list[α][qubit] in paulis
                setdiff!(paulis, operator_list[α][qubit])
                push!(indices, α)
            end
            if length(paulis) == 1
                break
            end
        end

        remaining = setdiff(remaining, indices)
        for α in remaining
            if operator_list[α][qubit] == 0
                continue
            end
            new_operator = pauli_product.(operator_list[α], operator_list[indices[1]])
            if new_operator[qubit] == 0
                operator_list[α] = deepcopy(new_operator)
                continue
            end
            if length(paulis) == 2
                continue
            end
            new_operator = pauli_product.(operator_list[α], operator_list[indices[2]])
            if new_operator[qubit] == 0
                operator_list[α] = deepcopy(new_operator)
                continue
            end
            new_operator = pauli_product.(operator_list[α], operator_list[indices[1]])
            new_operator = pauli_product.(new_operator, operator_list[indices[2]])
            if new_operator[qubit] == 0
                operator_list[α] = deepcopy(new_operator)
                continue
            end
        end

        for α in remaining
            if operator_list[α] == zeros(Int64, num_qubits)
                return false
            end
        end
    end

    return true
end

"""
    pauli_commutation(a::Int, b::Int) -> 0 or 1
    pauli_commutation(a::AbstractVector{Int}, b::AbstractVector{Int}) -> 0 or 1

Return the commutation relation of two Paulis, i.e. 0 if they commute and 1 if not.

# Examples
```jldoctest
julia> pauli_commutation(1, 1)  # X commutes with itself
0

julia> pauli_commutation(1, 3)  # X does not commute with Z
1
```

```jldoctest
julia> pauli_commutation([1, 3, 3, 1, 0], [1, 1, 1, 1, 1])  # XZZXI commutes with XXXXX
0

julia> stabilizers = [[1, 3, 3, 1, 0], [0, 1, 3, 3, 1], [1, 0, 1, 3, 3], [3, 1, 0, 1, 3]];

julia> error = [1, 1, 0, 0, 0];  # XXIII

julia> syndrome = pauli_commutation.(stabilizers, Ref(error))
4-element Vector{Int64}:
 1
 0
 0
 1
```
"""
function pauli_commutation(a::Int, b::Int)
    return a == 0 || b == 0 || a == b ? 0 : 1
end
function pauli_commutation(a::AbstractVector{Int}, b::AbstractVector{Int})
    return mod(sum(pauli_commutation(p1, p2) for (p1, p2) in zip(a, b); init=0), 2)
end

"""
    pauli_pow(a::Int, power::Int) -> Int

Return the Pauli raised to the given power.

# Examples
```jldoctest
julia> x_squared = pauli_pow(1, 2)
0

julia> z_inverse = pauli_pow(3, -1)
3
```
"""
function pauli_pow(a::Int, power::Int)
    return mod(power, 2) == 0 ? 0 : a
end

"""
    pauli_product(a::Int, b::Int) -> Int

Return the product of two Paulis.  Note that we are ignoring overall minus
signs.

# Examples
```jldoctest
julia> pauli_product(1, 2) # X.Y -> Z
3

julia> pauli_product.([1, 0, 3, 2], [1, 1, 1, 3]) # (XIZY).(XXXZ) -> (IXYX)
4-element Vector{Int64}:
 0
 1
 2
 1
```
"""
function pauli_product(a::Int, b::Int)
    if a == b  # pauli squared is identity
        return 0
    elseif a == 0 || b == 0  # pauli times identity is pauli
        return a + b
    else  # product of distinct paulis is the other pauli
        # X.Y=Z -> 6-(1+2)=3,  X.Z=Y -> 6-(1+3)=2,  Y.Z=X -> 6-(2+3)=1
        return 6 - (a + b)
    end
end

"""
    pauli_product(operators) -> Vector{Int}

Return the product of the Pauli operators, where `operators` is an iterable of
`AbstractVector{Int}`.

See also [`pauli_product_pow`](@ref).

# Examples
```jldoctest
julia> ops = [[3, 3, 0], [0, 1, 1], [2, 0, 2]];  # ZZI, IXX, YIY

julia> pauli_product(ops)
3-element Vector{Int64}:
 1
 2
 3
```
"""
function pauli_product(operators)
    return reduce((a, b)->pauli_product.(a, b), operators)
end

"""
    pauli_product_pow(operators, powers) -> Vector{Int}

Return the product of the Pauli operators each raised to the corresponding power, where
`operators` is an iterable of `AbstractVector{Int}` and `powers` is an iterable of `Int`.

Note: The product is evaluated to the length of the shorter of the two iterables.

See also [`pauli_product`](@ref), [`pauli_pow`](@ref).

# Examples
```jldoctest
julia> ops = [[3, 3, 0], [0, 1, 1], [2, 0, 2]];  # ZZI, IXX, YIY

julia> pauli_product_pow(ops, [1, 0, 1])
3-element Vector{Int64}:
 1
 3
 2
```
"""
function pauli_product_pow(operators, powers)
    return reduce(
        (a, b)->pauli_product.(a, b),
        (pauli_pow.(operator, power) for (operator, power) in zip(operators, powers))
    )
end

"""
    pauli_random_operator(n::Int, p::Real, rng::AbstractRNG=GLOBAL_RNG)
        -> AbstractVector{Int}

Return a random Pauli operator on `n` qubits with Paulis applied according to i.i.d.
depolarizing noise with probability `p`.

# Examples
```jldoctest
julia> using Random:MersenneTwister  # use RNG for reproducible example

julia> pauli_random_operator(5, 0.2, MersenneTwister(13))
5-element Vector{Int64}:
 3
 3
 0
 2
 0
```
"""
function pauli_random_operator(n::Int, p::Real, rng::AbstractRNG=GLOBAL_RNG)
    weights = ProbabilityWeights([1 - p, p / 3, p / 3, p / 3])
    return sample(rng, [0, 1, 2, 3], weights, n)
end

"""
    pauli_rep_change(pauli::Int) -> Char
    pauli_rep_change(pauli::Char) -> Int

Convert Pauli between alphabetical and numerical representation.

# Examples
```jldoctest
julia> pauli_rep_change(1)
'X': ASCII/Unicode U+0058 (category Lu: Letter, uppercase)

julia> pauli_rep_change('X')
1

julia> pauli_rep_change.((0, 1, 2, 3))
('I', 'X', 'Y', 'Z')

julia> pauli_rep_change.(('I', 'X', 'Y', 'Z'))
(0, 1, 2, 3)
```
"""
function pauli_rep_change(pauli::Int)
    return "IXYZ"[pauli+1]
end
const _pauli_rep_map = Dict(zip("IXYZ", 0:3))
function pauli_rep_change(pauli::Char)
    return _pauli_rep_map[pauli]
end

"""
    pauli_weight(operator::AbstractVector{Int}) -> Int

Return the weight of the Pauli operator, i.e. the number of non-identity elements.

# Examples
```jldoctest
julia> pauli_weight([3, 0, 2, 0, 1])  # ZIYIX
3
```
"""
function pauli_weight(operator::AbstractVector{Int})
    return count(!=(0), operator)
end


text:
"""
    bulkformulae(atemp,aqh,speed,sst,hu=10,ht=2,hq=2,zref=10,atmrho=1.2)

Units:
atemp  - mean air temperature (K)  at height ht (m)
aqh    - mean air humidity (kg/kg) at height hq (m)
speed  - mean wind speed (m/s)     at height hu (m)
sst    - sea surface temperature (K)

Bulk formulae formulation:
```
wind stress = (ust,vst) = rhoA * Cd * Ws * (del.u,del.v)
Sensib Heat flux = fsha = rhoA * Ch * Ws * del.T * CpAir
Latent Heat flux = flha = rhoA * Ce * Ws * del.Q * Lvap
                 = -Evap * Lvap
```

with Cd,Ch,Ce = drag coefficient, Stanton number and
Dalton number respectively [no-units], function of
height & stability; and

```
Ws = wind speed = sqrt(del.u^2 +del.v^2)
del.T = Tair - Tsurf ; del.Q = Qair - Qsurf
```
"""
function bulkformulae(atemp,aqh,speed,sst,hu=10,ht=2,hq=2,zref=10,atmrho=1.2)

umin=0.1;
karman=0.4;
gravity_mks=9.81;
cen2kel=273.150;
humid_fac=0.606;
cvapor_fac=640380;
#cvapor_fac_ice=11637800;
cvapor_exp=5107.4;
#cvapor_exp_ice=5897.8;
saltsat=0.980;
#atmrho=1.2;
gamma_blk=0.01;
cdrag_1=0.0027000;
cdrag_2=0.0001420;
cdrag_3=0.0000764;
cstanton_1 = 0.0327;
cstanton_2 = 0.0180;
cdalton = 0.0346;
niter_bulk=2;
psim_fac=5;
atmcp=1005;
flamb=2500000;
rhoConstFresh=999.8;
stefanBoltzmann = 5.670e-8;
ocean_emissivity=5.50e-8 / 5.670e-8;
albedo=0.1;

#-- Set surface parameters :
zwln = log(hu./zref);
ztln = log(ht./zref);
czol = hu*karman*gravity_mks;

wspeed=max(speed,umin);

#-   Surface Temp.
Tsf = sst + cen2kel;

#--- Compute turbulent surface fluxes
#-   Pot. Temp and saturated specific humidity

t0 = atemp.*(1 + humid_fac*aqh);
tmpbulk = cvapor_fac*exp(-cvapor_exp./Tsf);
ssq = saltsat*tmpbulk/atmrho;
deltap = atemp + gamma_blk*ht - Tsf;
delq   = aqh - ssq;

#--  initial guess for exchange coefficients:
#    take U_N = del.U ; stability from del.Theta ;
stable = 0.5 + 0.5*sign(deltap);

tmpbulk = cdrag_1./wspeed + cdrag_2 + cdrag_3.*wspeed; # CD

rdn = sqrt(tmpbulk);
rhn = (1-stable)*cstanton_1 + stable*cstanton_2;
ren = cdalton;
#--  calculate turbulent scales
ustar=rdn.*wspeed;
tstar=rhn.*deltap;
qstar=ren*delq;

rd=0;rh=0;re=0;huol=0;psimh=0;psixh=0;

#--- iterate with psi-functions to find transfer coefficients
for iter=1:niter_bulk

    huol = ( tstar./t0 + qstar./(1/humid_fac+aqh) ).*czol./(ustar.*ustar);

    #cph The following is different in Large&Pond1981 code:
    #cph huol = max(huol,zolmin) with zolmin = -100
    tmpbulk = min(abs(huol),10);
    huol   = tmpbulk.*sign(huol);
    htol   = huol.*ht./hu;
    hqol   = huol.*hq./hu;
    stable = 0.5 + 0.5*sign(huol);

    #                 Evaluate all stability functions assuming hq = ht.
    # The following is different in Large&Pond1981 code:
    # xsq = max(sqrt(abs(1 - 16.*huol)),1)
    xsq    = sqrt( abs(1 - huol*16) );
    x      = sqrt(xsq);
    psimh = -psim_fac*huol.*stable + (1-stable).*
        ( log( (1 + 2*x + xsq).*(1+xsq)*.125 ) - 2*atan(x) + 0.5*pi );

    # The following is different in Large&Pond1981 code:
    # xsq = max(sqrt(abs(1 - 16.*htol)),1)
    xsq   = sqrt( abs(1 - htol*16) );
    psixh = -psim_fac*htol.*stable + (1-stable).*( 2*log(0.5*(1+xsq)) );

    #-   Shift wind speed using old coefficient
    usn = wspeed./(1 + rdn/karman.*(zwln-psimh) );
    usm = max(usn, umin);

    #-   Update the 10m, neutral stability transfer coefficients
    tmpbulk = cdrag_1./usm + cdrag_2 + cdrag_3.*usm;

    rdn = sqrt(tmpbulk);
    rhn = (1-stable)*cstanton_1 + stable*cstanton_2;
    ren = cdalton;

    #-   Shift all coefficients to the measurement height and stability.
    rd = rdn./(1 + rdn.*(zwln-psimh)/karman);
    rh = rhn./(1 + rhn.*(ztln-psixh)/karman);
    re = ren./(1 + ren.*(ztln-psixh)/karman);

    #--  Update ustar, tstar, qstar using updated, shifted coefficients.
    ustar = rd.*wspeed;
    qstar = re.*delq;
    tstar = rh.*deltap;

    # end of iteration loop
end

#-   Coeff:
tau   = atmrho*rd.*wspeed;

#-   Turbulent Fluxes
hs = atmcp*tau.*tstar;
hl = flamb*tau.*qstar;
#   change sign and convert from kg/m^2/s to m/s via rhoConstFresh
evap = -tau.*qstar/rhoConstFresh;
ch=rd*rh
ce=rd*re

out = Dict("hl"=>hl,"hs"=>hs,"evap"=>evap,"ch"=>ch,"ce"=>ce,
"tau"=>tau,"ssq"=>ssq,"huol"=>huol,"rd"=>rd,"re"=>re,"rh"=>rh,
"ustar"=>ustar,"qstar"=>qstar,"tstar"=>tstar,"psimh"=>psimh,"psixh"=>psixh)

return out

end


text:
Base.isequal(x::NanoDate, y::NanoDate) =
    isequal(x.datetime, y.datetime) && isequal(x.nanosecs, y.nanosecs)

Base.isequal(x::NanoDate, y::Date) =
    isequal(x.datetime, y) && iszero(x.nanosecs)
Base.isequal(x::Date, y::NanoDate) =
    isequal(x, y.datetime) && iszero(y.nanosecs)

Base.isequal(x::NanoDate, y::DateTime) =
    isequal(x.datetime, y) && iszero(x.nanosecs)
Base.isequal(x::DateTime, y::NanoDate) =
    isequal(x, y.datetime) && iszero(y.nanosecs)

Base.isless(x::NanoDate, y::NanoDate) =
    isless(x.datetime, y.datetime) ||
    (isequal(x.datetime, y.datetime) && isless(x.nanosecs, y.nanosecs))

Base.isless(x::NanoDate, y::Date) = isless(x.datetime, y)
Base.isless(x::Date, y::NanoDate) =
    isless(x, y.datetime) ||
    (isequal(x, y.datetime) && !iszero(y.nanosecs))

Base.isless(x::NanoDate, y::DateTime) = isless(x.datetime, y)
Base.isless(x::DateTime, y::NanoDate) =
    isless(x, y.datetime) ||
    (isequal(x, y.datetime) && !iszero(y.nanosecs))

Base.:(==)(x::NanoDate, y::NanoDate) =
    (x.datetime == y.datetime) && (x.nanosecs == y.nanosecs)
Base.:(<)(x::NanoDate, y::NanoDate) =
    (x.datetime < y.datetime) || (x.datetime == y.datetime && x.nanosecs < y.nanosecs)
Base.:(<=)(x::NanoDate, y::NanoDate) =
    (x.datetime < y.datetime) || (x.datetime == y.datetime && x.nanosecs <= y.nanosecs)
Base.:(>)(x::NanoDate, y::NanoDate) =
    (x.datetime > y.datetime) || (x.datetime == y.datetime && x.nanosecs > y.nanosecs)
Base.:(>=)(x::NanoDate, y::NanoDate) =
    (x.datetime > y.datetime) || (x.datetime == y.datetime && x.nanosecs >= y.nanosecs)

for (T) in (:DateTime, :Date)
  @eval begin
    Base.:(==)(x::NanoDate, y::$T) = (x.datetime == y) && iszero(x.nanosecs)
    Base.:(==)(x::$T, y::NanoDate) = (x == y.datetime) && iszero(y.nanosecs)

    Base.:(<)(x::NanoDate, y::$T) = (x.datetime < y)
    Base.:(<)(x::$T, y::NanoDate) = (x < y.datetime)
    Base.:(>)(x::NanoDate, y::$T) = (x.datetime > y) || (x.datetime == y && !iszero(x.nanosecs))
    Base.:(>)(x::$T, y::NanoDate) = (x > y.datetime)

    Base.:(<=)(x::NanoDate, y::$T) = (x < y) || (x == y)
    Base.:(<=)(x::$T, y::NanoDate) = (x < y) || (x == y)
    Base.:(>=)(x::NanoDate, y::$T) = (x > y) || (x == y)
    Base.:(>=)(x::$T, y::NanoDate) = (x > y) || (x == y)
  end
end


text:
# Command: Get Element Tag Name
"""
    element_tag(element::Element)::String

Get Element Tag Name.
"""
function element_tag(element::Element)::String
    @unpack addr, id = element.session
    element_id = element.id
    response = HTTP.get(
        "$addr/session/$id/element/$element_id/name",
        [("Content-Type" => "application/json; charset=utf-8")],
    )
    @assert response.status == 200
    JSON3.read(response.body).value
end