text:
include("solver_options.jl")
import Base: copy, length, size, reset

"""
$(TYPEDEF)
    Type for solver states
"""
mutable struct SolverState
    constrained::Bool # Constrained solve
    minimum_time::Bool # Minimum time solve
    infeasible::Bool # Infeasible solve

    unconstrained_original_problem::Bool # Original problem is unconstrained but Solve converts to constrained for Infeasible or Minimum Time problem
    second_order_dual_update::Bool # Second order update for dual variables (Lagrange multipliers)
    fixed_constraint_jacobians::Bool # If no custom constraints are provided, all constraint Jacobians are fixed and only need to be updated once
    fixed_terminal_constraint_jacobian::Bool
    penalty_only::Bool  # initial phase where only penalty term is updated each outer loop

    function SolverState()
        new(false,false,false,false,false,false,false,false)
    end
end

function reset(state::SolverState)
    state.constrained = false
    state.minimum_time = false
    state.infeasible = false
    state.unconstrained_original_problem = false
    state.second_order_dual_update = false
    state.fixed_constraint_jacobians = false
    state.fixed_terminal_constraint_jacobian = false
    state.penalty_only = false
    return nothing
end

"""
$(TYPEDEF)
Responsible for storing all solve-dependent variables and solve parameters.
"""
struct Solver{M<:Model,O<:Objective}
    model::M             # Dynamics model
    obj::O               # Objective (cost function and constraints)
    opts::SolverOptions  # Solver options (iterations, method, convergence criteria, etc)
    state::SolverState   # Solver state
    dt::Float64          # Time step
    fd::Function         # Discrete in place dynamics function, `fd(_,x,u)`
    Fd::Function         # Jacobian of discrete dynamics, `fx,fu = F(x,u)`
    Fc::Function         # Jacobian of continuous dynamics, `fx, fu = F(x,u)`
    c_fun::Function
    c_jacobian::Function
    c_labels::Vector{String}  # Constraint labels
    N::Int64             # Number of time steps
    evals::Vector{Int}   # Evaluation counts
    integration::Symbol

    """
    $(SIGNATURES)
    Create a Solver from a model and objective. The user should specify either N (number of knot points) or dt (time step), but typically not both.

    Integration specifies the integration method for discretizing the continuous dynamics and must be (:midpoint, :rk3, or :rk4).

    Solver parameters can be passed in using opts.

    """
    function Solver(model::M, obj::O; integration::Symbol=:rk4, dt::Float64=NaN, N::Int=-1, opts::SolverOptions=SolverOptions()) where {M<:Model,O<:Objective}
        state = SolverState()

        # Check for minimum time
        if obj.tf == 0
            state.minimum_time = true
            dt = 0.
            if N==-1
                throw(ArgumentError("N must be specified for a minimum-time problem"))
            end
        else
            state.minimum_time = false

            # Handle combination of N and dt
            if isnan(dt) && N>0
                dt = obj.tf / (N-1)
            elseif ~isnan(dt) && N==-1
                N, dt = calc_N(obj.tf, dt)
            elseif isnan(dt) && N==-1
                @warn "Neither dt or N were specified. Setting N = 50"
                N = 50
                dt = obj.tf/N
            elseif ~isnan(dt) && N>0
                if dt !== obj.tf/(N-1)
                    throw(ArgumentError("Specified time step, number of knot points, and final time do not agree ($dt ≢ $(obj.tf)/$(N-1))"))
                end
            end
            if dt == 0
                throw(ArgumentError("dt must be non-zero for non-minimum time problems"))
            end
        end

        # Check N, dt for valid entries
        if N < 0
            err = ArgumentError("$N is not a valid entry for N. Number of knot points must be a positive integer.")
            throw(err)
        elseif dt < 0
            err = ArgumentError("$dt is not a valid entry for dt. Time step must be positive.")
            throw(err)
        end

        if O <: ConstrainedObjective
            state.constrained = true
        end

        # Evaluation counts
        evals = zeros(Int,4)  # [Fd,Fc,c,C]
        reset(model)

        # The dynamics function increments an eval counter each time it's called
        f!(xdot,x,u) = dynamics(model,xdot,x,u)

        n, m = model.n, model.m

        m̄ = m
        n̄ = n
        if state.minimum_time
            m̄ += 1
            n̄ += 1
            state.constrained = true
        end

        # Get integration scheme
        if isdefined(TrajectoryOptimization,integration)
            discretizer = eval(integration)
        else
            throw(ArgumentError("$integration is not a defined integration scheme"))
        end

        # Generate discrete dynamics equations
        fd! = discretizer(f!, dt)
        f_aug! = f_augmented!(f!, n, m)

        # Get continuous dynamics jacobian
        Jc = zeros(n+m,n+m)
        Sc = zeros(n+m)
        Scdot = zero(Sc)
        Fc!(Jc,dS,S) = ForwardDiff.jacobian!(Jc,f_aug!,dS,S)
        function fc_jacobians!(x,u)
            # infeasible = size(u,1) != m̄
            Sc[1:n] = x
            Sc[n+1:n+m] = u[1:m]
            Fc!(Jc,Scdot,Sc)
            return Jc[1:n,1:n], Jc[1:n,n+1:n+m] # fx, fu
        end

        """
        s = [x;u;h]
        x ∈ R^n
        u ∈ R^m
        h ∈ R, h = sqrt(dt_k)
        """
        fd_aug! = discretizer(f_aug!)
        nm1 = n + m + 1

        In = 1.0*Matrix(I,n,n)

        # Initialize discrete and continuous dynamics Jacobians
        Jd = zeros(nm1, nm1)
        Sd = zeros(nm1)
        Sdotd = zero(Sd)
        Fd!(Jd,Sdotd,Sd) = ForwardDiff.jacobian!(Jd,fd_aug!,Sdotd,Sd)

        # Discrete dynamics Jacobians
        function fd_jacobians!(fdx,fdu,x,u)
            # Check for infeasible solve
            infeasible = length(u) != m̄

            # Assign state, control (and dt) to augmented vector
            Sd[1:n] = x
            Sd[n+1:n+m] = u[1:m]
            state.minimum_time ? Sd[n+m+1] = u[m̄] : Sd[n+m+1] = √dt

            # Calculate Jacobian
            Fd!(Jd,Sdotd,Sd)

            # if infeasible
            #     return Jd[1:n,1:n], [Jd[1:n,n.+(1:m̄)] I] # fx, [fū I]
            # else
            #     return Jd[1:n,1:n], Jd[1:n,n.+(1:m̄)] # fx, fū
            # end

            fdx[1:n,1:n] = Jd[1:n,1:n]
            fdu[1:n,1:m̄] = Jd[1:n,n.+(1:m̄)]

            if infeasible
                fdu[1:n,m̄+1:m̄+n] = In
            end
            if state.minimum_time
                fdu[n̄,m̄] = 1.
            end

        end
        # Discrete dynamics Jacobian for stacked variables
        function fd_jacobians!(dz,z)
            Sd[1:n+m] = z
            Fd!(Jd,Sdotd,Sd)
            copyto!(dz,Jd[1:n,1:n+m])
        end

        # Generate constraint functions
        c!, c_jacobian!, c_labels = generate_constraint_functions(obj, max_dt = opts.max_dt, min_dt = opts.min_dt)

        # Copy solver options so any changes don't modify the options passed in
        options = copy(opts)

        new{M,O}(model, obj, options, state, dt, fd!, fd_jacobians!, fc_jacobians!, c!, c_jacobian!, c_labels, N, evals, integration)
    end
end

function Solver(solver::Solver; model=solver.model, obj=solver.obj,integration=solver.integration, dt=solver.dt, N=solver.N, opts=solver.opts)
     Solver(model, obj, integration=integration, dt=dt, N=N, opts=opts)
 end

""" $(SIGNATURES) Descriptive labels of the constraints
Any label with a "* " prefix is added by the solver and not in the original
constraints specified by the user
"""
function get_constraint_labels(solver::Solver)
     n,m,N = get_sizes(solver)
     c_labels = copy(solver.c_labels)
     if solver.state.infeasible
         lbl_inf = ["* infeasible control" for i = 1:n]
         append!(c_labels, lbl_inf)
     end
     if solver.state.minimum_time
         push!(c_labels, "* √dt (equality)")
     end
     return c_labels
 end

""" $(SIGNATURES)
Return boolean indices to the "original" constraints,
i.e. the constraints specified by the user (and not added by the solver).
"""
 function original_constraint_inds(solver::Solver)
     labels = get_constraint_labels(solver)
     map(x -> x[1:2] != "* ", labels)
 end

function calc_N(tf::Float64, dt::Float64)::Tuple
    N = convert(Int64,floor(tf/dt)) + 1
    dt = tf/(N-1)
    return N, dt
end

"""$(SIGNATURES) Reset the number of function evaluations and solver state"""
reset(solver::Solver) = begin reset(solver.state); reset_evals(solver) end
reset_evals(solver::Solver) = begin solver.evals .= zeros(4); reset(solver.model) end

"""$(SIGNATURES) Get the number of function evaluations for a given function """
function evals(solver::Solver, fun::Symbol)
    if fun in [:Fd, :discrete_dynamics_jacobian]
        return solver.evals[1]
    elseif fun in [:Fc, :continuous_dynamics_jacobian, :dynamics_jacobian]
        return solver.evals[2]
    elseif fun in [:c, :C, :constraints, :constraint_function, :c_fun]
        return solver.evals[3]
    elseif fun in [:∇c, :constraint_jacobian, :c_jacob, :c_jacobian]
        return solver.evals[4]
    elseif fun in [:f, :dynamics, :continuous_dynamics]
        return evals(solver.model)
    end
end

# Get functions from Solver (incrementing evals)
""" $(SIGNATURES) Return the constraint function, and optionally count the evaluations """
function constraint_function(solver::Solver, count::Bool=true)
    function c_fun_count!(c,x,u)
        solver.c_fun(c,x,u)
        count ? solver.evals[3] += 1 : nothing
    end
    function c_fun_count!(c,x)
        solver.c_fun(c,x)
        count ? solver.evals[3] += 1 : nothing
    end
    return c_fun_count!
end

""" $(SIGNATURES) Return the discrete dynamics jacobian function, and optionally count the evaluations """
function discrete_dynamics_jacobian(solver::Solver, count::Bool=true)
    function Fd!(A,B,x,u)
        solver.Fd(A,B,x,u)
        count ? solver.evals[1] += 1 : nothing
    end
    return Fd!
end

""" $(SIGNATURES) Return the constraint jacobian function, and optionally count the evaluations """
function constraint_jacobian(solver::Solver, count::Bool=true)
    function c_jacob!(Cx,Cu,x,u)
        solver.c_jacobian(Cx,Cu,x,u)
        count ? solver.evals[4] += 1 : nothing
    end
    function c_jacob!(Cx,x)
        solver.c_jacobian(Cx,x)
        count ? solver.evals[4] += 1 : nothing
    end
    return c_jacob!
end

""" $(SIGNATURES) Return the continuous dynamics jacobian function, and optionally count the evaluations """
function dynamics_jacobian(solver::Solver, count::Bool=true)
    function Fc!(x,u)
        count ? solver.evals[2] += 1 : nothing
        solver.Fc(x,u)
    end
    return Fc!
end


text:
# Most of the code here is copied from the Juno codebase

using REPL
using REPL.LineEdit

isREPL() = isdefined(Base, :active_repl) &&
           isdefined(Base.active_repl, :interface) &&
           isdefined(Base.active_repl.interface, :modes)

juliaprompt = "julia> "

current_prompt = juliaprompt

function hideprompt(f)
    isREPL() || return f()

    repl = Base.active_repl
    mistate = repl.mistate
    mode = mistate.current_mode

    buf = String(take!(copy(LineEdit.buffer(mistate))))

    # clear input buffer
    truncate(LineEdit.buffer(mistate), 0)
    LineEdit.refresh_multi_line(mistate)

    print(stdout, "\e[1K\r")
    r = f()

    flush(stdout)
    flush(stderr)
    sleep(0.05)

    # TODO Fix this
    # pos = @rpc cursorpos()
    pos = 1,1
    pos[1] != 0 && println()

    # restore prompt
    if applicable(LineEdit.write_prompt, stdout, mode)
      LineEdit.write_prompt(stdout, mode)
    elseif mode isa LineEdit.PrefixHistoryPrompt || :parent_prompt in fieldnames(typeof(mode))
      LineEdit.write_prompt(stdout, mode.parent_prompt)
    else
      printstyled(stdout, current_prompt, color=:green)
    end

    truncate(LineEdit.buffer(mistate), 0)

    # restore input buffer
    LineEdit.edit_insert(LineEdit.buffer(mistate), buf)
    LineEdit.refresh_multi_line(mistate)
    r
end

function withpath(f, path)
    tls = task_local_storage()
    hassource = haskey(tls, :SOURCE_PATH)
    hassource && (path′ = tls[:SOURCE_PATH])
    tls[:SOURCE_PATH] = path
    try
        return f()
    finally
        hassource ? (tls[:SOURCE_PATH] = path′) : delete!(tls, :SOURCE_PATH)
    end
end

ext: jl

max_stars_repo_name: caleb-allen/julia-vscode

max_issues_repo_name: caleb-allen/julia-vscode

max_forks_repo_name: caleb-allen/julia-vscode


text:
include( "transitionflux.jl" )
include( "nodeflux.jl" )
include( "nodepop.jl" )
include( "nodecomposition.jl" )
include( "subpop.jl" )
include( "subpopage.jl" )
include( "initpop.jl" )

const extensions = ["svg", "html", "png"]
const plotTypes = Dict{Symbol,NTuple{2,Function}}(
    :normal => (normalFluxPlot, normalCompPlot),
    :stacked => (stackedFluxPlot, stackedCompPlot) )
const labelInfix = Dict{Symbol,NTuple{2,String}}( :in => (" into ", " from "),
    :out => (" out of ", " to "), :within => (" within ", " to ") )


text:
#=
Created on Friday 4 December 2020
Last update: Saturday 26 December 2020

@author: Michiel Stock
michielfmstock@gmail.com

@author: Bram De Jaegher
bram.de.jaegher@gmail.com

Automatic grading of answers in our Pluto notebook.
Templates heavily based on the MIT course "Computational Thinking"

https://computationalthinking.mit.edu/Fall20/installation/
=#

mutable struct ProgressTracker
	correct::Int
	total::Int
	name::String
	email::String
	ProgressTracker(name, email) = return new(0, 0, name, email)
end

addQuestion!(t::ProgressTracker) = t.total += 1
accept!(t::ProgressTracker) =	t.correct += 1
Base.show(io::IO, t::ProgressTracker) = print(io, "Notebook of **$(t.name)** with a completion of **$(t.correct) out of $(t.total)** question(s).")

function check_answer(q::Question,  t::ProgressTracker)
	validators = q.validators
	if length(validators) < 1
		return md""
	end

	addQuestion!(t)
	all_valid = all(validators)
	some_valid = any(validators)
	if ismissing(all_valid)
		status = still_missing()
	elseif some_valid && !all_valid
		status = keep_working()
	elseif !all_valid
		status = keep_working()
	elseif all_valid
		accept!(t)
		status = correct()
	end
	return status
end

function check_answer(q::QuestionOptional,  t::ProgressTracker)
	return check_answer(q)
end

function validate(q::QuestionBlock, t::ProgressTracker)
	for (index, question) in enumerate(q.questions)
		q.questions[index].status = check_answer(question, t::ProgressTracker)
	end
	return q
end

function grade(fn)
	open(fn, "a") do f
    write(f, "\ntracker\n")
	end
	include(fn)
	return tracker
end


text:
@doc raw"""
# Planar Point Vortices

"""
module PointVortices

    using GeometricBase
    using GeometricEquations

    export point_vortices_ode, point_vortices_iode, point_vortices_idae,
           point_vortices_dg, point_vortices_formal_lagrangian,
           hamiltonian, angular_momentum, ϑ1, ϑ2, ϑ3, ϑ4

    export compute_energy_error, compute_angular_momentum_error

    const γ₁ = +0.5
    const γ₂ = +0.5
    const X0 = +0.5
    const Y0 = +0.1
    const X1 = +0.5
    const Y1 = -0.1

    function S(x::T, y::T) where {T}
        1 + x^2 + y^2
    end

    function S2(x::T, y::T) where {T}
        1 + 2x^2 + 2y^2
    end

    function dSdx(x::T, y::T) where {T}
        2x
    end

    function dSdy(x::T, y::T) where {T}
        2y
    end

    function hamiltonian(t,q)
        γ₁ * γ₂ * S(q[1],q[2]) * S(q[3],q[4]) * log( (q[1] - q[3])^2 + (q[2] - q[4])^2 ) / (2π)
    end

    function angular_momentum(t,q)
        q[1] * ϑ2(t, q) - q[2] * ϑ1(t, q) +
        q[3] * ϑ4(t, q) - q[4] * ϑ3(t, q)
    end

    function ϑ1(t, q)
        - γ₁ * q[2] * S(q[1], q[2]) / 2
    end

    function ϑ2(t, q)
        + γ₁ * q[1] * S(q[1], q[2]) / 2
    end

    function ϑ3(t, q)
        - γ₂ * q[4] * S(q[3], q[4]) / 2
    end

    function ϑ4(t, q)
        + γ₂ * q[3] * S(q[3], q[4]) / 2
    end

    const q₀ = [X0, Y0, X1, Y1]
    const p₀ = [ϑ1(0, q₀), ϑ2(0, q₀), ϑ3(0, q₀), ϑ4(0, q₀)]

    function dϑ1d1(t, q)
        - γ₁ * q[2] * dSdx(q[1], q[2]) / 2
    end

    function dϑ1d2(t, q)
        - γ₁ * q[2] * dSdy(q[1], q[2]) / 2 - γ₁ * S(q[1], q[2]) / 2
    end

    function dϑ1d3(t, q)
        zero(eltype(q))
    end

    function dϑ1d4(t, q)
        zero(eltype(q))
    end

    function dϑ2d1(t, q)
        + γ₁ * q[1] * dSdx(q[1], q[2]) / 2 + γ₁ * S(q[1], q[2]) / 2
    end

    function dϑ2d2(t, q)
        + γ₁ * q[1] * dSdy(q[1], q[2]) / 2
    end

    function dϑ2d3(t, q)
        zero(eltype(q))
    end

    function dϑ2d4(t, q)
        zero(eltype(q))
    end

    function dϑ3d1(t, q)
        zero(eltype(q))
    end

    function dϑ3d2(t, q)
        zero(eltype(q))
    end

    function dϑ3d3(t, q)
        - γ₂ * q[4] * dSdx(q[3], q[4]) / 2
    end

    function dϑ3d4(t, q)
        - γ₂ * q[4] * dSdy(q[3], q[4]) / 2 - γ₂ * S(q[3], q[4]) / 2
    end

    function dϑ4d1(t, q)
        zero(eltype(q))
    end

    function dϑ4d2(t, q)
        zero(eltype(q))
    end

    function dϑ4d3(t, q)
        + γ₂ * q[3] * dSdx(q[3], q[4]) / 2 + γ₂ * S(q[3], q[4]) / 2
    end

    function dϑ4d4(t, q)
        + γ₂ * q[3] * dSdy(q[3], q[4]) / 2
    end

    function ϑ(t, q, p)
        p[1] = ϑ1(t, q)
        p[2] = ϑ2(t, q)
        p[3] = ϑ3(t, q)
        p[4] = ϑ4(t, q)
    end

    function ω(t, q, Ω)
        Ω[1,1] = 0
        Ω[1,2] = dϑ1d2(t,q) - dϑ2d1(t,q)
        Ω[1,3] = dϑ1d3(t,q) - dϑ3d1(t,q)
        Ω[1,4] = dϑ1d4(t,q) - dϑ4d1(t,q)

        Ω[2,1] = dϑ2d1(t,q) - dϑ1d2(t,q)
        Ω[2,2] = 0
        Ω[2,3] = dϑ2d3(t,q) - dϑ3d2(t,q)
        Ω[2,4] = dϑ2d4(t,q) - dϑ4d2(t,q)

        Ω[3,1] = dϑ3d1(t,q) - dϑ1d3(t,q)
        Ω[3,2] = dϑ3d2(t,q) - dϑ2d3(t,q)
        Ω[3,3] = 0
        Ω[3,4] = dϑ3d4(t,q) - dϑ4d3(t,q)

        Ω[4,1] = dϑ4d1(t,q) - dϑ1d4(t,q)
        Ω[4,2] = dϑ4d2(t,q) - dϑ2d4(t,q)
        Ω[4,3] = dϑ4d3(t,q) - dϑ3d4(t,q)
        Ω[4,4] = 0

        nothing
    end

    function f1(t, q, v)
        γ₁ * ( dSdx(q[1],q[2]) * (q[1] * v[2] - q[2] * v[1]) + v[2] * S(q[1], q[2]) ) / 2
    end

    function f2(t, q, v)
        γ₁ * ( dSdy(q[1],q[2]) * (q[1] * v[2] - q[2] * v[1]) - v[1] * S(q[1], q[2]) ) / 2
    end

    function f3(t, q, v)
        γ₂ * ( dSdx(q[3],q[4]) * (q[3] * v[4] - q[4] * v[3]) + v[4] * S(q[3], q[4]) ) / 2
    end

    function f4(t, q, v)
        γ₂ * ( dSdy(q[3],q[4]) * (q[3] * v[4] - q[4] * v[3]) - v[3] * S(q[3], q[4]) ) / 2
    end

    function dHd1(t, q)
        + γ₁ * γ₂ * dSdx(q[1],q[2]) * S(q[3],q[4]) * log( (q[1] - q[3])^2 + (q[2] - q[4])^2 ) / (2π) +
          γ₁ * γ₂ * S(q[1],q[2]) * S(q[3],q[4]) * (q[1] - q[3]) / ( (q[1] - q[3])^2 + (q[2] - q[4])^2 ) / π
    end

    function dHd2(t, q)
        + γ₁ * γ₂ * dSdy(q[1],q[2]) * S(q[3],q[4]) * log( (q[1] - q[3])^2 + (q[2] - q[4])^2 ) / (2π) +
          γ₁ * γ₂ * S(q[1],q[2]) * S(q[3],q[4]) * (q[2] - q[4]) / ( (q[1] - q[3])^2 + (q[2] - q[4])^2 ) / π
    end

    function dHd3(t, q)
        + γ₁ * γ₂ * dSdx(q[3],q[4]) * S(q[1],q[2]) * log( (q[1] - q[3])^2 + (q[2] - q[4])^2 ) / (2π) -
          γ₁ * γ₂ * S(q[1],q[2]) * S(q[3],q[4]) * (q[1] - q[3]) / ( (q[1] - q[3])^2 + (q[2] - q[4])^2 ) / π
    end

    function dHd4(t, q)
        + γ₁ * γ₂ * dSdy(q[3],q[4]) * S(q[1],q[2]) * log( (q[1] - q[3])^2 + (q[2] - q[4])^2 ) / (2π) -
          γ₁ * γ₂ * S(q[1],q[2]) * S(q[3],q[4]) * (q[2] - q[4]) / ( (q[1] - q[3])^2 + (q[2] - q[4])^2 ) / π
    end

    function dH(t, q, dH)
        dH[1] = dHd1(t, q)
        dH[2] = dHd2(t, q)
        dH[3] = dHd3(t, q)
        dH[4] = dHd4(t, q)
        nothing
    end

    function point_vortices_p₀(q₀, t₀=0)
        p₀ = zero(q₀)
        tq = zeros(eltype(q₀), size(q₀,1))
        tp = zeros(eltype(p₀), size(p₀,1))

        if ndims(q₀) == 1
            ϑ(t₀, q₀, p₀)
        else
            for i in 1:size(q₀,2)
                simd_copy_xy_first!(tq, q₀, i)
                ϑ(t₀, tq, tp)
                simd_copy_yx_first!(tp, p₀, i)
            end
        end
        p₀
    end

    function point_vortices_v(t, q, v)
        denominator1 = 1 / (γ₁ * S2(q[1], q[2]))
        denominator2 = 1 / (γ₂ * S2(q[3], q[4]))

        v[1] = - dHd2(t,q) * denominator1
        v[2] = + dHd1(t,q) * denominator1
        v[3] = - dHd4(t,q) * denominator2
        v[4] = + dHd3(t,q) * denominator2

        nothing
    end

    function point_vortices_v(t, q, p, v)
        point_vortices_v(t, q, v)
    end

    function point_vortices_ϑ(t, q, p)
        p[1] = ϑ1(t,q)
        p[2] = ϑ2(t,q)
        p[3] = ϑ3(t,q)
        p[4] = ϑ4(t,q)
        nothing
    end

    function point_vortices_ϑ(t, q, v, p)
        point_vortices_ϑ(t, q, p)
    end

    function point_vortices_f(t, q, v, f)
        f[1] = f1(t,q,v) - dHd1(t,q)
        f[2] = f2(t,q,v) - dHd2(t,q)
        f[3] = f3(t,q,v) - dHd3(t,q)
        f[4] = f4(t,q,v) - dHd4(t,q)
        nothing
    end

    function point_vortices_f(t, q, v, p, f)
        point_vortices_f(t, q, v, f)
    end

    function point_vortices_g(t, q, λ, g)
        g[1] = f1(t,q,λ)
        g[2] = f2(t,q,λ)
        g[3] = f3(t,q,λ)
        g[4] = f4(t,q,λ)
        nothing
    end

    function point_vortices_g(t, q, p, λ, g)
        point_vortices_g(t, q, λ, g)
    end

    function point_vortices_u(t, q, v, u)
        u[1] = v[1]
        u[2] = v[2]
        u[3] = v[3]
        u[4] = v[4]
        nothing
    end

    function point_vortices_u(t, q, p, v, u)
        point_vortices_u(t, q, v, u)
    end

    function point_vortices_ϕ(t, q, p, ϕ)
        ϕ[1] = p[1] - ϑ1(t,q)
        ϕ[2] = p[2] - ϑ2(t,q)
        ϕ[3] = p[3] - ϑ3(t,q)
        ϕ[4] = p[4] - ϑ4(t,q)
        nothing
    end

    function point_vortices_ode(q₀=q₀)
        ODE(point_vortices_v, q₀)
    end

    function point_vortices_iode(q₀=q₀, p₀=point_vortices_p₀(q₀))
        IODE(point_vortices_ϑ, point_vortices_f,
             point_vortices_g, q₀, p₀;
             v̄=point_vortices_v)
    end

    function point_vortices_idae(q₀=q₀, p₀=point_vortices_p₀(q₀), λ₀=zero(q₀))
        IDAE(point_vortices_ϑ, point_vortices_f,
             point_vortices_u, point_vortices_g,
             point_vortices_ϕ, q₀, p₀, λ₀;
             v̄=point_vortices_v)
    end

    function point_vortices_dg(q₀=q₀)
        IODE(point_vortices_ϑ, point_vortices_f,
             point_vortices_g, q₀, q₀;
             v=point_vortices_v)
    end

    function point_vortices_formal_lagrangian(q₀=q₀)
        p₀ = point_vortices_p₀(q₀)
        LODE(ϑ, point_vortices_f, point_vortices_g, q₀, p₀;
             v̄=point_vortices_v, Ω=ω, ∇H=dH)
    end

    function compute_energy_error(t, q::AbstractDataSeries{T}) where {T}
        h = DataSeries(T, q.nt)
        e = DataSeries(T, q.nt)

        for i in axes(q,2)
            h[i] = hamiltonian(t[i], q[:,i])
            e[i] = (h[i] - h[0]) / h[0]
        end

        (h, e)
    end

    function compute_angular_momentum_error(t, q::AbstractDataSeries{T}) where {T}
        m = DataSeries(T, q.nt)
        e = DataSeries(T, q.nt)

        for i in axes(q,2)
            m[i] = angular_momentum(t[i], q[:,i])
            e[i] = (m[i] - m[0]) / m[0]
        end

        (m, e)
    end

end


text:
"""
    _fst_permutation(population_1::T, population_2::T) where T<:AbstractMatrix
Returns two matrices with rows (samples) shuffled between them. Respects the
number of rows of the original matrices (i.e. population sizes).
"""
function _fst_permutation(data::T, n1::Integer, n2::Integer) where T<:AbstractMatrix
    # total number of samples
    n_range = 1:(n1 +n2)
    # generate permutation indices for the first "population"
    perm_matrix = @inbounds data[shuffle(n_range), :]
    #partition shuffled matrix into original sample sizes
    new_pop_1, new_pop_2 = partitionarray(perm_matrix, [n1, n2])
    return new_pop_1, new_pop_2
end

function _fst_permutation(data::T) where T<:AbstractMatrix
    @inbounds data[shuffle(1:size(data,1)), :]
end

function _permuted_Hudson(data::PopData, iterations::Int64)
    !isbiallelic(data) && throw(error("Data must be biallelic to use the Hudson estimator"))
    idx_pdata = groupby(data.genodata, :population)
    pops = getindex.(keys(idx_pdata), :population)
    npops = data.metadata.populations
    n_loci = data.metadata.loci
    results = zeros(npops, npops)
    perm_vector = zeros(iterations-1)
    p = Progress(sum(1:(npops - 1)) * (iterations - 1), dt = 1, color = :blue, barglyphs = BarGlyphs("|=> |"), barlen = 30)
    @inbounds for i in 2:npops
        pop1 = reshape(idx_pdata[i].genotype, :, n_loci)
        n_pop1 = size(pop1, 1)
        @inbounds for j in 1:(i-1)
            pop2 = reshape(idx_pdata[j].genotype, :, n_loci)
            n_pop2 = size(pop2, 1)
            merged = vcat(pop1, pop2)
            fst_val = _hudson_fst(pop1,pop2)
            @inbounds @sync for iter in 1:iterations-1
                Base.Threads.@spawn begin
                    perm_p1, perm_p2 = _fst_permutation(merged, n_pop1, n_pop2)
                    @inbounds perm_vector[iter] = _hudson_fst(perm_p1, perm_p2)
                    pops_text = string(pops[i]) * " & " * string(pops[j])
                    ProgressMeter.next!(p; showvalues = [(:Populations, pops_text), (:Iteration, "$iter")])
                end
            end
            @inbounds results[i,j] = fst_val
            @inbounds results[j,i] = (sum(fst_val .<= perm_vector) + 1) / iterations
        end
    end
    println("Below diagonal: FST values | Above diagonal: P values")
    return PairwiseFST(DataFrame(results, Symbol.(pops)), "Hudson et al. 1992 (with p-values)")
end

function _permuted_Nei(data::PopData, iterations::Int64)
    idx_pdata = groupby(data.genodata, :population)
    pops = getindex.(keys(idx_pdata), :population)
    npops = data.metadata.populations
    n_loci = data.metadata.loci
    results = zeros(npops, npops)
    perm_vector = zeros(iterations-1)
    p = Progress(sum(1:(npops - 1)) * (iterations - 1), dt = 1, color = :blue, barglyphs = BarGlyphs("|=> |"), barlen = 30)
    @inbounds for i in 2:npops
        pop1 = reshape(idx_pdata[i].genotype, :, n_loci)
        n_pop1 = size(pop1, 1)
        @inbounds for j in 1:(i-1)
            pop2 = reshape(idx_pdata[j].genotype, :, n_loci)
            n_pop2 = size(pop2, 1)
            merged = vcat(pop1, pop2)
            fst_val = _nei_fst(pop1,pop2)
            @inbounds @sync for iter in 1:iterations-1
                Base.Threads.@spawn begin
                    perm_p1, perm_p2 = _fst_permutation(merged, n_pop1, n_pop2)
                    @inbounds perm_vector[iter] = _nei_fst(perm_p1, perm_p2)
                    pops_text = string(pops[i]) * " & " * string(pops[j])
                    ProgressMeter.next!(p; showvalues = [(:Populations, pops_text), (:Iteration, "$iter")])
                end
            end
            @inbounds results[i,j] = fst_val
            @inbounds results[j,i] = (sum(fst_val .<= perm_vector) + 1) / iterations
        end
    end
    println("Below diagonal: FST values | Above diagonal: P values")
    return PairwiseFST(DataFrame(results, Symbol.(pops)), "Nei 1987 (with p-values)")
end

function _permuted_WeirCockerham(data::PopData, iterations::Int64)
    idx_pdata = groupby(data.genodata, :population)
    pops = getindex.(keys(idx_pdata), :population)
    npops = data.metadata.populations
    n_loci = data.metadata.loci
    results = zeros(npops, npops)
    perm_vector = zeros(iterations-1)
    p = Progress(sum(1:(npops - 1)) * (iterations - 1), dt = 1, color = :blue, barglyphs = BarGlyphs("|=> |"), barlen = 30)
    @inbounds for i in 2:npops
        pop1 = reshape(idx_pdata[i].genotype, :, n_loci)
        n_pop1 = size(pop1, 1)
        @inbounds for j in 1:(i-1)
            pop2 = reshape(idx_pdata[j].genotype, :, n_loci)
            n_pop2 = size(pop2, 1)
            merged = vcat(pop1, pop2)
            fst_val = _weircockerham_fst(pop1,pop2)
            @inbounds @sync for iter in 1:iterations-1
                Base.Threads.@spawn begin
                    perm_p1, perm_p2 = _fst_permutation(merged, n_pop1, n_pop2)
                    @inbounds perm_vector[iter] = _weircockerham_fst(perm_p1, perm_p2)
                    pops_text = string(pops[i]) * " & " * string(pops[j])
                    ProgressMeter.next!(p; showvalues = [(:Populations, pops_text), (:Iteration, "$iter")])
                end
            end
            @inbounds results[i,j] = fst_val
            @inbounds results[j,i] = (sum(fst_val .<= perm_vector) + 1) / iterations
        end
    end
    println("Below diagonal: FST values | Above diagonal: P values")
    return PairwiseFST(DataFrame(results, Symbol.(pops)), "Weir & Cockerham (with p-values)")
end


text:
# NOTE: Tests are broken where commented out.

# For testing: include(Pkg.dir("CodeGen", "test/runtests.jl"))

using Test
using CodeGen
using LLVM

Base.CoreLogging.global_logger(Base.CoreLogging.SimpleLogger(STDERR))
# Base.CoreLogging.global_logger(Base.CoreLogging.SimpleLogger(STDERR, Base.CoreLogging.Debug))

"""
    @cgtest fun(args...)

Test if `fun(args...)` is equal to `CodeGen.run(fun, args...)`
"""
macro cgtest(e)
    f = e.args[1]
    args = length(e.args) > 1 ? e.args[2:end] : Any[]
    quote
        @test $(esc(e)) == @jitrun($(esc(f)), $(esc.(args)...))
    end
end

f(x) = 3
m = codegen(f, Tuple{Int})
@cgtest f(2)

array_max(x) = maximum(Int[4,3x])
m = codegen(array_max, Tuple{Int})
verify(m)
@test @jitrun(array_max, -1) == array_max(-1)
@cgtest array_max(2)
@cgtest array_max(-1)

function variables(i, j, k)
    i = 2i + k
    l = j + i
    return l
end
@cgtest variables(3,4,5)
@test variables(3,4,5) == @jitrun(variables,3,4,5)

function variables2(i)
    x = 5
    y = 9i
    x = x + y
    z = 2x
    return z
end
@cgtest variables2(3)

f() = 99.0
mod = codegen(f, Tuple{})
optimize!(mod)
# @cgtest f()                # both @cgtest and @jitrun are broken with no args
# @test f() == @jitrun(f)

fx(x) = 2x + 50
mod = codegen(fx, Tuple{Int})
optimize!(mod)
@test @jitrun(fx, 10) == fx(10)
# The following is the same test:
@cgtest fx(10)
@cgtest fx(10.0)

z = codegen(abs, Tuple{Float64})
optimize!(z)
verify(z)
@cgtest abs(-10.0)
@cgtest abs(10.0)
@cgtest abs(-10)
@cgtest abs(10)

sum_tuple(x) = sum((x, x, 1.0))
codegen(sum_tuple, Tuple{Float64})
codegen(sum_tuple, Tuple{Float32})
codegen(sum_tuple, Tuple{Int64})
codegen(sum_tuple, Tuple{Int32})
# codegen(sum_tuple, Tuple{Complex128})    # can't box
@cgtest sum_tuple(5)
@cgtest sum_tuple(5.5)

array_max(x) = maximum(Int[3,x])
@cgtest array_max(1)
@cgtest array_max(4)

array_sum(x) = sum(Int[3,x])
# m = codegen(array_sum, Tuple{Int})
# verify(m)
# @cgtest array_sum(1)

function an_if(x)
    return x == 3 ? 4 : 5
end
@cgtest an_if(3)
@cgtest an_if(5)

function an_if2(x)
    if x < 3
        x += 3
    end
    return x+2
end
@cgtest an_if2(5)
@cgtest an_if2(1)

function ifs(x)
    if x < 7
        if x == 5
            x += 2
        end
        x += 2
    end
    return x+2
end
@cgtest ifs(0)
@cgtest ifs(5)
@cgtest ifs(7)
@cgtest ifs(8)
@cgtest ifs(0.0)
@cgtest ifs(5.0)
@cgtest ifs(7.0)
@cgtest ifs(8.0)

function ifs2(x)
    if x > 3
        z = 2x
    else
        z = x-4
    end
    return z+2
end
@cgtest ifs2(5)
@cgtest ifs2(3)
@cgtest ifs2(1)

function while_loop(i)
    x = 5
    while i < 4
        i += 1
        x += 2
    end
    return x
end
@cgtest while_loop(0)

function for_loop(x)
    a = 3
    for i in 1:5
        x += i * a
    end
    x
end
@cgtest for_loop(2)
@cgtest for_loop(2.2)

@noinline f(x) = 2x
function call_another(x)
    y = f(x)
    return x + y
end
@cgtest call_another(3)

function check_identity(x)
    return x === 3
end
codegen(check_identity, Tuple{Int})

ccall_cos(x) = ccall(:cos, Float64, (Float64,), x)
@cgtest ccall_cos(1.1)

abstract type AbstractMA end
mutable struct MA <: AbstractMA
    xx::Int
    yy::Float64
end
@noinline f(x::MA) = x.yy
function newmdt(x)
    a = MA(1, x)
    f(a)
end
m = codegen(newmdt, Tuple{Float64})
@cgtest newmdt(1.1)

abstract type AbstractA end
struct A <: AbstractA
    xx::Int
    yy::Float64
end
@noinline f(x::A) = x.yy + x.xx
function newdt(x)
    a = A(3, 2x)
    f(a)
end
@cgtest newdt(1.1)

@noinline tf(x) = x[2]
function tuple_fun(x)
    t = (x, 2x, 1.0)
    return tf(t)
end
@cgtest tuple_fun(1)

@noinline tf2(x) = x[1][1] + x[2]
function tuple_fun2(x)
    t = ((x,),3)
    return tf2(t)
end
@cgtest tuple_fun(2)

function test_arrays(x)
    y = fill(2pi, 5)
    push!(y, 3x)
    z = reverse(y)
    zz = y .+ z.^2
    return maximum(zz)
end
@cgtest test_arrays(2.2)

function array_len(x)
    a = Int[3,4,5]
    return length(a)
end
@cgtest array_len(2)

rfib(n) = n < 2 ? n : rfib(n-1) + rfib(n-2)
@cgtest rfib(9)

@noinline mdisp(x) = 4x
@noinline mdisp(x,y) = x + y
test_dispatch(x,y) = mdisp(x, y) + mdisp(x) + mdisp(y)
@cgtest test_dispatch(1,2)
@cgtest test_dispatch(1,2.0)

@cgtest Base.shl_int(4096,2)
@cgtest Base.shl_int(UInt(4096),2)
@cgtest Base.shl_int(UInt16(4096),2)
@cgtest Base.shl_int(4096,UInt8(2))
@cgtest Base.lshr_int(4096,2)
@cgtest Base.lshr_int(UInt(4096),2)
@cgtest Base.lshr_int(UInt16(4096),2)
@cgtest Base.lshr_int(4096,UInt8(2))
@cgtest Base.ashr_int(4096,2)
@cgtest Base.ashr_int(UInt(4096),2)
@cgtest Base.ashr_int(UInt16(4096),2)
@cgtest Base.ashr_int(4096,UInt8(2))

mutable struct X
    x::Int
end
@noinline f(x) = X(x)
g(x) = f(x).x + x
m = codegen(g, Tuple{Int})
# verify(m)
@cgtest g(2)

function varargs_fun(x, y::Int...)
    a = y[1]
    b = 3
    return x + a + b
end
@cgtest varargs_fun(2, 3)
@cgtest varargs_fun(2, 3, 4)

nothing


text:
@testset "Overdetermined systems" begin
    @polyvar x y z
    f = [x * z - y^2, y - z^2, x - y * z, x + y + z + 1]

    tracker, starts = pathtracker_startsolutions(f; seed = 213412)
    @test count(is_success, track.(tracker, starts)) == 3
    @test count(r -> r.return_code == :excess_solution, track.(tracker, starts)) == 2

        # check that it also works with polyhedral homotopy
    tracker, starts = pathtracker_startsolutions(
        f;
        seed = 213412,
        start_system = :polyhedral,
    )
    HC.prepare!(tracker, starts)
    @test count(is_success, track.(tracker, starts)) == 3

    minors = include(joinpath(@__DIR__, "examples", "3_by_5_minors.jl"))
    tracker, starts = pathtracker_startsolutions(minors)
    @test count(is_success, track.(tracker, starts)) == 80

        # chemical reaction network examples
    g = let
        @polyvar x[1:3]
        p = [0.04, 0.04, 1.0, 1.0, 10.0, 0.0, 0.04, 35.0, 0.1, 0.04]
        [
         -x[1] * x[3] * p[3] - x[1] * p[2] + p[1],
         x[1] * x[2] * x[3] * p[8] * p[9] + x[1] * x[3] * p[7] * p[8] * p[9] -
         x[2] * x[3] * p[5] * p[6] - x[2] * p[5] * p[6] * p[10] - x[2] * x[3] * p[4] -
         x[2] * p[4] * p[10],
         -x[1] * x[2] * x[3] * p[8] * p[9] - x[1] * x[3] * p[7] * p[8] * p[9] +
         x[2] * x[3] * p[5] * p[6] + x[2] * p[5] * p[6] * p[10] + x[2] * x[3] * p[4] +
         x[2] * p[4] * p[10],
         x[2] + x[3] - 1.0,
        ]
    end

    tracker, starts = pathtracker_startsolutions(g; seed = 213412)
    @test count(is_success, track.(tracker, starts)) == 3
end


text:
using SimpleNeRF: NeRFDataset, CameraView, FileNeRF, ModelMetadata, iterate_batches, RGB, channelview
using Random
using Statistics
using LinearAlgebra
using Test

cv1 = CameraView(camera_direction=(0.0, 1.0, 0.0),
    camera_origin=(2.0, 2.0, 2.0),
    x_axis=(-1.0, 0.0, 0.0),
    y_axis=(0.0, 0.0, 1.0),
    x_fov=60.0,
    y_fov=60.0,)

cv2 = CameraView(camera_direction=(1.0, 0.0, 0.0),
    camera_origin=(-2.0, 2.0, 2.0),
    x_axis=(-0.0, 0.0, -1.0),
    y_axis=(0.0, 1.0, 0.0),
    x_fov=60.0,
    y_fov=60.0,)

file1 = FileNeRF(joinpath(@__DIR__, "testimg1.png"), cv1)
file2 = FileNeRF(joinpath(@__DIR__, "testimg2.png"), cv2)
ds = NeRFDataset(ModelMetadata((0.0, 0.0, 0.0), (1.0, 1.0, 1.0)), [file1, file2])
batches = mktempdir() do tmp_dir
    [d for d in iterate_batches(ds, tmp_dir; rng=MersenneTwister(42), batch_size=51, repeat=false)]
end

@test length(batches) == 4
@test size(batches[end])[end] == 200 - 51*3

combined = cat(batches...; dims=3)

for view in ds.views
    # samples w.r.t origin
    origin = collect(view.cameraview.camera_origin)
    origins = combined[:, 1, :]
    view_mask = sum(abs, origins .- origin; dims=1) .< 1e-5
    count = sum(view_mask)
    num_pixels = prod(size(SimpleNeRF.image(view)))
    @test round(Int, count) == num_pixels

    view_rays = combined[:, :, vec(view_mask)]

    # camera direction
    directions = view_rays[:, 2, :]
    mean_direction = mean(directions; dims=ndims(directions))
    mean_direction /= norm(mean_direction)
    camera_dot = sum(mean_direction .* collect(view.cameraview.camera_direction))
    @test abs(Float32(camera_dot) - 1) < 1e-5

    # color validity
    colors = view_rays[:, 3, :]
    mean_color = mean(colors; dims=ndims(colors))
    actual_mean = mean(channelview(RGB.(SimpleNeRF.image(view))); dims=(2, 3))
    diff = mean(abs, mean_color .- actual_mean)
    @test Float32(diff) < 1e-5
end


text:
# Default is to fallback to the inner call
#
# Do this fallback scheme to specialize for certain combinations of arguments if needed.
factory(args...; kw...) = _factory(args...; kw...)

# Ratio optimizers go through a refinement step
function factory(
        backend::nGraph.Backend{nGraph.CPU},
        func,
        opt::AbstractOptimizer{Rational{Int64}};
        kw...
    )
    return ratiosearch(_factory, backend, func, opt; kw...)
end

struct CompilerRetry <: Exception end
function _factory(
        backend::nGraph.Backend,
        func,
        opt::T;
        # Useful for the GPU case
        adjust_io = false,
        defrag = true,
        just_profile = false,
        profile_kw...
    ) where {T <: AbstractOptimizer}

    # Get the function, arguments, and keyword arguments from the provided function
    Profiler.enable_passes()
    A = func()::Actualizer

    # add a callback that will populate a reference to a `FunctionData` type
    frame_ref = Ref{Frame}()
    limits_ref = Ref{Vector{Int}}()
    creation_times = Float64[]
    optimization_times = Float64[]
    remote_args_ref = Ref{Set{TensorDescriptor}}()

    # A callback that profiles the ngraph function
    function cb(f::nGraph.NFunction)
        # Do some minor editing the order of nodes in the graph to hopefully yield slightly
        # better memory characteristics
        priority_pass!(f)
        data = profile(f, backend; profile_kw...)
        just_profile && throw(CompilerExit())

        # Initialize the node dram limits if needed
        if !isdefined(limits_ref, :x)
            # Get the limit from the optimizer
            # Find the input and output size of the function and subtract that from the
            # limit along with a fudge factor so we can fit the model on the GPU
            if adjust_io
                io_size = sum(sizeof, input_tensors(f)) + sum(sizeof, output_tensors(f))
                @show io_size
                limit = max(getlimit(opt) - io_size, 0)
                opt_adjusted = _optimizer(opt, limit)
                modeltype = opt_adjusted(data, backend)
            else
                modeltype = opt(data, backend)
            end

            # Save the limits to the reference for reuse across defragmentation.
            limits_ref[] = modeltype.dram_limits
        else
            modeltype = opt(data, backend)
            modeltype.dram_limits = limits_ref[]
        end

        # Record statistics about how long solving takes.
        creation_time = @elapsed(frame = create_model(modeltype, data))
        optimization_time = @elapsed(optimize!(frame))

        # update the frame with the local args
        empty!(frame.local_args)
        for tensor in tensors(data)
            islocalarg(frame, tensor) && push!(frame.local_args, tensor)
        end

        push!(creation_times, creation_time)
        push!(optimization_times, optimization_time)

        remote_args = configure!(f, frame)
        frame_ref[] = frame
        remote_args_ref[] = remote_args

        return nothing
    end

    # Defrag callback - if a function needs defragging, throws a `CompilerExit` exception to
    # avoid nGraph trying to allocate too much GPU memory
    function defrag_cb(f::nGraph.NFunction)
        if exceeds_limit(f, frame_ref[].modeltype, frame_ref[].local_args)
            # This is pretty ugly - sorry about that.
            modeltype = update(
                frame_ref[].modeltype,
                frame_ref[].local_args,
                profile(f, backend; profile_kw...)
            )
            limits_ref[] = modeltype.dram_limits

            throw(CompilerRetry())
        end
    end

    # Create a function to let the nGraph.jl compiler know if a function parameter or
    # output is supposed to be remote
    isremote(x::nGraph.Node) = in(first(outputs(x)), remote_args_ref[])

    # Compile the function to a ngraph executable
    local fex
    retry = true
    while retry
        retry = false

        # Setup callbacks
        #
        # If the function needs defragging, a `CompilerExit` exception will be thrown and we
        # will have to try again.
        callbacks = CallbackChain()
        callback!(callbacks, cb)
        defrag && callback!(callbacks, defrag_cb)

        try
            #pre_args = nGraph._compile_snoop(backend, f, args...; isremote = isremote)

            fex = nGraph.compile(
                backend,
                A.f,
                A.args...;
                callback = callbacks,
                emit_timing = true,
                #isremote = isremote,
                A.kw...
            )
        catch e
            isa(e, CompilerExit) && return nothing
            isa(e, CompilerRetry) || rethrow(e)
            retry = true
        end
    end

    # Following compilation, set up the appropriate arguments to live in persistent memory
    # or not.
    remote_args = remote_args_ref[]

    metadata = Dict(
        :creation_times => creation_times,
        :optimization_times => optimization_times,
    )

    return fex, frame_ref[], metadata
end

#####
##### Ratio Searching
#####

"""
    ratiosearch(f, backend, func, opt; search_ratio = true, refinements = 7, kw...)

`f` -
`func` - An AutoTM compatible function constructor

Further keywords get passed to the inner call to `f`.
"""
function ratiosearch(f, backend, func, opt; search_ratio = true, refinements = 7, kw...)
    @info "Trying Ratio $(getratio(opt))"

    # Just return the inner factory if we aren't interesting in performing a binary search
    # for the actual ratio to input that will return the desired ratio
    search_ratio || return f(backend, func, opt; kw...)

    # Perform a binary search
    ret = f(backend, func, opt; kw...)
    fex = first(ret)
    args = Base.tail(ret)

    # If we're within the desired tolerance, just return like normal
    checkmargin(fex, opt) && return (fex, args...)

    # For now, I'm assuming that the ratio generated by the optimized graph will always
    # be greater than the desired ratio due to defragmentation.
    mul = getratio(fex) > getratio(opt) ? 1 : -1

    # Start doing a grid search
    #
    # The function desired_ratio -> actual_ratio is not necessarily monotonic, so we can't
    # use a binary search. Instead this grid search thing seems to work alright.
    best_ratio = getratio(opt)
    best_err = geterr(fex, opt)
    current_ratio = best_ratio

    best_fex = fex
    best_args = args

    for i in 1:refinements
        # Use a step size starting with 1 and increasing or decreasing by the step size
        # until the ratio crosses the boundary of what we want.
        step = 1 // (2 ^ (i))
        @info """
        ------------------------
        Performing Refinement Iteration $i
        Step: $step
        ------------------------
        """

        for _ in 1:1
            current_ratio = current_ratio - (mul * step)
            current_ratio < 0 && break

            @info "Trying Ratio: $(convert(Float64, current_ratio))"
            ret = f(backend, func, _optimizer(opt, current_ratio); kw...)
            fex = first(ret)
            args = Base.tail(ret)

            @show typeof(fex)
            @show typeof(args)

            # If the ratios switch sign, time to exit
            if mul == 1
                getratio(fex) <= getratio(opt) && break
            elseif mul == -1
                getratio(fex) >= getratio(opt) && break
            end

            current_err = geterr(fex, opt)
            @info """
            Current Ratio: $(convert(Float64, current_ratio))
            Best Ratio: $(convert(Float64, best_ratio))

            Current Error: $(convert(Float64, current_err))
            Best Error: $(convert(Float64, best_err))
            """
            if current_err < best_err
                best_ratio = current_ratio
                best_err = current_err
                best_fex = fex
                best_args = args
            end
        end
        current_ratio = best_ratio
    end

    return (best_fex, best_args...)
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
module Disease
    export epidemic

    function epidemic(tm::Number, n::Number, s0::Number, i0::Number, b::Number, a::Number)
        dt::Float64 = tm/n
        dta::Float64 = dt * a
        dtb::Float64 = -dt * b
        s::Float64 = s0
        i::Float64 = i0
        r::Float64 = 0
        for t::Float64 = 0:dt:tm
            ds = dtb * s * i
            dr = dta * i
            di = -ds -dr
            if di < 0
                break
            end
            s += ds
            r += dr
            i += di
        end
        return i
    end
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
# This file is a part of Actors.jl, licensed under the MIT License (MIT).

__precompile__(true)

module Actors

include("actor.jl")

end # module


text:
using PyPlot
using NCDatasets
using Formatting
plt = PyPlot

schemes = ["AGA2020", ]#"KSC2018", "CO2012"]

d = Dict()

println("Loading data...")
for scheme in schemes
    global d, ϵ, γ

    filename = "eps_gam_logpost2D_$(scheme).nc"
    println("Load file $(filename)")
    ds = Dataset(filename, "r")

    ϵ   = nomissing(ds["epsilon"][:], NaN) * 86400
    γ   = nomissing(ds["gamma"][:], NaN)
    post_w = transpose(nomissing(ds["post_w"][:], NaN))
    post_u = transpose(nomissing(ds["post_u"][:], NaN))
    post_uw = transpose(nomissing(ds["post_uw"][:], NaN))

    d[scheme] = Dict(
        :post_w => post_w,
        :post_u => post_u,
        :post_uw => post_uw,
    )
end

ϵ_vec  = collect(Float64, range(0.5, 1.0, length=11)) # day^-1
γ_vec = collect(Float64, range(0.0, 0.2, length=11)) # day^-1
#day_vec = collect(Float64, range(0.6, 2.0, length=7))
#day_on_ϵ_vec = 1 ./ (day_vec * 86400.0)

println("Plotting")
fig, ax = plt.subplots(1, 3, figsize=(8, 6), constrained_layout=true)

for _ax in ax
    _ax.plot([-10, 10], [-10, 10], "k--")
end

for scheme in schemes
    data = d[scheme]
    ax[1].contour(ϵ, γ, data[:post_w])
    ax[2].contour(ϵ, γ, data[:post_u])
    ax[3].contour(ϵ, γ, data[:post_uw])

end

ax[1].set_title("\$ w \$")
ax[2].set_title("\$\\vec{v}\$")
ax[3].set_title("\$w\$ and \$\\vec{v}\$")

for _ax in ax
    _ax.set_xticks(ϵ_vec)
    _ax.set_yticks(γ_vec)
    _ax.set_xlabel("\$ \\epsilon \$")
    _ax.set_ylabel("\$ \\gamma \$")

    _ax.set_xlim([0.5, 1])
    _ax.set_ylim([0, 0.2])
#    _ax.set_aspect(1.0)
end
#ax.set_xticks(day_on_ϵ_vec)
#ax.set_xticklabels([ format("{:.2f}", day) for day in day_vec])

#ax.invert_xaxis()

plt.show()

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
export Groups

function Groups(;groupName=[],groupSize=[])
  groups=Dict()
  groups[:N]=sum(groupSize)
  groups[:index]=[]
  startIndex=cumsum(groupSize)-groupSize+1
  for (idx,name) in enumerate(groupName)
    groups[name]=Dict()
    groups[name][:localindex]=startIndex[idx]:startIndex[idx]+groupSize[idx]-1
    groups[name][:processor]=idx
    push!(groups[:index],groups[name][:localindex])
  end
  return groups
end


text:
module Internal

include("base.jl")
include("array_slices.jl")
include("array.jl")

end