text:
##############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# cluster.jl
#
# From "Approximating K-means-type clustering via semidefinite programming"
# By Jiming Peng and Yu Wei
#
# Given a set of points a_1, ..., a_m  in R_n, allocate them to k clusters
#############################################################################

using JuMP
using Mosek

# Data points
n = 2
m = 6
a = Any[[2.0, 2.0], [2.5, 2.1], [7.0, 7.0],
        [2.2, 2.3], [6.8, 7.0], [7.2, 7.5]]
k = 2

# Weight matrix
W = zeros(m,m)
for i = 1:m
    for j = i+1:m
        dist = exp(-norm(a[i] - a[j])/1.0)
        W[i,j] = dist
        W[j,i] = dist
    end
end

mod = Model()

# Z >= 0, PSD
@variable(mod, Z[1:m,1:m], SDP)
@constraint(mod, Z .>= 0)

# min Tr(W(I-Z))
@objective(mod, Min, trace(W * (eye(m,m) - Z)))

# Z e = e
@constraint(mod, Z*ones(m) .== ones(m))
#for i = 1:m
#    addConstraint(mod, sum([ Z[i,j] for j = 1:m]) == 1)
#end

# Tr(Z) = k
@constraint(mod, trace(Z) == k)

solve(mod)

Z_val = getvalue(Z)[:,:]
println("Raw solution")
println(round(Z_val,4))

# A simple rounding scheme
which_cluster = zeros(Int,m)
num_clusters = 0
for i = 1:m
    Z_val[i,i] <= 1e-6 && continue

    if which_cluster[i] == 0
        num_clusters += 1
        which_cluster[i] = num_clusters
        for j = i+1:m
            if norm(Z_val[i,j] - Z_val[i,i]) <= 1e-6
                which_cluster[j] = num_clusters
            end
        end
    end
end

# Print results
for cluster = 1:k
    println("Cluster $cluster")
    for i = 1:m
        if which_cluster[i] == cluster
            println(a[i])
        end
    end
end


text:
module TestSystems
    using ..HomConBase

    export cyclic5, cyclic5Solutions, cyclic7, cyclic7Solutions

    """
        cyclical(iter, n)

    Cycles through an `iter` with n consecutive elements.

    ###Example:
        cyclical([1, 2, 3, 4], 2) == [[1, 2], [2, 3], [3, 4], [4, 1]]
    """
    function cyclical(iter, n)
        m = length(iter)
        if n == m
            return [collect(iter)]
        end
        values = collect(iter);
        res = Vector{Vector{eltype(values)}}()
        for k=0:m-1
            push!(res, [values[((k + i) % m) + 1] for i=0:n-1])
        end
        res
    end

    function monompoly(monomials::Vector{Vector{Int}}, totalvars)
        exps = zeros(Int, totalvars, length(monomials))
        for i in eachindex(monomials)
            for var in monomials[i]
                exps[var, i] = 1
            end
        end
        coeffs = ones(Int, length(monomials))
        Poly(exps, coeffs)
    end

    function cyclical_polys(n)
        F = map(k -> monompoly(cyclical(1:n, k), n), 1:n-1)
        # now we have to construct x₁x₂x₃...x_n - 1
        exps = zeros(Int, n, 2)
        exps[:,1] = 1
        push!(F, Poly(exps, [1, -1]))
        PolySystem(F, [Symbol("z$i") for i=1:n])
    end

    """
        cyclic5()

    Cyclic5 example problem folowing [^1]

    [^1]: A faster way to count the solutions of inhomogeneous systems of algebraic equations,
    with applications to cyclic n-roots
    """
    function cyclic5()
        cyclical_polys(5)
    end

    """
        cyclic5Solutions()

    Solutions to the Cyclic5 example.
    """
    function cyclic5Solutions()
        # (1,w,w^{2k},w^{3k},w^{4k}) with w = exp(2πi/5) permuted cyclical, 1 ≦ k ≦ 4
        w = exp(2.0*π*im/5.0)
        classic_sols = vec([ map(i -> w^(i*k), perm) for k in 1:4, perm in cyclical(0:4, 5)])

        #solution of ɛ^2+3ɛ+1=0
        ɛ = -1.5 - 0.5 * √5.0 + 0*im
        extended_sols = vec([ rev ? reverse!(w^k * perm) : w^k * perm for perm in cyclical([ɛ, 1/ɛ, 1, 1, 1], 5),
                                                                            k in 0:4,
                                                                            rev in [false, true]])
        vcat(classic_sols, extended_sols)
    end
    """
        cyclic7()

    Cyclic7 problem folowing [^1]

    [^1]: A faster way to count the solutions of inhomogeneous systems of algebraic equations,
    with applications to cyclic n-roots
    """
    function cyclic7()
        cyclical_polys(7)
    end

    """
        cyclic7Solutions()

    Solutions to the Cyclic7 example.
    """
    function cyclic7Solutions()
        # (1,w,w^{2k},w^{3k},w^{4k},w^{5k},w^{6k}) with w = exp(2πi/7) permuted cyclical, 1 ≦ k ≦ 6
        w = exp(2.0*π*im/7.0)
        classic_sols = vec([ map(i -> w^(i*k), perm) for k in 1:6, perm in cyclical(0:6, 7)])

        #solution of ɛ^2+5ɛ+1=0
        ɛ = -2.5 - 0.5 * √21.0 + 0*im
        extended_sols = vec([ rev ? reverse!(w^k * perm) : w^k * perm for perm in cyclical([ɛ, 1/ɛ, 1, 1, 1, 1, 1], 7),
                                                                            k in 0:6,
                                                                            rev in [false, true]])

        # index two solutions
        d = 0.25*im*(√7 + 3*im)
        index_two = vec([ rev ? reverse!(w^k * perm) : w^k * perm for perm in cyclical([d, conj(d), d, 1, conj(d), 1, 1], 7),
                                                                            k in 0:6,
                                                                            rev in [false, true]])
        index_two_conj = vec([ rev ? reverse!(w^k * perm) : w^k * perm for perm in cyclical([conj(d), d, conj(d), 1, d, 1, 1], 7),
                                                                            k in 0:6,
                                                                            rev in [false, true]])

        # index three solutions
        a = 2.738895317095 + 0im
        b = 3.436680125767 + 0im
        c = -0.1298393513967 + 0im
        sol_1 = [a, b, c, 1, 1/c, 1/b, 1/a]
        sol_2 = [a*b,c,1/(b*c),1,b*c,1/c, 1/(a*b)]
        sol_3 = [a*b*c,1/(b*c),b,1,1/b,b*c,1/(a*b*c)]

        perms = [cyclical(sol_1, 7); cyclical(sol_2, 7); cyclical(sol_3, 7)]
        index_three_first = vec([ rev ? reverse!(w^k * perm) : w^k * perm for perm in perms, k in 0:6, rev in [false, true]])

        a = exp(4.3128389787245im)
        b = exp(1.356227956787im)
        c = exp(1.900668281165im)
        sol_1 = [a, b, c, 1, 1/c, 1/b, 1/a]
        sol_2 = [a*b,c,1/(b*c),1,b*c,1/c, 1/(a*b)]
        sol_3 = [a*b*c,1/(b*c),b,1,1/b,b*c,1/(a*b*c)]
        perms = [cyclical(sol_1, 7); cyclical(sol_2, 7); cyclical(sol_3, 7)]
        index_three_second = vec([ rev ? reverse!(w^k * perm) : w^k * perm for perm in perms, k in 0:6, rev in [false, true]])

        vcat(classic_sols, extended_sols, index_two, index_two_conj, index_three_first, index_three_second)
    end
end

ext: jl


text:
mutable struct CategoricalVector{T}
    items::Vector{T}
    cdf::Vector{Float64}

    CategoricalVector{T}(item::T, weight::Float64) where T = new(T[item], Float64[weight])
end

CategoricalVector(item::T, weight::Float64) where T = CategoricalVector{T}(item, weight)

n_items(d::CategoricalVector) = length(d.items)

function insert!(c::CategoricalVector{T}, item::T, weight::Float64) where T
    push!(c.items, item)
    push!(c.cdf, c.cdf[end]+weight)
end

function rand(rng::AbstractRNG, d::CategoricalVector)
    t = rand(rng)*d.cdf[end]
    large = length(d.cdf) # index of cdf value that is bigger than t
    small = 0 # index of cdf value that is smaller than t
    while large > small + 1
        new = div(small + large, 2)
        if t < d.cdf[new]
            large = new
        else
            small = new
        end
    end
    return d.items[large]
end

function mean(d::CategoricalVector)
    sum = first(d.items)*first(d.cdf)
    for i in 2:length(d.items)
        sum += d.items[i]*(d.cdf[i]-d.cdf[i-1])
    end
    return sum/last(d.cdf)
end

function first_mean(d::CategoricalVector{Tuple{S,Float64}}) where S
    sum = first(first(d.items))*first(d.cdf)
    for i in 2:length(d.items)
        sum += first(d.items[i])*(d.cdf[i]-d.cdf[i-1])
    end
    return sum/last(d.cdf)
end

#=
function rand(rng::AbstractRNG, d::CategoricalVector)
    t = rand(rng) * d.weight_sum
    i = 1
    cw = d.weights[1]
    while cw < t && i < length(d.weights)
        i += 1
        @inbounds cw += d.weights[i]
    end
    return d.items[i]
end
=#


text:
using FSimBase
using DifferentialEquations
using ComponentArrays
using Test
using LinearAlgebra
using DataFrames

function main()
    state0 = [1.0, 2.0]
    tf = 10
    Δt = 1
    @Loggable function dynamics!(dx, x, p, t)
        @log t
        @log x
        dx .= 0.99*x
    end
    simulator = Simulator(
                          state0, dynamics!;
                          tf=tf,
                          Problem=DiscreteProblem,
                          solver=FunctionMap(),
                         )
    # solve approach (automatically reinitialised)
    @time _df = solve(simulator; savestep=Δt)
    @show _df
    # interactive simulation
    ## step_until!
    df = DataFrame()
    reinit!(simulator)
    ts = 0:Δt:tf
    @time for t in ts
        step_until!(simulator, t)
        push!(simulator, df)  # flag=true is default
        # safer way:
        # flag = step_until!(simulator, t)
        # push!(simulator, df, flag)
        # or, equivalently,
        # push!(simulator, df, step_until!(simulator, t))  # compact form
    end
    @test df == _df
end


text:
using MLStyle

rmlines = @λ begin
    e :: Expr           -> Expr(e.head, filter(x -> x !== nothing, map(rmlines, e.args))...)
      :: LineNumberNode -> nothing
    a                   -> a
end

macroexpand(@__MODULE__, :(@match 1 begin
        1 => 1
end)) |> rmlines

@testset "call / expr" begin
@test @match Expr(:call, :f, :a) begin
    Expr(:call, tail...) => collect(tail) == [:f, :a]
end
end

@testset "call / ast" begin
ast = :(f(a, b))

@test @match Expr(:call, :f, :a, :b) begin
  :($f($a, $b)) => (a, b) == (:a, :b)
end
end

# https://github.com/JuliaLang/julia/pull/35138
const LN_OFF = VERSION >= v"1.5.0" ? 1 : 0

@testset "Ast Pattern" begin
instance = :(function f(a, b, c)
                   a + b + c
             end)
@test @match instance begin
    Expr(:function, Expr(:call, funcname, args...), block) => begin
        (funcname, collect(args)) == (:f, [:a, :b, :c]) &&
        block.head == :block                            &&
        block.args[1 + LN_OFF] isa LineNumberNode                &&
        block.args[2 + LN_OFF] ==  :(a + b + c)
    end
end

@testset "function" begin
@test @match instance begin
        :(function $funcname($(args...)) $(block...) end) => begin
        (funcname, collect(args)) == (:f, [:a, :b, :c]) &&
        block[1 + LN_OFF] isa LineNumberNode                     &&
        block[2 + LN_OFF] ==  :(a + b + c)
    end
end
end

let_expr = :(let a = 10 + 20, b = 20
              20a
             end)

@testset "let binding" begin

    @test @match let_expr begin
            :(let $bind_name = $fn($left, $right), $(other_bindings...)
                $(block...)
            end) => begin
            bind_name == :a   &&
            fn        == :(+) &&
            left      == 10   &&
            right     == 20   &&
            block[1]  isa LineNumberNode &&
            block[2]  == :(20a)
        end
    end
    end
end

@testset "only head" begin
    @test @match Expr(:f) begin
        Expr(a) => a === :f
        _       => false
    end
end

@testset "only pack" begin
    @test @match Expr(:f) begin
        Expr(a...) => a == [:f]
        _          => false
    end
end

@testset "pack all" begin
    @test @match Expr(:f, :a) begin
        Expr(a...) => a == [:f, :a]
        _          => false
    end
end

module ASTSamples

    node_fn1 = :(function f(a, b) a + b end)
    node_fn2 = :(function f(a, b, c...) c end)
    node_let = :(let x = a + b
                2x
            end)
    node_chain = :(subject.method(arg1, arg2))
    node_struct = :(
            struct name <: base
                field1 :: Int
                field2 :: Float32
            end
    )

    node_const = :(
            const a = value
    )
    node_assign = :(a = b + c)

end

@testset "case from Match.jl" begin

extract_name = @λ begin
        e ::Symbol                         -> e
        Expr(:<:, a, _)                    -> extract_name(a)
        Expr(:struct, _, name, _)          -> extract_name(name)
        Expr(:call, f, _...)               -> extract_name(f)
        Expr(:., subject, attr, _...)      -> extract_name(subject)
        Expr(:function, sig, _...)         -> extract_name(sig)
        Expr(:const, assn, _...)           -> extract_name(assn)
        Expr(:(=), fn, body, _...)         -> extract_name(fn)
        Expr(expr_type,  _...)             -> error("Can't extract name from ",
                                                    expr_type, " expression:\n",
                                                    "    $e\n")
end

@test extract_name(ASTSamples.node_fn1) == :f
@test extract_name(ASTSamples.node_fn2) == :f
@test extract_name(ASTSamples.node_chain) == :subject
@test extract_name(ASTSamples.node_struct) == :name
@test extract_name(ASTSamples.node_const) == :a
@test extract_name(ASTSamples.node_assign) == :a
@test_skip extract_name(:(1 + 1))

end


text:
using SgdlDev
using Test

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


text:
import REPL.TerminalMenus
import REPL.Terminals

struct MenuItem{V}
    name::String
    label::String
    value::V
end

name(m::MenuItem) = m.name
label(m::MenuItem) = isempty(m.label) ? m.name : m.label
value(m::MenuItem) = m.value
text(m::MenuItem) = begin
    l = label(m)
    w = Terminals.width(TerminalMenus.terminal)
    #HACK: length(l) contains escape sequences, so actual line may look shorter
    v = repr(value(m); context=:maxlength => w - length(l))
    isempty(l) ? v : "$l $(misc_color("=")) $v"
end

title(m::MenuItem{<:System}) = "$(system_color(name(m)))"
title(m::MenuItem{<:State}) = "$(state_color(name(m)))"
title(m::MenuItem) = "$(non_state_color(name(m)))"
title(t::Vector{<:MenuItem}) = begin
    sep = " $(misc_color(">")) "
    join(title.(t), sep)
end

dive(s::System, t) = dive(map(zip(fieldnamesalias(s), s)) do ((n, a), v)
    k = string(n)
    cn = string(uncanonicalname(n))
    ca = string(uncanonicalname(a))
    l = "$(variable_color(cn))"
    !isnothing(a) && (l *= " $(misc_color("($ca)"))")
    MenuItem(k, l, v)
end, t)
dive(s::Vector{<:System}, t) = dive(map(t -> MenuItem(string(t[1]), "", t[2]), enumerate(s)), t)
dive(s::State{<:System}, t) = dive([MenuItem("1", "", s')], t)
dive(s::State{<:Vector}, t) = dive(map(t -> MenuItem(string(t[1]), "", t[2]), enumerate(s')), t)
dive(s::State, t) = dive(t) do io
    look(io, value(t[end-1]), Symbol(name(t[end])))
end
dive(l::Vector{<:MenuItem}, t) = begin
    isempty(l) && return
    term = TerminalMenus.terminal
    o = term.out_stream
    i = 1
    while true
        println(o, title(t))
        N = length(l)
        M = TerminalMenus.RadioMenu(text.(l); charset=:unicode, pagesize=N)
        i = TerminalMenus.request(M; cursor=i)
        n = min(N, M.pagesize)
        #HACK: for single option menu?
        n == 1 && (n += 1)
        print(o, repeat("\x1b[9999D\x1b[1A", n+1)) # move up
        print(o, "\x1b[J") # clear lines below
        if i > 0
            v = l[i]
            dive(value(v), [t..., v])
        else
            break
        end
    end
end
dive(m::MenuItem, t) = dive(value(m), t)
dive(v, t) = dive(t) do io
    show(IOContext(io, :limit => true), MIME("text/plain"), v)
end
dive(f::Function, t) = begin
    term = TerminalMenus.terminal
    i = term.in_stream
    o = term.out_stream
    b = IOBuffer()
    #HACK: dive() assumes color terminal
    x = IOContext(b, :color => get(o, :color, true))
    println(x, title(t))
    f(x)
    println(x)
    n = countlines(seekstart(b))
    print(o, String(take!(b)))
    Terminals.raw!(term, true) && print(o, "\x1b[?25l") # hide cursor
    c = TerminalMenus.readkey(i)
    Terminals.raw!(term, false) && print(o, "\x1b[?25h") # unhide cursor
    print(o, repeat("\x1b[9999D\x1b[1A", n)) # move up
    print(o, "\x1b[J") # clear lines below
    if c == 13 # enter
        throw(value(t[end]))
    elseif c == 3 # ctrl-c
        throw(InterruptException())
    end
end

"""
    dive(s)

Inspect an instance of system `s` by navigating hierarchy of variables displayed in a tree structure.

Pressing up/down arrow keys allows navigation. Press 'enter' to dive into a deeper level and press 'q' to come back. A leaf node of the tree shows an output of `look` regarding the variable. Pressing 'enter' again would return a variable itself and exit to REPL.

Only works in a terminal environment; not working on Jupyter Notebook.

See also: [`look`](@ref)

# Arguments
- `s::System`: instance of target system.

# Examples
```julia-repl
julia> @system S(Controller) begin
           a => 1 ~ preserve(parameter)
       end;

julia> s = instance(S);

julia> dive(s)
S
 → context = <Context>
   config = <Config>
   a = 1.0
```
"""
dive(s::System) = begin
    if isdefined(Main, :IJulia) && Main.IJulia.inited
        return look(s)
    end
    try
        dive(s, [MenuItem(string(namefor(s)), "", s)])
    catch e
        !isa(e, InterruptException) && return e
    end
    nothing
end

export dive


text:
export ss, ss2tf

# State Space Constructor
# =======================

"""
### function StateSpace(A::Matrix, B::Matrix, C::Matrix, D::Matrix)

Initialize a state space with the following form:

x = Ax + Bu
y = Cx + Du

##### Args

* A: Matrix A.
* B: Matrix B.
* C: Matrix C.
* D: Matrix D.

##### Returns

The new state space.

"""

function ss(A::Array, B::Array, C::Array, D::Array)
    StateSpace(A,B,C,D)
end

# Print Functions
# ===============

"""
### function Base.show(io::IO, Sys::StateSpace)

Print space state.

##### Args

* io: IO in which the state space will be printed.
* Sys: Space state.

"""

function Base.show(io::IO, Sys::StateSpace)
    print(io, "\n")
    print(io, "A = ")
    print(io, Sys.A, "\n\n")
    print(io, "B = ")
    print(io, Sys.B, "\n\n")
    print(io, "C = ")
    print(io, Sys.C, "\n\n")
    print(io, "D = ")
    print(io, Sys.D, "\n\n")
end

Base.print(io::IO, Sys::StateSpace) = Base.show(io, Sys)

# Conversion from State Space to Transfer Function
# ================================================

"""
### function ss2tf(sys::StateSpace)

Convert a state space representation of a system to a transfer function.

##### Args

* sys: State space representation.

##### Returns

The transfer function.

##### Remarks

The output will always be of type `Float64` due to the operations to convert
between the representations.

"""

function ss2tf(sys::StateSpace)
    ( !issiso(sys) ) && error("ss2tf() is only supported for SISO systems.")

    # Get the invariant zeros of the system.
    z = izero(sys)

    # Get the poles of the system.
    p = pole(sys)

    # Get the gain of the system.
    K = gain(sys)

    # Compute the transfer function.
    zpk(z, p, K)
end


text:
using WebDriver

ENV["WEBDRIVER_HOST"] = get(ENV, "WEBDRIVER_HOST", "localhost")
ENV["WEBDRIVER_PORT"] = get(ENV, "WEBDRIVER_PORT", "4444")

capabilities = Capabilities("chrome")
wd = RemoteWebDriver(
    capabilities,
    host = ENV["WEBDRIVER_HOST"],
    port = parse(Int, ENV["WEBDRIVER_PORT"]),
    )
# New Session
session = Session(wd)

navigate!(session, "https://sutra.oslpr.org/osl/esutra/MedidaBus.aspx")
using Base64
using Base64
ss = write(joinpath(@__DIR__, "img.png"), base64decode(screenshot(session)))

#how about?
fi = "c:\\temp\\myfile.png"
function screenshot(fi,session)
    ss = write(fi, base64decode(screenshot(session)))
    return ss
end
using Base64
ss = write(joinpath("img.png"), base64decode(screenshot(session)))

max_stars_repo_name: Nosferican/PRLegVotos

max_issues_repo_name: Nosferican/PRLegVotos

max_forks_repo_name: Nosferican/PRLegVotos


text:
struct Resource
    data::Any
    name::Symbol
    persistent::Bool
end

"""
Frame-global structure that holds all data needed in the frame.

Its linear allocator is used for allocating lots of small objects, like material parameters and vertex data.

Other resources that require a global descriptor set (bind-once strategy) are put into a `ResourceDescriptors`.
This includes general image data & samplers, with the corresponding descriptors.

The index list is used to append index data and is turned into an index buffer before initiating the render sequence.
"""
struct GlobalData
    allocator::LinearAllocator
    resources::ResourceDescriptors
    index_list::Vector{UInt32}
    index_buffer::Ref{BufferBlock{MemoryBlock}}
end

GlobalData(device) = GlobalData(
    LinearAllocator(device, 1_000_000), # 1 MB
    ResourceDescriptors(device),
    [],
    Ref{BufferBlock{MemoryBlock}}(),
)

function populate_descriptor_sets!(gd::GlobalData)
    state = gd.resources.gset.state
    types = [Vk.DESCRIPTOR_TYPE_SAMPLED_IMAGE, Vk.DESCRIPTOR_TYPE_STORAGE_IMAGE, Vk.DESCRIPTOR_TYPE_SAMPLER, Vk.DESCRIPTOR_TYPE_COMBINED_IMAGE_SAMPLER]
    writes = types |> Enumerate() |> Map() do (i, type)
        infos = state[type]
        !isempty(infos) || return nothing
        Vk.WriteDescriptorSet(gd.resources.gset.set.handle, i-1, 0, type, infos, [], [])
    end |> Filter(!isnothing) |> collect
    !isempty(writes) || return
    Vk.update_descriptor_sets(
        device(gd.resources.gset.set),
        writes,
        [],
    )
end

function allocate_index_buffer(gd::GlobalData, device::Device)
    gd.index_buffer[] = buffer(device, convert(Vector{UInt32}, gd.index_list .- 1); usage = Vk.BUFFER_USAGE_INDEX_BUFFER_BIT)
end

struct Frame
    resources::Dictionary{Symbol,Resource}
    gd::GlobalData
end

Frame(device) = Frame(Dictionary(), GlobalData(device))

function register(frame::Frame, resource_name::Symbol, resource_data; persistent = true)
    insert!(frame.resources, resource_name, Resource(resource_data, resource_name, persistent))
end

function Base.get(frame::Frame, ::Type{T}, symbol::Symbol) where {T}
    frame.resources[symbol].data::T
end

#=

struct FrameSynchronization
    image_acquired::Vk.Semaphore
    image_rendered::Vk.Semaphore
    has_rendered::Vk.Fence
end

FrameSynchronization(device) = FrameSynchronization(Vk.Semaphore(device), Vk.Semaphore(device), Vk.Fence(device; flags = FENCE_CREATE_SIGNALED_BIT))

struct FrameState
    device::Device
    swapchain::Ref{Swapchain}
    frames::Vector{Frame}
    current_frame::Ref{Frame}
    frame_count::Ref{Int}
    syncs::Dictionary{Frame,FrameSynchronization}
end

device(fs::FrameState) = fs.device

function FrameState(device, swapchain::Swapchain)
    max_in_flight = info(swapchain).min_image_count
    fs = FrameState(device, Ref(swapchain), [], Ref{Frame}(), Ref(0), Dictionary())
    update!(fs)
end

function Vk.SurfaceCapabilitiesKHR(fs::FrameState)
    unwrap(get_physical_device_surface_capabilities_khr(device(fs).physical_device, info(fs.swapchain[]).surface))
end

function recreate_swapchain!(fs::FrameState, new_extent::NTuple{2,Int})
    swapchain = fs.swapchain[]
    swapchain_info = setproperties(info(swapchain), old_swapchain = handle(swapchain), image_extent = Vk.Extent2D(new_extent...))
    swapchain_handle = unwrap(create_swapchain_khr(device(fs), swapchain_info))
    fs.swapchain[] = Created(swapchain_handle, swapchain_info)
    fs
end

function get_frames(fs::FrameState)
    swapchain = fs.swapchain[]
    extent = info(swapchain).image_extent
    map(unwrap(get_swapchain_images_khr(device(fs), swapchain))) do img
        view = View(
            device(fs),
            img,
            IMAGE_VIEW_TYPE_2D,
            info(swapchain).image_format,
            ComponentMapping(fill(COMPONENT_SWIZZLE_IDENTITY, 4)...),
            ImageSubresourceRange(IMAGE_ASPECT_COLOR_BIT, 0, 1, 0, 1),
        )

        fb = Framebuffer(device(fs), fs.render_pass, [view], extent.width, extent.height, 1)
        Frame(img, view, fb)
    end
end

function update!(fs::FrameState)
    empty!(fs.frames)
    empty!(fs.syncs)

    # TODO: make fields of returned-only structs high-level in Vulkan.jl
    _extent = SurfaceCapabilitiesKHR(fs).current_extent.vks
    new_extent = Extent2D(_extent.width, _extent.height)
    if new_extent ≠ info(fs.swapchain[]).image_extent
        recreate_swapchain!(fs, new_extent)
    end

    append!(fs.frames, get_frames(fs))
    fs.current_frame[] = first(fs.frames)

    foreach(fs.frames) do frame
        insert!(fs.syncs, frame, FrameSynchronization(device(fs.render_pass)))
    end
    fs
end

"""
    command_buffers(frame_state, app)

Command buffers generated by the application.
"""
function command_buffers end

function next_frame!(fs::FrameState, dispatch::QueueDispatch, app)
    swapchain = fs.swapchain[]

    # acquire next image
    image_acquired = fs.syncs[fs.current_frame[]].image_acquired
    status = @timeit to "Acquire next image" acquire_next_image_khr(device(fs), swapchain, typemax(UInt64); semaphore = image_acquired)
    if !iserror(status)
        idx, result = unwrap(status)
        @assert result in (SUCCESS, SUBOPTIMAL_KHR) "$result: Could not retrieve next swapchain image"
        frame = fs.frames[idx + 1]
        sync = fs.syncs[frame]

        # submit rendering commands
        @timeit to "Submit rendering commands" begin
            cbuffs = @timeit to "Create command buffers" command_buffers(fs, app)
            img_acquired_info = SemaphoreSubmitInfoKHR(image_acquired, 0, 0; stage_mask = PIPELINE_STAGE_2_COLOR_ATTACHMENT_OUTPUT_BIT_KHR)
            img_rendered_info = SemaphoreSubmitInfoKHR(sync.has_rendered, 0, 0; stage_mask = PIPELINE_STAGE_2_COLOR_ATTACHMENT_OUTPUT_BIT_KHR)
            render_info = SubmitInfo2KHR([img_acquired_info], CommandBufferSubmitInfoKHR.(cbuffs, 0), [img_rendered_info])
            @timeit to "Wait for frame being rendered" wait_for_fences(device(fs), [sync.has_rendered], false, typemax(UInt64))
            reset_fences(device(fs), [sync.has_rendered])
            submit(dispatch, QUEUE_GRAPHICS_BIT, [render_info]; fence = sync.has_rendered)
        end

        # submit presentation commands
        @timeit to "Submit presentation commands" begin
            present_info = PresentInfoKHR([sync.image_rendered], [swapchain], [idx])
            unwrap(present(dispatch, present_info))
        end

        fs.frame_count[] += 1
        fs.current_frame[] = frame
    else
        err = unwrap_error(status)
        if err.code == ERROR_OUT_OF_DATE_KHR
            # recreate swapchain and start over
            @timeit to "Recreate swapchain" update!(fs)
            next_frame!(fs, rdr, app)
        end
    end
end

function wait_hasrendered(fs::FrameState)
    wait_for_fences(device(fs), getproperty.(fs.syncs, :has_rendered), true, typemax(UInt64))
end

function next_frame!(ctx::FrameContext)
    next_frame!(ctx.state)
    foreach(ctx.pending) do resource
        finalize(resource)
    end
    empty!(ctx.pending)
end

function stage!(ctx::FrameContext, resource)
    push!(ctx.pending, resource)
end

=#


text:
#=
    The File contains `BloomSparsik<T>` and related funcs implementation

    `BloomSparsik` implements `AbstractSparseObject` interface
=#

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

include("../typics.jl")
include("../gizmos.jl")

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

import Base: ==, !=, +, -, *
import Nemo: GaloisField, GaloisFmpzField, fmpq, ZZ
import Distributions: Bernoulli
import AbstractAlgebra: elem_type, Field, FieldElem,
                        characteristic, FracField
import StaticArrays: SA

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

_powersof2 = SA[UInt64(1), UInt64(2), UInt64(4), UInt64(8), UInt64(16),
                UInt64(32), UInt64(64), UInt64(128), UInt64(256), UInt64(512),
                UInt64(1024), UInt64(2048), UInt64(4096), UInt64(8192), UInt64(16384),
                UInt64(32768), UInt64(65536), UInt64(131072), UInt64(262144), UInt64(524288),
                UInt64(1048576), UInt64(2097152), UInt64(4194304), UInt64(8388608), UInt64(16777216),
                UInt64(33554432), UInt64(67108864), UInt64(134217728), UInt64(268435456),
                UInt64(536870912), UInt64(1073741824), UInt64(2147483648), UInt64(4294967296),
                UInt64(8589934592), UInt64(17179869184), UInt64(34359738368), UInt64(68719476736),
                UInt64(137438953472), UInt64(274877906944), UInt64(549755813888), UInt64(1099511627776),
                UInt64(2199023255552), UInt64(4398046511104), UInt64(8796093022208), UInt64(17592186044416),
                UInt64(35184372088832), UInt64(70368744177664), UInt64(140737488355328), UInt64(281474976710656),
                UInt64(562949953421312), UInt64(1125899906842624), UInt64(2251799813685248), UInt64(4503599627370496),
                UInt64(9007199254740992), UInt64(18014398509481984), UInt64(36028797018963968), UInt64(72057594037927936),
                UInt64(144115188075855872), UInt64(288230376151711744), UInt64(576460752303423488), UInt64(1152921504606846976),
                UInt64(2305843009213693952), UInt64(4611686018427387904), UInt64(9223372036854775808)]

_defaultbloom(x) = _powersof2[(x % 64) + 1]

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

mutable struct BloomSparsik{T<:Field} <: AbstractSparseVector{T}
    dim::Int
    field::T
    nonzero::Vector{Int}
    data::Dict{Int, <:FieldElem}
    nnzbloom::UInt64
    funcbloom
end

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

# returns the density of `v`,
# i.e the percentage of nonzero-valued elements in `v`
density(v::BloomSparsik) = length(v.nonzero) / v.dim

# returns the number of nonzero-valued elements in `v`
nnz(v::BloomSparsik) = length(v.nonzero)

# ground field!
Nemo.base_ring(v::BloomSparsik) = v.field

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

# hitherto `Sparsik(other::Sparik)` has stood for the `deepcopy`,
# from now on we override Base.deepcopy
function Base.deepcopy_internal(x::BloomSparsik, stackdict::IdDict)
    y = Sparsik(x.dim, x.field,
            Base.deepcopy_internal(x.nonzero, stackdict),
            Base.deepcopy_internal(x.data, stackdict),
            nnzhash, funcbloom)
    stackdict[x] = y
    return y
end

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

# returns the index of first nonzero or -1 if not present
# O(1)
function first_nonzero(v::BloomSparsik)
    if iszero(v)
        return -1
    end
    return first(v.nonzero)
end

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

# returns v * c
# does modify `v`
# if `v` has k nonzeroes
# O(k)
function scale!(v::BloomSparsik, c)
    if iszero(base_ring(v)(c))
        return empty!(v)
    end

    for idx in v.nonzero
        v.data[idx] *= c
    end
    return v
end

# returns v * c
function scale(v::BloomSparsik, c)
    return scale!(deepcopy(v), c)
end

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

# returns v + c * u
# mutates v
# if the vectors `v` and `u` have k and r nonzeroes respectively
# O(k + r)
function reduce!(v::BloomSparsik{T}, u::BloomSparsik{T}, c) where {T}
    new_nonzero = Int[]
    i, j = 1, 1

    v_nnz = deepcopy(v.nonzero)
    u_nnz = u.nonzero
    empty!(v.nonzero)

    while i <= length(v_nnz) || j <= length(u_nnz)
        new_idx = 0

        if i > length(v_nnz)
            new_idx = u_nnz[j]
            j += 1
        elseif j > length(u_nnz)
            new_idx = v_nnz[i]
            i += 1
        else
            new_idx = min(v_nnz[i], u_nnz[j])
            if v_nnz[i] < u_nnz[j]
                i += 1
            elseif v_nnz[i] > u_nnz[j]
                j += 1
            else
                i += 1
                j += 1
            end
        end
        v.data[new_idx] = v[new_idx] + c * u[new_idx]
        if !iszero(v[new_idx])
            push!(v.nonzero, new_idx)
            v.nnzbloom |= v.funcbloom(new_idx)
        else
            delete!(v.data, new_idx)
        end
    end

    return v
end

# returns v + c * u
# O(k + r)
function Base.reduce(v::BloomSparsik{T}, u::BloomSparsik{T}, c) where {T}
    return reduce!(deepcopy(v), u, c)
end

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

# returns standard inner product of v and u
# if the vectors `v` and `u` have k and r nonzeroes respectively
# O(min(k, r))
#
# This should not be called in case of non-overlapping
function inner(v::BloomSparsik{T}, u::BloomSparsik{T}) where {T}
    ans = zero(v.field)

    if length(v) > length(u)
        inner(u, v)
    end

    for i in v.nonzero
        # zero values can be skipped instantly
        # TODO: fix
        if haskey(u.data, i)
            ans += v.data[i] * u.data[i]
        end
    end

    return ans
end

function inner(A::BloomSparsik{T}, c::FieldElem) where {T}
    if iszero(c)
        return zero(c)
    end

    error("not implemented!")
end

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

# constructs a `BloomSparsik` instance from a dense Vector,
# the provided field and the function for bloom hashing
# each `a` element is explicitly converted to its field representation
# O(n) where n = length(a)
function from_dense(a::AbstractVector, field, funcbloom)
    new_nonzero = Int[]
    nnzbloom = UInt64(0)
    new_data = Dict{Int, elem_type(field)}()

    for idx in findall(!(iszero ∘ field), a)
        push!(new_nonzero, idx[1])
        nnzbloom |= funcbloom(idx)
        new_data[idx[1]] = field(a[idx])
    end

    return BloomSparsik(length(a), field, new_nonzero, new_data, nnzbloom, funcbloom)
end

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

# empties the `v`
# if the vector `v` has k nonzeroes
# O(k)
function Base.empty!(v::BloomSparsik)
    empty!(v.nonzero)
    empty!(v.data)
    v.nnzbloom = zero(v.nnzbloom)
    return v
end

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

# function for Sparsiks compatibility
# O(1)
function zero_sparsik(dim, field, funcbloom)
    return BloomSparsik(dim, field, Int[], Dict{Int, elem_type(field)}(), UInt64(0), funcbloom)
    # return zero(field)
end

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

Base.zero(v::BloomSparsik) = zero_sparsik(v.dim, v.field, _default_bloom)
Base.iszero(v::BloomSparsik{T}) where {T} = length(v) == 0

Base.get(v::BloomSparsik, i::Int) = get(v.data, i, zero(v.field))

function Base.haskey(A::BloomSparsik, i::Int)
    return haskey(A.data, i)
end

Base.getindex(v::BloomSparsik, i::Int) = get(v, i)

Base.size(v::BloomSparsik) = (v.dim, )
Base.size(v::BloomSparsik, i::Int) = size(v)[i]
Base.length(v::BloomSparsik) = length(v.nonzero)
dim(v::BloomSparsik) = v.dim

Base.isempty(v::BloomSparsik) = length(v.nonzero) == 0

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

# Strategy of types promotion:
# explicit restriction to the operands types
+(v::BloomSparsik{T}, u::BloomSparsik{T}) where {T} = reduce(v, u, one(v.field))
-(v::BloomSparsik{T}, u::BloomSparsik{T}) where {T} = reduce(v, u, - one(v.field))
-(v::BloomSparsik) = scale(v, - one(v.field))

*(v::BloomSparsik, c::Number) = scale(v, c)
*(c::Number, v::BloomSparsik) = v * c

==(v::BloomSparsik{T}, u::BloomSparsik{T}) where {T} = (v.dim == u.dim &&
            v.field == u.field &&
            v.data == u.data &&
            v.nonzero == u.nonzero)
!=(v::BloomSparsik{T}, u::BloomSparsik{T}) where {T} = !(v == u)

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

Base.repr(::MIME"text/plain", A::BloomSparsik) = "bloom: $(A.nnzbloom)\n($(join(map(idx -> A[idx], 1 : A.dim), ", ")))"
Base.show(io::IO, A::BloomSparsik) = print(io, repr(MIME("text/plain"), A))

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

# iterates over all nonzeroes in `v`
# each returned element is a pair (idx, value)
# representing a nonzero element from `v`
# linear in number of nonzeroes
function Base.iterate(v::BloomSparsik, state=1)
    if state <= length(v)
        idx = v.nonzero[state]
        state += 1
        return ((idx, v[idx]), state)
    else
        return nothing
    end
end

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

Base.valtype(v::BloomSparsik) = elem_type(v.field)
Base.eltype(v::BloomSparsik) = (Int, valtype(v))

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

# Returns a new Sparsik object consisting of elements
# from the given `v` each converted to the `field`
# Throws if any of the denominators reduces to zero
#
# O(k) if k is the number of nonzero in `v`
function modular_reduction(v::BloomSparsik, field)
    new_nonzero = Int[]
    new_data = Dict{Int, elem_type(field)}()
    nnzbloom = UInt64(0)

    for i in v.nonzero
        # wrapped into a function
        y = modular_reduction(v.data[i], field)
        if !iszero(y)
            push!(new_nonzero, i)
            nnzbloom |= v.funcbloom(i)
            new_data[i] = y
        end
    end

    return BloomSparsik(dim(v), field, new_nonzero, new_data, nnzbloom, v.funcbloom)
end

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

# returns a new Sparsik object consisting of elements
# from the `v` each reconstructed from v.field to QQ
#
# O†(k) if k is the number of nonzero in `v`
# † reconstruction cost not included
function rational_reconstruction(v::BloomSparsik)
    new_nonzero = Int[]
    new_data = Dict{Int, fmpq}()
    ground = base_ring(v)
    nnzbloom = UInt64(0)

    tmp_type = BigInt
    if isa(ground, GaloisField)
        tmp_type = Int
    end
    ch = convert(tmp_type, characteristic(ground))

    for (i, x) in v
        y = rational_reconstruction(convert(tmp_type, x), ch)
        if !iszero(y)
            push!(new_nonzero, i)
            nnzbloom |= v.funcbloom(i)
            new_data[i] = y
        end
    end

    return BloomSparsik(dim(v), QQ, new_nonzero, new_data, nnzbloom, v.funcbloom)
end

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

function unit_sparsik(dim, i, field, funcbloom)
    return BloomSparsik(dim..., field, [i], Dict(i => one(field)), _powersof2[i], funcbloom)
end

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

# returns a Sparsik object consisting of O(sz * density) nnz entries,
# each entry is generated uniformly and independently
function random_sparsik(sz::Int, field::T, funcbloom; density=0.1) where {T<:Field}
    λ = density * sz

    nnz = Int[]
    nnzbloom = UInt64(0)
    data = Dict{Int, elem_type(field)}()
    for _ in 1:λ
        while (idx = rand(1:sz)) != nothing
            !haskey(data, idx) && break
        end
        X = rand(field)
        push!(nnz, idx)
        nnzbloom |= funcbloom(idx)
        data[idx] = X
    end

    sort!(nnz)

    return BloomSparsik(sz, field, nnz, data, nnzbloom, funcbloom)
end

function random_sparsik(sz::Tuple{Int}, field::T, funcbloom; density=0.1) where {T<:Field}
    random_sparsik(sz..., field, funcbloom, density=density)
end

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

# Inserts `value` to the `idx`-th position in `v`
# If idx is already present, rewrites it
# if k = length(v)
# O(klogk)
#
# it is not used
function Base.insert!(v::BloomSparsik, idx::Int, value)
    if iszero(value)
        return v
    end
    if iszero(v[idx])
        pos = searchsortedfirst(v.nonzero, idx)
        insert!(v.nonzero, pos, idx)
    end
    v.data[idx] = value
    return v
end

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

# Appends `value` to the vector treated as an array
# and increments the dimension of the sparsik by one
# O(1)
function Base.append!(v::BloomSparsik, value)
    v.dim += 1
    if iszero(value)
        return v
    end

    v.data[v.dim] = value
    push!(v.nonzero, v.dim)
    nnzbloom |= funcbloom(v.dim)

    return v
end

# Erases the value at the last vector index
# and decrements the dimension of the sparsik by one
# Throws if v is 0-dim
# O(1)
function Base.pop!(v::BloomSparsik)
    if v.dim == 0
        error("sparsik $v is very very empty")
    end

    if haskey(v.data, v.dim)
        pop!(v.nonzero)
        nnzbloom ^= funcbloom(v.dim)
    end
    delete!(v.data, v.dim)
    v.dim -= 1
    return v
end

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

function restrict(v::BloomSparsik, coords)
    field = base_ring(v)
    nonzero = Int[]
    data = Dict{Int, elem_type(field)}()
    nnzbloom = UInt64(0)

    for (i, col) in enumerate(coords)
        if ! iszero(v[col])
            push!(nonzero, i)
            data[i] = v[col]
            v.nnzbloom |= _powerof2[idx]
        end
    end

    return BloomSparsik(length(coords), field, nonzero, data)
end

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

function to_dense(A::BloomSparsik)
    ans = fill(zero(base_ring(A)), size(A))
    for (i, x) in A
        ans[i] = x
    end
    return ans
end

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

# v[idx] += c
# mutates v
# O(k)
# ( but O(1) generically )
#
function reduce!(v::BloomSparsik{T}, idx::Int, c) where {T}
    if iszero(c)
        return v
    end

    if haskey(v.data, idx)
        v.data[idx] += c
        if !iszero(v.data[idx])

        else
            v.nnzbloom ^= _powerof2[idx]
            popat!(v.nonzero, idx)
            delete!(v.data, idx)
        end
    else
        v.data[idx] = c
        pos = searchsortedfirst(v.nonzero, idx)
        v.nnzbloom |= _powerof2[idx]
        insert!(v.nonzero, pos, idx)
    end
    return v
end


text:
#=
Author: Andrew Ning

A general blade element momentum (BEM) method for propellers/fans and turbines.

Some unique features:
- a simple yet very robust solution method ideal for use with optimization
- designed for compatibility with algorithmic differentiation tools
- allows for arbitrary inflow conditions, including reversed flow, hover, etc.
- convenience methods for common wind turbine inflow scenarios

=#

module CCBlade

import FLOWMath

export Rotor, Section, OperatingPoint, Outputs
export simple_op, windturbine_op
export solve, thrusttorque, nondim

include("cs_safe.jl")

include("airfoils.jl")  # all the code related to airfoil data

# --------- structs -------------

"""
    Rotor(Rhub, Rtip, B; precone=0.0, turbine=false,
        mach=nothing, re=nothing, rotation=nothing, tip=PrandtlTipHub())

Parameters defining the rotor (apply to all sections).

**Arguments**
- `Rhub::Float64`: hub radius (along blade length)
- `Rtip::Float64`: tip radius (along blade length)
- `B::Int64`: number of blades
- `precone::Float64`: precone angle
- `turbine::Bool`: true if using wind turbine conventions
- `mach::MachCorrection`: correction method for Mach number
- `re::ReCorrection`: correction method for Reynolds number
- `rotation::RotationCorrection`: correction method for blade rotation
- `tip::TipCorrection`: correction method for hub/tip loss
"""
struct Rotor{TF, TI, TB,
        T1 <: Union{Nothing, MachCorrection}, T2 <: Union{Nothing, ReCorrection},
        T3 <: Union{Nothing, RotationCorrection}, T4 <: Union{Nothing, TipCorrection}}
    Rhub::TF
    Rtip::TF
    B::TI
    precone::TF
    turbine::TB
    mach::T1
    re::T2
    rotation::T3
    tip::T4
end

# convenience constructor with keyword parameters
Rotor(Rhub, Rtip, B; precone=0.0, turbine=false, mach=nothing, re=nothing, rotation=nothing, tip=PrandtlTipHub()
    ) = Rotor(Rhub, Rtip, B, precone, turbine, mach, re, rotation, tip)

"""
    Section(r, chord, theta, af)

Define sectional properties for one station along rotor

**Arguments**
- `r::Float64`: radial location along blade
- `chord::Float64`: corresponding local chord length
- `theta::Float64`: corresponding twist angle (radians)
- `af::Function or AFType`: if function form is: `cl, cd = af(alpha, Re, Mach)`
"""
struct Section{TF1, TF2, TF3, TAF}
    r::TF1  # different types b.c. of dual numbers.  often r is fixed, while chord/theta vary.
    chord::TF2
    theta::TF3
    af::TAF
end

# convenience function to access fields within an array of structs
function Base.getproperty(obj::Vector{Section{TF1, TF2, TF3, TAF}}, sym::Symbol) where {TF1, TF2, TF3, TAF}
    return getfield.(obj, sym)
end # This is not always type stable b/c we don't know if the return type will be float or af function.

"""
    OperatingPoint(Vx, Vy, rho; pitch=0.0, mu=1.0, asound=1.0)

Operation point for a rotor.
The x direction is the axial direction, and y direction is the tangential direction in the rotor plane.
See Documentation for more detail on coordinate systems.
`Vx` and `Vy` vary radially at same locations as `r` in the rotor definition.

**Arguments**
- `Vx::Float64`: velocity in x-direction along blade
- `Vy::Float64`: velocity in y-direction along blade
- `pitch::Float64`: pitch angle (radians)
- `rho::Float64`: fluid density
- `mu::Float64`: fluid dynamic viscosity (unused if Re not included in airfoil data)
- `asound::Float64`: fluid speed of sound (unused if Mach not included in airfoil data)
"""
struct OperatingPoint{TF1, TF2, TF3, TF4, TF5}
    Vx::TF1
    Vy::TF1
    rho::TF2  # different type to accomodate ReverseDiff
    pitch::TF3
    mu::TF4
    asound::TF5
end

# convenience constructor when Re and Mach are not used.
OperatingPoint(Vx, Vy, rho) = OperatingPoint(Vx, Vy, rho; pitch=zero(rho), mu=one(rho), asound=one(rho))

# convenience function to access fields within an array of structs
function Base.getproperty(obj::Vector{OperatingPoint{TF1, TF2, TF3, TF4, TF5}}, sym::Symbol) where {TF1, TF2, TF3, TF4, TF5}
    return getfield.(obj, sym)
end

"""
    Outputs(Np, Tp, a, ap, u, v, phi, alpha, W, cl, cd, cn, ct, F, G)

Outputs from the BEM solver along the radius.

**Arguments**
- `Np::Float64`: normal force per unit length
- `Tp::Float64`: tangential force per unit length
- `a::Float64`: axial induction factor
- `ap::Float64`: tangential induction factor
- `u::Float64`: axial induced velocity
- `v::Float64`: tangential induced velocity
- `phi::Float64`: inflow angle
- `alpha::Float64`: angle of attack
- `W::Float64`: inflow velocity
- `cl::Float64`: lift coefficient
- `cd::Float64`: drag coefficient
- `cn::Float64`: normal force coefficient
- `ct::Float64`: tangential force coefficient
- `F::Float64`: hub/tip loss correction
- `G::Float64`: effective hub/tip loss correction for induced velocities: `u = Vx * a * G, v = Vy * ap * G`
"""
struct Outputs{TF}
    Np::TF
    Tp::TF
    a::TF
    ap::TF
    u::TF
    v::TF
    phi::TF
    alpha::TF
    W::TF
    cl::TF
    cd::TF
    cn::TF
    ct::TF
    F::TF
    G::TF
end

# convenience constructor to initialize
Outputs() = Outputs(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

# convenience function to access fields within an array of structs
function Base.getproperty(obj::Vector{Outputs{TF}}, sym::Symbol) where TF
    return getfield.(obj, sym)
end

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

# ------------ BEM core ------------------

"""
(private) residual function
"""
function residual(phi, rotor, section, op)

    # unpack inputs
    r = section.r
    chord = section.chord
    theta = section.theta
    af = section.af

    Rhub = rotor.Rhub
    Rtip = rotor.Rtip
    B = rotor.B

    Vx = op.Vx
    Vy = op.Vy
    rho = op.rho
    pitch = op.pitch

    if real(chord) < zero(real(chord))
        @warn "chord = $chord m"
        chord = 0.1*0.0254*one(chord)
    end

    # constants
    sigma_p = B*chord/(2.0*pi*r)
    sphi = sin(phi)
    cphi = cos(phi)

    # angle of attack
    alpha = (theta + pitch) - phi

    # Reynolds/Mach number
    W0 = sqrt(Vx^2 + Vy^2)  # ignoring induction, which is generally a very minor difference and only affects Reynolds/Mach number
    Re = rho * W0 * chord / op.mu
    Mach = W0/op.asound  # also ignoring induction

    # airfoil cl/cd
    if rotor.turbine
        cl, cd = afeval(af, -alpha, Re, Mach)
        cl *= -1
    else
        cl, cd = afeval(af, alpha, Re, Mach)
    end

    # airfoil corrections
    if !isnothing(rotor.re)
        cl, cd = re_correction(rotor.re, cl, cd, Re)
    end
    if !isnothing(rotor.mach)
        cl, cd = mach_correction(rotor.mach, cl, cd, Mach)
    end
    if !isnothing(rotor.rotation)
        cl, cd = rotation_correction(rotor.rotation, cl, cd, chord/r, r/Rtip, Vy/Vx*Rtip/r, alpha, phi)
    end

    # resolve into normal and tangential forces
    cn = cl*cphi - cd*sphi
    ct = cl*sphi + cd*cphi

    # hub/tip loss
    F = one(chord)
    if !isnothing(rotor.tip)
        F = tip_correction(rotor.tip, r, Rhub, Rtip, phi, B)
    end

    # sec parameters
    k = cn*sigma_p/(4.0*F*sphi*sphi)
    kp = ct*sigma_p/(4.0*F*sphi*cphi)

    # --- solve for induced velocities ------
    if isapprox(real(Vx), 0.0, atol=1e-6)

        u = sign(phi)*kp*cn/ct*Vy
        v = zero(phi)
        a = zero(phi)
        ap = zero(phi)
        R = sign(phi) - k

    elseif isapprox(real(Vy), 0.0, atol=1e-6)

        u = zero(phi)
        v = k*ct/cn*abs(Vx)
        a = zero(phi)
        ap = zero(phi)
        R = sign(Vx) + kp

    else

        if real(phi) < 0
            k *= -1
        end

        if isapprox(real(k), 1.0, atol=1e-6)  # state corresopnds to Vx=0, return any nonzero residual
            return 1.0, Outputs()
        end

        if real(k) >= -2.0/3  # momentum region
            a = k/(1 - k)

        else  # empirical region
            g1 = F*(2*k - 1) + 10.0/9
            g2 = F*(F - 2*k - 4.0/3)
            g3 = 2*F*(1 - k) - 25.0/9

            if isapprox(real(g3), 0.0, atol=1e-6)  # avoid singularity
                a = 1.0/(2.0*sqrt(g2)) - 1
            else
                a = (g1 + sqrt(g2)) / g3
            end
        end

        u = a * Vx

        # -------- tangential induction ----------
        if real(Vx) < 0
            kp *= -1
        end

        if isapprox(real(kp), -1.0, atol=1e-6)  # state corresopnds to Vy=0, return any nonzero residual
            return 1.0, Outputs()
        end

        ap = kp/(1 + kp)
        v = ap * Vy

        # ------- residual function -------------
        R = sin(phi)/(1 + a) - Vx/Vy*cos(phi)/(1 - ap)
    end

    # ------- loads ---------
    W = sqrt((Vx + u)^2 + (Vy - v)^2)
    Np = cn*0.5*rho*W^2*chord
    Tp = ct*0.5*rho*W^2*chord

    # The BEM methodology applies hub/tip losses to the loads rather than to the velocities.
    # This is the most common way to implement a BEM, but it means that the raw velocities are misleading
    # as they do not contain any hub/tip loss corrections.
    # To fix this we compute the effective hub/tip losses that would produce the same thrust/torque.
    # In other words:
    # CT = 4 a (1 + a) F = 4 a G (1 + a G)\n
    # This is solved for G, then multiplied against the wake velocities.

    if isapprox(real(Vx), 0.0, atol=1e-6)
        G = sqrt(F)
    elseif isapprox(real(Vy), 0.0, atol=1e-6)
        G = F
    else
        G = (-1.0 + sqrt(1.0 + 4*a*(1.0 + a)*F))/(2*a)
    end
    u *= G
    v *= G

    if rotor.turbine
        return R, Outputs(-Np, -Tp, -a, -ap, -u, -v, phi, -alpha, W, -cl, cd, -cn, -ct, F, G)
    else
        return R, Outputs(Np, Tp, a, ap, u, v, phi, alpha, W, cl, cd, cn, ct, F, G)
    end

end

"""
(private) Find a bracket for the root closest to xmin by subdividing
interval (xmin, xmax) into n intervals.

Returns found, xl, xu.
If found = true a bracket was found between (xl, xu)
"""
function firstbracket(f, xmin, xmax, n, backwardsearch=false)

    xvec = range(xmin, xmax, length=n)
    if backwardsearch  # start from xmax and work backwards
        xvec = reverse(xvec)
    end

    fprev = f(xvec[1])
    for i = 2:n
        fnext = f(xvec[i])
        if fprev*fnext < 0  # bracket found
            if backwardsearch
                return true, xvec[i], xvec[i-1]
            else
                return true, xvec[i-1], xvec[i]
            end
        end
        fprev = fnext
    end

    return false, 0.0, 0.0
end

"""
    solve(rotor, section, op)

Solve the BEM equations for given rotor geometry and operating point.

**Arguments**
- `rotor::Rotor`: rotor properties
- `section::Section`: section properties
- `op::OperatingPoint`: operating point

**Returns**
- `outputs::Outputs`: BEM output data including loads, induction factors, etc.
"""
function solve(rotor, section, op)

    # error handling
    if typeof(section) <: Vector
        error("You passed in an vector for section, but this funciton does not accept an vector.\nProbably you intended to use broadcasting (notice the dot): solve.(Ref(rotor), sections, ops)")
    end

    # check if we are at hub/tip
    if isapprox(section.r, rotor.Rhub, atol=1e-6) || isapprox(section.r, rotor.Rtip, atol=1e-6)
        return Outputs()  # no loads at hub/tip
    end

    # parameters
    npts = 10  # number of discretization points to find bracket in residual solve

    # unpack
    Vx = op.Vx
    Vy = op.Vy
    theta = section.theta + op.pitch

    # ---- determine quadrants based on case -----
    Vx_is_zero = isapprox(Vx, 0.0, atol=1e-6)
    Vy_is_zero = isapprox(Vy, 0.0, atol=1e-6)

    # quadrants
    epsilon = 1e-6
    q1 = [epsilon, pi/2]
    q2 = [-pi/2, -epsilon]
    q3 = [pi/2, pi-epsilon]
    q4 = [-pi+epsilon, -pi/2]

    if Vx_is_zero && Vy_is_zero
        return Outputs()

    elseif Vx_is_zero

        startfrom90 = false  # start bracket at 0 deg.

        if Vy > 0 && theta > 0
            order = (q1, q2)
        elseif Vy > 0 && theta < 0
            order = (q2, q1)
        elseif Vy < 0 && theta > 0
            order = (q3, q4)
        else  # Vy < 0 && theta < 0
            order = (q4, q3)
        end

    elseif Vy_is_zero

        startfrom90 = true  # start bracket search from 90 deg

        if Vx > 0 && abs(theta) < pi/2
            order = (q1, q3)
        elseif Vx < 0 && abs(theta) < pi/2
            order = (q2, q4)
        elseif Vx > 0 && abs(theta) > pi/2
            order = (q3, q1)
        else  # Vx < 0 && abs(theta) > pi/2
            order = (q4, q2)
        end

    else  # normal case

        startfrom90 = false

        if Vx > 0 && Vy > 0
            order = (q1, q2, q3, q4)
        elseif Vx < 0 && Vy > 0
            order = (q2, q1, q4, q3)
        elseif Vx > 0 && Vy < 0
            order = (q3, q4, q1, q2)
        else  # Vx[i] < 0 && Vy[i] < 0
            order = (q4, q3, q2, q1)
        end

    end

    # ----- solve residual function ------

    # # wrapper to residual function to accomodate format required by fzero
    R(phi) = residual(phi, rotor, section, op)[1]

    success = false
    for j = 1:length(order)  # quadrant orders.  In most cases it should find root in first quadrant searched.
        phimin, phimax = order[j]

        # check to see if it would be faster to reverse the bracket search direction
        backwardsearch = false
        if !startfrom90
            if phimin == -pi/2 || phimax == -pi/2  # q2 or q4
                backwardsearch = true
            end
        else
            if phimax == pi/2  # q1
                backwardsearch = true
            end
        end

        # force to dual numbers if necessary
        phimin = phimin*one(section.chord)
        phimax = phimax*one(section.chord)

        # find bracket
        success, phiL, phiU = firstbracket(R, phimin, phimax, npts, backwardsearch)

        # once bracket is found, solve root finding problem and compute loads
        if success
            phistar, _ = FLOWMath.brent(R, phiL, phiU)
            _, outputs = residual(phistar, rotor, section, op)
            return outputs
        end
    end

    # it shouldn't get to this point.  if it does it means no solution was found
    # it will return empty outputs
    # alternatively, one could increase npts and try again

    @warn "Invalid data (likely) for this section.  Zero loading assumed."
    return Outputs()
end

# ------------ inflow ------------------

"""
    simple_op(Vinf, Omega, r, rho; pitch=0.0, mu=1.0, asound=1.0, precone=0.0)

Uniform inflow through rotor.  Returns an OperatingPoint object.

**Arguments**
- `Vinf::Float`: freestream speed (m/s)
- `Omega::Float`: rotation speed (rad/s)
- `r::Float`: radial location where inflow is computed (m)
- `rho::Float`: air density (kg/m^3)
- `pitch::Float`: pitch angle (rad)
- `mu::Float`: air viscosity (Pa * s)
- `asounnd::Float`: air speed of sound (m/s)
- `precone::Float`: precone angle (rad)
"""
function simple_op(Vinf, Omega, r, rho; pitch=zero(rho), mu=one(rho), asound=one(rho), precone=zero(Vinf))

    # error handling
    if typeof(r) <: Vector
        error("You passed in an vector for r, but this function does not accept an vector.\nProbably you intended to use broadcasting")
    end

    Vx = Vinf * cos(precone)
    Vy = Omega * r * cos(precone)

    return OperatingPoint(Vx, Vy, rho, pitch, mu, asound)

end

"""
    windturbine_op(Vhub, Omega, pitch, r, precone, yaw, tilt, azimuth, hubHt, shearExp, rho, mu=1.0, asound=1.0)

Compute relative wind velocity components along blade accounting for inflow conditions
and orientation of turbine.  See Documentation for angle definitions.

**Arguments**
- `Vhub::Float64`: freestream speed at hub (m/s)
- `Omega::Float64`: rotation speed (rad/s)
- `pitch::Float64`: pitch angle (rad)
- `r::Float64`: radial location where inflow is computed (m)
- `precone::Float64`: precone angle (rad)
- `yaw::Float64`: yaw angle (rad)
- `tilt::Float64`: tilt angle (rad)
- `azimuth::Float64`: azimuth angle to evaluate at (rad)
- `hubHt::Float64`: hub height (m) - used for shear
- `shearExp::Float64`: power law shear exponent
- `rho::Float64`: air density (kg/m^3)
- `mu::Float64`: air viscosity (Pa * s)
- `asound::Float64`: air speed of sound (m/s)
"""
function windturbine_op(Vhub, Omega, pitch, r, precone, yaw, tilt, azimuth, hubHt, shearExp, rho, mu=one(rho), asound=one(rho))
    sy = sin(yaw)
    cy = cos(yaw)
    st = sin(tilt)
    ct = cos(tilt)
    sa = sin(azimuth)
    ca = cos(azimuth)
    sc = sin(precone)
    cc = cos(precone)

    # coordinate in azimuthal coordinate system
    x_az = -r*sin(precone)
    z_az = r*cos(precone)
    y_az = 0.0  # could omit (the more general case allows for presweep so this is nonzero)

    # get section heights in wind-aligned coordinate system
    heightFromHub = (y_az*sa + z_az*ca)*ct - x_az*st

    # velocity with shear
    V = Vhub*(1 + heightFromHub/hubHt)^shearExp

    # transform wind to blade c.s.
    Vwind_x = V * ((cy*st*ca + sy*sa)*sc + cy*ct*cc)
    Vwind_y = V * (cy*st*sa - sy*ca)

    # wind from rotation to blade c.s.
    Vrot_x = -Omega*y_az*sc
    Vrot_y = Omega*z_az

    # total velocity
    Vx = Vwind_x + Vrot_x
    Vy = Vwind_y + Vrot_y

    # operating point
    return OperatingPoint(Vx, Vy, rho, pitch, mu, asound)

end

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

# -------- convenience methods ------------

"""
    thrusttorque(rotor, sections, outputs::Vector{TO}) where TO

integrate the thrust/torque across the blade,
including 0 loads at hub/tip, using a trapezoidal rule.

**Arguments**
- `rotor::Rotor`: rotor object
- `sections::Vector{Section}`: rotor object
- `outputs::Vector{Outputs}`: output data along blade

**Returns**
- `T::Float64`: thrust (along x-dir see Documentation).
- `Q::Float64`: torque (along x-dir see Documentation).
"""
function thrusttorque(rotor, sections, outputs::Vector{TO}) where TO

    # add hub/tip for complete integration.  loads go to zero at hub/tip.
    rvec = [s.r for s in sections]
    rfull = [rotor.Rhub; rvec; rotor.Rtip]
    Npfull = [0.0; outputs.Np; 0.0]
    Tpfull = [0.0; outputs.Tp; 0.0]

    # integrate Thrust and Torque (trapezoidal)
    thrust = Npfull*cos(rotor.precone)
    torque = Tpfull.*rfull*cos(rotor.precone)

    T = rotor.B * FLOWMath.trapz(rfull, thrust)
    Q = rotor.B * FLOWMath.trapz(rfull, torque)

    return T, Q
end

"""
    thrusttorque(rotor, sections, outputs::Matrix{TO}) where TO

Integrate the thrust/torque across the blade given an array of output data.
Generally used for azimuthal averaging of thrust/torque.
`outputs[i, j]` corresponds to `sections[i], azimuth[j]`.  Integrates across azimuth
"""
function thrusttorque(rotor, sections, outputs::Matrix{TO}) where TO

    T = 0.0
    Q = 0.0
    nr, naz = size(outputs)

    for j = 1:naz
        Tsub, Qsub = thrusttorque(rotor, sections, outputs[:, j])
        T += Tsub / naz
        Q += Qsub / naz
    end

    return T, Q
end

"""
    nondim(T, Q, Vhub, Omega, rho, rotor, rotortype)

Nondimensionalize the outputs.

**Arguments**
- `T::Float64`: thrust (N)
- `Q::Float64`: torque (N-m)
- `Vhub::Float64`: hub speed used in turbine normalization (m/s)
- `Omega::Float64`: rotation speed used in propeller normalization (rad/s)
- `rho::Float64`: air density (kg/m^3)
- `rotor::Rotor`: rotor object
- `rotortype::String`: normalization type

**Returns**

`if rotortype == "windturbine"`
- `CP::Float64`: power coefficient
- `CT::Float64`: thrust coefficient
- `CQ::Float64`: torque coefficient

`if rotortype == "propeller"`
- `eff::Float64`: efficiency
- `CT::Float64`: thrust coefficient
- `CQ::Float64`: torque coefficient

`if rotortype == "helicopter"`
- `FM::Float64`: figure of merit
- `CT::Float64`: thrust coefficient
- `CQ or CP::Float64`: torque/power coefficient (they are identical)
"""
function nondim(T, Q, Vhub, Omega, rho, rotor, rotortype)

    P = Q * Omega
    Rp = rotor.Rtip*cos(rotor.precone)

    if rotortype == "windturbine"  # wind turbine normalizations

        q = 0.5 * rho * Vhub^2
        A = pi * Rp^2

        CP = P / (q * A * Vhub)
        CT = T / (q * A)
        CQ = Q / (q * Rp * A)

        return CP, CT, CQ

    elseif rotortype == "propeller"

        n = Omega/(2*pi)
        Dp = 2*Rp

        if real(T) < 0
            eff = zero(typeof(T))  # creating drag not thrust
        else
            eff = T*Vhub/P
        end
        CT = T / (rho * n^2 * Dp^4)
        CQ = Q / (rho * n^2 * Dp^5)

        return eff, CT, CQ

    elseif rotortype == "helicopter"

        A = pi * Rp^2

        CT = T / (rho * A * (Omega*Rp)^2)
        CP = P / (rho * A * (Omega*Rp)^3)  # note that CQ = CP
        if real(T) < 0
            FM = zero(typeof(T))
        else
            FM = CT^(3.0/2)/(sqrt(2)*CP)
        end

        return FM, CT, CP
    end

end

end  # module


text:
#Test Values
test_values = (vec([852 845 846 851 846 846 840 825 823 821 836 854 854 850 822 805 802 804 804 823 815 799 802 793 792 807 821 823 809 809 823 840 858 847 823 799 788 788 789 799 805 805 817 853 873 896 900 886 862 840 845 854 858 861 830]))
# Ideal Case (No changes will be performed)
@testset "Ideal Case" begin
    fixed_signal, outliers, nonnormal = adaptiveHRVFilter(test_values)
    @test test_values == fixed_signal
    @test sum(outliers) == 0
    @test sum(nonnormal) == 0
end
fixed_signal, outliers, nonnormal = adaptiveHRVFilter(test_values)

# Changes but possibly physiological
changed_idxs = [20, 21]
changed_values = [630, 1340]
changed_test_values = copy(test_values)
changed_test_values[changed_idxs] = changed_values

@testset "Physiological Changes" begin
    fixed_signal, outliers, nonnormal = adaptiveHRVFilter(changed_test_values)
    notchanged = broadcast(!, nonnormal)
    @test sum(outliers) == 0
    @test changed_test_values[notchanged] == fixed_signal[notchanged]
end

# Changes to non physiological values
changed_idxs = [20, 21]
changed_values = [180, 2500]
changed_test_values = copy(test_values)
changed_test_values[changed_idxs] = changed_values

@testset "Non Physiological Changes" begin
    fixed_signal, outliers, nonnormal = adaptiveHRVFilter(changed_test_values)
    notchanged = broadcast(!, outliers)
    @test sum(nonnormal) == 0
    @test sum(outliers) == 2
    @test length(fixed_signal) + sum(outliers) == length(changed_test_values)
end


text:
#=
- curve http://jsxgraph.org/docs/symbols/Curve.html
- datacurve ? (via curve)
=#

mutable struct FunctionGraph <: Object
    f::JSFun
    a::Union{Nothing,<:Real}
    b::Union{Nothing,<:Real}
    # --- other board options --- jsxgraph.org/docs/symbols/Functiongraph.html
    opts::Option{LittleDict{Symbol,Any}}
end

"""
    functiongraph(f; a, b)

* `f` - a function of a single parameter `x` ex: (x -> 3x)
* `a` - optionally set the minimum value of `x`
* `b` - optionall set the maximum value of `x`
"""
function functiongraph(f::JSFun; a=nothing, b=nothing, kw...)
    return FunctionGraph(f, a, b, dict(;kw...))
end

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

function str(fg::FunctionGraph, b::Board)
    xmin = ifelse(isnothing(fg.a), b.xlim[1], fg.a)
    xmax = ifelse(isnothing(fg.b), b.xlim[2], fg.b)
    opts = get_opts(fg)
    fs, fss, frp = strf(fg.f, "FN", b)
    jss = js".create('functiongraph', [$fss, $xmin, $xmax], $opts);"
    return fs * b.name * replacefn(jss.s, frp)
end

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

abstract type Curve{X,Y,A,B} <: Object end

mutable struct ParametricCurve{X<:FR,Y<:FR,A<:FR,B<:FR} <: Curve{X,Y,A,B}
    x::X
    y::Y
    a::A
    b::B
    # -- opts -- jsxgraph.org/docs/symbols/Curve.html
    opts::Option{LittleDict{Symbol,Any}}
end
function parametriccurve(f::FR, g::FR; a::Option{FR}=nothing,
                         b::Option{FR}=nothing, kw...)
    isnothing(a) && (a = -10)
    isnothing(b) && (b = 10)
    return ParametricCurve(f, g, a, b, dict(;kw...))
end
pcurve = parametriccurve

mutable struct DataPlot{X<:AR,Y<:AFR} <: Curve{X,Y,Nothing,Nothing}
    x::X
    y::Y
    # -- opts -- jsxgraph.org/docs/symbols/Curve.html
    opts::Option{LittleDict{Symbol,Any}}
end
dataplot(x::AR, y::JSFun; kw...) = DataPlot(x, y, dict(;kw...))
dataplot(x::AR, y::AR; kw...) = (check_dims(x, y); DataPlot(x, y, dict(;kw...)))

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

function str(pc::ParametricCurve, b::Board)
    isnothing(pc.a) && (pc.a = -10)
    isnothing(pc.b) && (pc.b = xmax)
    xs, xss, xrp = strf(pc.x, "FX", b)
    ys, yss, yrp = strf(pc.y, "FY", b)
    as, ass, arp = strf(pc.a, "FA", b)
    bs, bss, brp = strf(pc.b, "FB", b)
    opts = get_opts(pc)
    jss = js".create('curve', [$xss, $yss, $ass, $bss], $opts);"
    return xs * ys * as * bs * b.name * replacefn(jss.s, xrp, yrp, arp, brp)
end

function str(dp::DataPlot, b::Board)
    opts = get_opts(dp)
    jss = js".create('curve', [$(dp.x), $(dp.y)], $opts);"
    return b.name * jss.s
end

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

plot(f::JSFun; kw...) = functiongraph(f; kw...)
plot(f::FR, g::FR; kw...) = parametriccurve(f, g; kw...)
plot(x::AR, y::AFR; kw...) = dataplot(x, y; kw...)

# # {curveType: 'polar'}
# mutable struct PolarCurve{R<:Function,O<:AFR,A<:FR,B<:FR} <: Curve{R,O,A,B}
#     r::R
#     offset::O
#     a::A
#     b::B
# end


text:
@testset "contractions" begin

D = 10
d = 3
sites = 5
T = Float32

ψ = CUDA.randn(CuMPS{T}, sites, D, d)
ϕ = CUDA.randn(CuMPS{T}, sites, D, d)

Id = [CuMatrix{T}(I(d)) for _ ∈ 1:length(ϕ)]

@testset "dot products" begin
    @test dot(ψ, ψ) ≈ dot(ψ, ψ)
    @test dot(ψ, ϕ) ≈ conj(dot(ϕ, ψ))

    @test dot(ψ, Id, ϕ) ≈ conj(dot(ϕ, Id, ψ))
    @test dot(ψ, Id, ϕ) ≈ dot(ψ, ϕ)

    @test norm(ψ) ≈ sqrt(abs(dot(ψ, ψ)))
end

@testset "left environment" begin
    L = left_env(ϕ, ψ)
    @test tr(L[end]) ≈ dot(ϕ, ψ)
end

@testset "Cauchy-Schwarz inequality" begin
    @test abs(dot(ϕ, ψ)) <= norm(ϕ) * norm(ψ)
end
end


text:
using PackageCompiler
create_sysimage([
:Pluto,
:PlutoUI,
:FFTW,
:IterativeSolvers,
:LinearAlgebra,
:MRIReco,
:PlotlyJS,
:Plots,
:Pluto,
:PlutoUI
];

precompile_execution_file="warmup.jl",
replace_default=true,
cpu_target = PackageCompiler.default_app_cpu_target())


text:
#####
##### Viscous fluxes
#####

@inline viscous_flux_ux(i, j, k, grid, clock, νᶜᶜᶜ::Number, u) = νᶜᶜᶜ * ℑxᶜᵃᵃ(i, j, k, grid, Axᵃᵃᶜ) * ∂xᶜᵃᵃ(i, j, k, grid, u)
@inline viscous_flux_uy(i, j, k, grid, clock, νᶠᶠᶜ::Number, u) = νᶠᶠᶜ * ℑyᵃᶠᵃ(i, j, k, grid, Ayᵃᵃᶜ) * ∂yᵃᶠᵃ(i, j, k, grid, u)
@inline viscous_flux_uz(i, j, k, grid, clock, νᶠᶜᶠ::Number, u) = νᶠᶜᶠ * ℑzᵃᵃᶠ(i, j, k, grid, Azᵃᵃᵃ) * ∂zᵃᵃᶠ(i, j, k, grid, u)

@inline viscous_flux_vx(i, j, k, grid, clock, νᶠᶠᶜ::Number, v) = νᶠᶠᶜ * ℑxᶠᵃᵃ(i, j, k, grid, Axᵃᵃᶜ) * ∂xᶠᵃᵃ(i, j, k, grid, v)
@inline viscous_flux_vy(i, j, k, grid, clock, νᶜᶜᶜ::Number, v) = νᶜᶜᶜ * ℑyᵃᶜᵃ(i, j, k, grid, Ayᵃᵃᶜ) * ∂yᵃᶜᵃ(i, j, k, grid, v)
@inline viscous_flux_vz(i, j, k, grid, clock, νᶜᶠᶠ::Number, v) = νᶜᶠᶠ * ℑzᵃᵃᶠ(i, j, k, grid, Azᵃᵃᵃ) * ∂zᵃᵃᶠ(i, j, k, grid, v)

@inline viscous_flux_wx(i, j, k, grid, clock, νᶠᶜᶠ::Number, w) = νᶠᶜᶠ * ℑxᶠᵃᵃ(i, j, k, grid, Axᵃᵃᶜ) * ∂xᶠᵃᵃ(i, j, k, grid, w)
@inline viscous_flux_wy(i, j, k, grid, clock, νᶜᶠᶠ::Number, w) = νᶜᶠᶠ * ℑyᵃᶠᵃ(i, j, k, grid, Ayᵃᵃᶜ) * ∂yᵃᶠᵃ(i, j, k, grid, w)
@inline viscous_flux_wz(i, j, k, grid, clock, νᶜᶜᶜ::Number, w) = νᶜᶜᶜ * ℑzᵃᵃᶜ(i, j, k, grid, Azᵃᵃᵃ) * ∂zᵃᵃᶜ(i, j, k, grid, w)

# Viscosities-as-functions

@inline viscous_flux_ux(i, j, k, grid, clock, ν::Function, u) =
        viscous_flux_ux(i, j, k, grid, clock, ν(xnode(Cell, i, grid), ynode(Cell, j, grid), znode(Cell, k, grid), clock.time), u)

@inline viscous_flux_uy(i, j, k, grid, clock, ν::Function, u) =
        viscous_flux_uy(i, j, k, grid, clock, ν(xnode(Cell, i, grid), ynode(Face, j, grid), znode(Face, k, grid), clock.time), u)

@inline viscous_flux_uz(i, j, k, grid, clock, ν::Function, u) =
        viscous_flux_uz(i, j, k, grid, clock, ν(xnode(Face, i, grid), ynode(Cell, j, grid), znode(Face, k, grid), clock.time), u)

@inline viscous_flux_vx(i, j, k, grid, clock, ν::Function, v) =
        viscous_flux_vx(i, j, k, grid, clock, ν(xnode(Face, i, grid), ynode(Face, j, grid), znode(Cell, k, grid), clock.time), v)

@inline viscous_flux_vy(i, j, k, grid, clock, ν::Function, v) =
        viscous_flux_vy(i, j, k, grid, clock, ν(xnode(Cell, i, grid), ynode(Cell, j, grid), znode(Cell, k, grid), clock.time), v)

@inline viscous_flux_vz(i, j, k, grid, clock, ν::Function, v) =
        viscous_flux_vz(i, j, k, grid, clock, ν(xnode(Cell, i, grid), ynode(Face, j, grid), znode(Face, k, grid), clock.time), v)

@inline viscous_flux_wx(i, j, k, grid, clock, ν::Function, w) =
        viscous_flux_wx(i, j, k, grid, clock, ν(xnode(Face, i, grid), ynode(Cell, j, grid), znode(Face, k, grid), clock.time), w)

@inline viscous_flux_wy(i, j, k, grid, clock, ν::Function, w) =
        viscous_flux_wy(i, j, k, grid, clock, ν(xnode(Cell, i, grid), ynode(Face, j, grid), znode(Face, k, grid), clock.time), w)

@inline viscous_flux_wz(i, j, k, grid, clock, ν::Function, w) =
        viscous_flux_wz(i, j, k, grid, clock, ν(xnode(Cell, i, grid), ynode(Cell, j, grid), znode(Cell, k, grid), clock.time), w)

#####
##### Viscous dissipation operators
#####

"""
    ∂ⱼνᵢⱼ∂ᵢu(i, j, k, grid, νˣ, νʸ, νᶻ, u)

Calculates viscous dissipation for the u-velocity via

    1/V * [δxᶠᵃᵃ(νˣ * ℑxᶜᵃᵃ(Ax) * δxᶜᵃᵃ(u)) + δyᵃᶜᵃ(νʸ * ℑyᵃᶠᵃ(Ay) * δyᵃᶠᵃ(u)) + δzᵃᵃᶜ(νᶻ * ℑzᵃᵃᶠ(Az) * δzᵃᵃᶠ(u))]

which will end up at the location `fcc`.
"""
@inline function ∂ⱼνᵢⱼ∂ᵢu(i, j, k, grid, clock, νˣ, νʸ, νᶻ, u)
    return 1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, viscous_flux_ux, clock, νˣ, u) +
                                    δyᵃᶜᵃ(i, j, k, grid, viscous_flux_uy, clock, νʸ, u) +
                                    δzᵃᵃᶜ(i, j, k, grid, viscous_flux_uz, clock, νᶻ, u))
end

"""
    ∂ⱼνᵢⱼ∂ᵢu(i, j, k, grid, νˣ, νʸ, νᶻ, v)

Calculates viscous dissipation for the v-velocity via

    1/Vᵛ * [δxᶜᵃᵃ(νˣ * ℑxᶠᵃᵃ(Ax) * δxᶠᵃᵃ(v)) + δyᵃᶠᵃ(νʸ * ℑyᵃᶜᵃ(Ay) * δyᵃᶜᵃ(v)) + δzᵃᵃᶜ(νᶻ * ℑzᵃᵃᶠ(Az) * δzᵃᵃᶠ(v))]

which will end up at the location `cfc`.
"""
@inline function ∂ⱼνᵢⱼ∂ᵢv(i, j, k, grid, clock, νˣ, νʸ, νᶻ, v)
    return 1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, viscous_flux_vx, clock, νˣ, v) +
                                    δyᵃᶠᵃ(i, j, k, grid, viscous_flux_vy, clock, νʸ, v) +
                                    δzᵃᵃᶜ(i, j, k, grid, viscous_flux_vz, clock, νᶻ, v))
end

"""
    ∂ⱼνᵢⱼ∂ᵢu(i, j, k, grid::Grid, νˣ, νʸ, νᶻ, w)

Calculates viscous dissipation for the w-velocity via

    1/Vʷ * [δxᶜᵃᵃ(νˣ * ℑxᶠᵃᵃ(Ax) * δxᶠᵃᵃ(w)) + δyᵃᶜᵃ(νʸ * ℑyᵃᶠᵃ(Ay) * δyᵃᶠᵃ(w)) + δzᵃᵃᶠ(νᶻ * ℑzᵃᵃᶜ(Az) * δzᵃᵃᶜ(w))]

which will end up at the location `ccf`.
"""
@inline function ∂ⱼνᵢⱼ∂ᵢw(i, j, k, grid, clock, νˣ, νʸ, νᶻ, w)
    return 1/Vᵃᵃᶠ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, viscous_flux_wx, clock, νˣ, w) +
                                    δyᵃᶜᵃ(i, j, k, grid, viscous_flux_wy, clock, νʸ, w) +
                                    δzᵃᵃᶠ(i, j, k, grid, viscous_flux_wz, clock, νᶻ, w))
end

#####
##### Viscous dissipation for isotropic viscosity
#####

@inline ∂ⱼνᵢⱼ∂ᵢu(i, j, k, grid, clock, ν, u) = ∂ⱼνᵢⱼ∂ᵢu(i, j, k, grid, clock, ν, ν, ν, u)
@inline ∂ⱼνᵢⱼ∂ᵢv(i, j, k, grid, clock, ν, v) = ∂ⱼνᵢⱼ∂ᵢv(i, j, k, grid, clock, ν, ν, ν, v)
@inline ∂ⱼνᵢⱼ∂ᵢw(i, j, k, grid, clock, ν, w) = ∂ⱼνᵢⱼ∂ᵢw(i, j, k, grid, clock, ν, ν, ν, w)


text:
module ParameterSchedulers

using Flux

export Lambda, Sequence, Loop,
       Step, Exp, Poly, Inv,
       Tri, TriDecay2, TriExp,
       Sin, SinDecay2, SinExp,
       Cos,
       ScheduleIterator, ScheduledOptim, next!

include("schedule.jl")
include("decay.jl")
include("cyclic.jl")
include("optimizers.jl")

end


text:
module DummyPool

# dummy allocator that passes through any requests, calling into the GC if that fails.

using ..CuArrays: @pool_timeit, actual_alloc, actual_free

using CUDAdrv

init() = return

deinit() = @assert isempty(allocated) "Cannot deinitialize memory pool with outstanding allocations"

const allocated = Dict{CuPtr{Nothing},Int}()

function alloc(sz)
    ptr = nothing
    for phase in 1:3
        if phase == 2
            @pool_timeit "$phase.0 gc (incremental)" GC.gc(VERSION >= v"1.4.0-DEV.257" ? GC.Incremental : false)
        elseif phase == 3
            @pool_timeit "$phase.0 gc (full)" GC.gc(VERSION >= v"1.4.0-DEV.257" ? GC.Full : true)
        end

        @pool_timeit "$phase.1 alloc" begin
            ptr = actual_alloc(sz)
        end
        ptr === nothing || break
    end

    if ptr !== nothing
        allocated[ptr] = sz
        return ptr
    else
        return nothing
    end
end

function free(ptr)
    sz = allocated[ptr]
    delete!(allocated, ptr)
    actual_free(ptr)
    return
end

used_memory() = isempty(allocated) ? 0 : sum(sizeof, values(allocated))

cached_memory() = 0

dump() = return

end


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb03ld(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    NMAX = 50
    LDA = NMAX÷2
    LDB = NMAX÷2
    LDDE = NMAX÷2
    LDFG = NMAX÷2
    LDQ = 2*NMAX
    LDWORK = 8*NMAX*NMAX + max( 8*NMAX + 32, NMAX÷2 + 168, 272 )
    LIWORK = max( 32, NMAX + 12, NMAX*2 + 3 )
    A = Array{Float64,2}(undef, LDA,NMAX÷2)
    DE = Array{Float64,2}(undef, LDDE,NMAX÷2+1)
    B = Array{Float64,2}(undef, LDB,NMAX÷2)
    FG = Array{Float64,2}(undef, LDFG,NMAX÷2+1)
    Q = Array{Float64,2}(undef, LDQ,2*NMAX)
    ALPHAR = Array{Float64,1}(undef, NMAX÷2)
    ALPHAI = Array{Float64,1}(undef, NMAX÷2)
    BETA = Array{Float64,1}(undef, NMAX÷2)
    BWORK = Array{BlasBool,1}(undef, NMAX÷2)
    IWORK = Array{BlasInt,1}(undef, LIWORK)
    DWORK = Array{Float64,1}(undef, LDWORK)
f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    COMPQ =  vs[1][1]
    ORTH =  vs[2][1]
    N = parse(BlasInt, vs[3])
    if ( N<0 || N>NMAX )
        @error "Illegal N=$N"
    end

    M = N÷2

    vs = String[]
    _isz,_jsz = (M,M)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (M,M+1)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       DE[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (M,M)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       B[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (M,M+1)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       FG[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
    close(f)
# interp call 1

    NEIG, INFO = SLICOT.mb03ld!(COMPQ, ORTH, N, A, DE, B, FG, Q, ALPHAR, ALPHAI, BETA, LIWORK)
    @test INFO == 0
    INFO == 0 || return
    println(io, "NEIG = $NEIG")

# interp output 1
    println(io, "A:")
    _nc = M
    _nr = M
    show(io, "text/plain", A[1:_nr,1:_nc])
    println(io)

# interp output 2
    println(io, "DE:")
    _nc = M+1
    _nr = M
    show(io, "text/plain", DE[1:_nr,1:_nc])
    println(io)

# interp output 3
    println(io, "B:")
    _nc = M
    _nr = M
    show(io, "text/plain", B[1:_nr,1:_nc])
    println(io)

# interp output 4
    println(io, "FG:")
    show(io, "text/plain", FG[1:M,2:M+1])
    println(io)

# interp output 5
    println(io, "ALPHAR:")
    _nr = M
    show(io, "text/plain", ALPHAR[1:_nr])
    println(io)

# interp output 6
    println(io, "ALPHAI:")
    _nr = M
    show(io, "text/plain", ALPHAI[1:_nr])
    println(io)

# interp output 7
    println(io, "BETA:")
    _nr = M
    show(io, "text/plain", BETA[1:_nr])
    println(io)

if ( LSAME( COMPQ, 'C' ) && NEIG>0 )
# interp output 8
    println(io, "Q:")
    _nc = NEIG
    _nr = N
    show(io, "text/plain", Q[1:_nr,1:_nc])
    println(io)

end # if

end # run_mb03ld()


text:
struct DiracDelta{D} <: GridapType
  Γ::Triangulation{D}
  dΓ::Measure
end

function DiracDelta{D}(
  model::DiscreteModel,
  face_to_bgface::AbstractVector{<:Integer},
  degree::Integer) where D

  @assert 0 <= D && D < num_cell_dims(model) """\n
  Incorrect value of D=$D for building a DiracDelta{D} on a model with $(num_cell_dims(model)) cell dims.

  D should be in [0,$(num_cell_dims(model))).
  """

  topo = get_grid_topology(model)
  bgface_grid = Grid(ReferenceFE{D},model)
  face_grid = view(bgface_grid,face_to_bgface)
  cell_grid = get_grid(model)
  bgface_to_lcell = Fill(1,num_faces(model,D))
  glue = Geometry.FaceToCellGlue(topo,cell_grid,face_grid,face_to_bgface,bgface_to_lcell)
  trian = BodyFittedTriangulation(model,face_grid,face_to_bgface)
  Γ = BoundaryTriangulation(trian,glue)
  dΓ = Measure(Γ,degree)
  DiracDelta(Γ,dΓ)
end

function DiracDelta{D}(
  model::DiscreteModel,
  bgface_to_mask::AbstractVector{Bool},
  degree::Integer) where D

  face_to_bgface = findall(bgface_to_mask)
  DiracDelta{D}(model,face_to_bgface,degree)
end

function DiracDelta{D}(
  model::DiscreteModel,
  labeling::FaceLabeling,
  degree::Integer;tags) where D

  @assert 0 <= D && D < num_cell_dims(model) """\n
  Incorrect value of D=$D for building a DiracDelta{D} on a model with $(num_cell_dims(model)) cell dims.

  D should be in [0,$(num_cell_dims(model))).
  """
  face_to_mask = get_face_mask(labeling,tags,D)
  DiracDelta{D}(model,face_to_mask,degree)
end

function DiracDelta{D}(model::DiscreteModel,degree::Integer;tags) where D
  labeling = get_face_labeling(model)
  DiracDelta{D}(model,labeling,degree,tags=tags)
end

function DiracDelta{0}(model::DiscreteModel,labeling::FaceLabeling;tags)
  degree = 0
  DiracDelta{0}(model,labeling,degree;tags=tags)
end

function DiracDelta{0}(model::DiscreteModel;tags)
  degree = 0
  DiracDelta{0}(model,degree;tags=tags)
end

function (d::DiracDelta)(f)
  evaluate(d,f)
end

function evaluate!(cache,d::DiracDelta,f)
  ∫(f)*d.dΓ
end

function get_triangulation(d::DiracDelta)
  d.Γ
end


text:
# # Example: long, single end reads
#
# This is an example generated from this source
# file: [`se-example.jl`](@__REPO_ROOT_URL__/examples/sequencing/se-example.jl)
# You are seeing the
#md # online documentation version. The corresponding notebook can be found
#md # here: [`se-example.ipynb`](@__NBVIEWER_ROOT_URL__/man/sequencing/examples/se-example.ipynb)
#nb # jupyter notebook version. The corresponding online documentation page can
#nb # be found here: [`se-example.html`](https://bioinfologics.github.io/Pseudoseq.jl/latest/examples/sequencing/se-example),
# and the script can be found here: [`se-example.jl`](./se-example.jl)

# Let's see how you might simulate something like an Oxford Nanopore sequencing
# experiment.
#
# For the simulation we are going to:
#
# 1. Create a pool of 5000 copies of a reference genome.
# 2. Fragment the DNA molecules in the pool, to an average length of 40,000bp.
# 3. Subsample the molecules in the pool to achieve approximatly 30x coverage.
# 4. Create a set of single-end reads, the enitre length of each molecule.
# 5. Apply errors to the reads at a rate of 0.10 (1 error every 10bp).
# 6. Generate an output FASTQ file.

using Pseudoseq.Sequencing

# Starting with a FASTA formatted file containing the genome we want to sequence,
# we create a pool with a copy of the genome.

m = Molecules("ecoli-ref.fasta")

# We want to first have 5000 full copies of the genome, so we will make an amplifier.

amp = amplify(5000)

# Next we want a fragmenter to make a pool of shorter DNA molecules, with an average
# length of 40,000bp.

frag = fragment(40000bp)

# Next we create a subsampler which will randomly sample molecules, to give us
# a desired expected coverage. We specify 40000bp as our predicted average read
# length.

ssmpl = subsample(50X, 40000bp)

# We then want a read-maker that will give us single end reads.
# Since we call `makereads` without providing a read length, the function will
# generate reads from the entire length of each molecule in the pool. We do this
# to emulate what Nanopore sequencing is supposed to do: It takes an entire DNA
# fragment, feeds it through an electrically charged pore, producing a read for
# the entire fragment.

readmaker = makereads()

# Once we have reads, we will mark positions in the reads that are incorrectly
# detected by the sequencer: errors.
# We will construct a `FixedProbSubstitutions` function with a per base error
# probability of 0.1 and pass it to the `make_substitutions` method.
# This will make errors fall totally randomly over each read.

errmaker = make_substitutions(FixedProbSubstitutions(0.1))

# Now we can push our molecules through a pipeline of these processors, and out
# to a FASTQ file:

m |> amp |> frag |> ssmpl |> readmaker |> errmaker |> generate("se-reads.fastq")

# You can also compose the processors together into one whole function.
# Typing \circ in the julia repl and then hitting tab gives you the circular
# composition symbol. Note how pipelining above progresses from left to right,
# but composition is right to left in order.

my_protocol = errmaker ∘ readmaker ∘ ssmpl ∘ frag ∘ amp

m |> my_protocol |> generate("se-reads.fastq")


text:
module Spark

export
    SparkConf,
    SparkContext,
    RDD,
    JuliaRDD,
    JavaRDD,
    text_file,
    parallelize,
    map,
    map_pair,
    map_partitions,
    map_partitions_pair,
    map_partitions_with_index,
    reduce,
    filter,
    collect,
    count,
    id,
    num_partitions,
    close,
    @attach,
    share_variable,
    @share,
    flat_map,
    flat_map_pair,
    cartesian,
    group_by_key,
    reduce_by_key,
    cache,
    repartition,
    coalesce,
    pipe,
    # SQL
    SparkSession,
    Dataset,
    sql,
    read_json,
    write_json,
    read_parquet,
    write_parquet,
    read_df,
    write_df

include("core.jl")

end


text:
using Distributed;

nheads = @distributed (+) for i = 1:200000000
    Int(rand(Bool))
end

println("nheads = $nheads")


text:
module FoldsChainRulesBenchmarks

using BenchmarkTools: Benchmark, BenchmarkGroup

include("bench_seq_simple_sum.jl")
include("bench_seq_squared_sum.jl")

function setup()
    suite = BenchmarkGroup()
    suite["SeqSimpleSum"] = BenchSeqSimpleSum.setup()
    suite["SeqSquaredSum"] = BenchSeqSquaredSum.setup()
    return suite
end

function set_smoke_params!(bench)
    bench.params.seconds = 0.001
    bench.params.evals = 1
    bench.params.samples = 1
    bench.params.gctrial = false
    bench.params.gcsample = false
    return bench
end

foreach_benchmark(f!, bench::Benchmark) = f!(bench)
function foreach_benchmark(f!, group::BenchmarkGroup)
    for x in values(group)
        foreach_benchmark(f!, x)
    end
end

function setup_smoke()
    suite = setup()
    foreach_benchmark(set_smoke_params!, suite)
    return suite
end

function clear()
    BenchSeqSimpleSum.clear()
    BenchSeqSquaredSum.clear()
end

end  # module FoldsChainRulesBenchmarks


text:
include("testSAforDivSigGrad.jl")
include("testGMGRAPforPoisson.jl")
include("testGMG.jl")
include("testLinSolveMGWrapper.jl");


text:
#### Functions for quasi-canard finder
#### "Slow organisms exhibit sudden population disappearances in a reddened world" by Greyson-Gaito, Gellner, & McCann.

function orientation(p1, p2, p3)
    val = (p2[2] - p1[2]) * (p3[1] - p2[1]) - (p3[2] - p2[2]) * (p2[1] - p1[1])
    if (val == 0)
        return 0 #colinear
    elseif (val > 0)
        return 1 #clockwise
    else
        return -1 #anticlockwise
    end
end # Part 1 function from https://www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/

function dointersect(p1, p2, q1, q2)
    o1 = orientation(p1, p2, q1)
    o2 = orientation(p1, p2, q2)
    o3 = orientation(q1, q2, p1)
    o4 = orientation(q1, q2, p2)

    if (o1 != o2 && o3 != o4)
        return true
    else
        return false
    end
end # Part 2 function from https://www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/

function cf_returnmap_check(sol, pass_points, rm_point1, rm_point2, which_pass)
    if which_pass ∈ ["first", "second"] == false
        error("which_pass must be either first or second")
    end
    rm_pass_points = []
    for j in 1:length(pass_points)
        for l in Int64(pass_points[j][1]):length(sol)-1
            if sol.u[l+1][1] < rm_point1[1] && sol.u[l][1] > rm_point1[1]
                if dointersect(sol.u[l],sol.u[l+1],rm_point1,rm_point2) == true #Stricter condition for intersection of return map
                    if which_pass == "first"
                        append!(rm_pass_points, [[l , sol.u[l][1], sol.u[l][2]]])
                    else
                        append!(rm_pass_points, [[l , sol.u[l][1], sol.u[l][2]]])
                        break
                    end
                end
            end
        end
    end
    if length(rm_pass_points) < 1
        return false
    else
        return [true, rm_pass_points]
    end
end #First (and last) check in canard finder (Steps 1 and 5 in SI). Whether trajectory intersects small line at the maximum of the resource isocline.

function cf_resaxis_check(sol, pass_points, res_Hopf_point, res_lims, con_lims)
    resaxis_pass_points = []
    for j in 1:length(pass_points)
        for l in Int64(pass_points[j][1]):length(sol)-1
            if sol.u[l+1][1] < res_Hopf_point
                if res_lims[1] < sol.u[l][1] < res_lims[2] && con_lims[1] < sol.u[l][2] < con_lims[2]
                    append!(resaxis_pass_points, [[l , sol.u[l][1], sol.u[l][2]]])
                    break
                end
            else
                break
            end
        end
    end
    if length(resaxis_pass_points) < 1
        return false
    else
        return [true, resaxis_pass_points]
    end
end #Second check in canard finder (Step 2 in SI). Whether trajectory enters thin box on the consumer axis from just above the intersection of the resource isocline and the consumer axis up to the maximum of the resource isocline

function cf_box_check(sol, pass_points, res_lims, con_lims)
    new_pass_points = []
    for j in 1:length(pass_points)
        for l in Int64(pass_points[j][1]):length(sol)-1
            if res_lims[1] < sol.u[l][1] < res_lims[2] && con_lims[1] < sol.u[l][2] < con_lims[2]
                append!(new_pass_points, [[l , sol.u[l][1], sol.u[l][2]]])
                break
            end
        end
    end
    if length(new_pass_points) < 1
        return false
    else
        return [true, new_pass_points]
    end
end #Third check in canard finder (Step 3 in SI). Whether trajectory enters thin box on the consumer axis from 0.0 up to specified point on consumer axis.

function res_iso(model, R, p)
    if model == "RozMac"
        @unpack a, k, r, h = p
        return r * (a * h * k * R - a * h * R^2 + k - R) / (a * k)
    else
        @unpack x, y, δ, R₀, k  = p
        return (R*(R*δ - R + R₀*δ - R₀) + k*(-R*δ + R - R₀*δ + R₀))/(k*x*y)
    end
end #function to help calculate the resource isocline for both the Rosenzweig-MacArthur and Yodzis-Innes models

function cf_ressiocline_check(model, sol, pass_points, res_Hopf_point, par) #TODO maybe do it so checks whether stays within bubble for certain period of time?
    new_pass_points = []
    for j in 1:length(pass_points)
        for l in Int64(pass_points[j][1]):length(sol)-5
            if sol.u[l][1] > res_Hopf_point && isapprox(sol.u[l][2], res_iso(model, sol.u[l][1], par); atol = 1e-1)
                append!(new_pass_points, [[l , sol.u[l][1], sol.u[l][2]]])
                break
            end
        end
    end

    if length(new_pass_points) < 1
        return false
    else
        return [true, new_pass_points]
    end
end #Fourth check in canard finder (Step 4 in SI). Whether trajectory crosses the resource isocline.

function axial_checker(sol)
    if isapprox(sol.u[end][1], 3.0; atol = 1e-2)  && isapprox(sol.u[end][2], 0.00; atol = 1e-2)
            return "axial"
        else
            return "nothing"
    end
end #Function to check whether trajectory lands on the axial solution or not

function hopfpointsfinder(model, effR0)
    if model == "RozMac"
        @unpack k, a, h, r, m = RozMacPar(e = effR0)
        hopf_R = k / 2 - 1 / (2 * a * h)
        hopf_C = r * (a * h * k * hopf_R - a * h * hopf_R^2 + k - hopf_R) / (a * k)
        zero_C = r * k  / (a * k)
    elseif model == "YodInn"
        @unpack k, R₀, δ, x, y = YodInnScalePar(R₀ = effR0)
        hopf_R = -R₀/2 + k/2
        hopf_C = (hopf_R*(hopf_R*δ - hopf_R + R₀*δ - R₀) + k*(-hopf_R*δ + hopf_R - R₀*δ + R₀))/(k*x*y)
        zero_C = (k*(- R₀*δ + R₀))/(k*x*y)
    else
        error("model must be either RozMac or YodInn")
    end
    return vcat(hopf_R, hopf_C, zero_C)
end #function to find the resource and consumer values of the maximum point on the resource isocline. Also returns the consumer value of the intersection of the resource isocline with the consumer axis.

function cf_returnmap(model, ep, effR0, freq, r, seed, tsend, tvals)
    if model == "RozMac"
        sol = RozMac_pert(ep, effR0, freq, r, seed, tsend, tvals)
        par = RozMacPar(e = effR0, ε = ep)
    elseif model == "YodInn"
        sol = YodInn_pert(ep, effR0, freq, r, seed, tsend, tvals)
        par = YodInnScalePar(R₀ = effR0, ε = ep)
    else
        error("model must be either RozMac or YodInn")
    end
    hopfpoints = hopfpointsfinder(model, effR0)
    res_Hopf_point = hopfpoints[1]
    rm_point1 = [hopfpoints[1], hopfpoints[2]-(hopfpoints[2]*0.1)]
    rm_point2 = [hopfpoints[1], hopfpoints[2]+(hopfpoints[2]*0.1)]
    resaxisboxcheckpoints = [0.0, hopfpoints[3]*0.8]
    rm_pass_points1 = cf_returnmap_check(sol, [[1 , sol.u[1][1], sol.u[1][2]]], rm_point1, rm_point2, "first")
    if rm_pass_points1 == false
        return axial_checker(sol)
    else
        resaxis_pass_points = cf_resaxis_check(sol, rm_pass_points1[2], res_Hopf_point, [0.0, 0.1], [hopfpoints[3], hopfpoints[2]])
    end
    if resaxis_pass_points == false
        return axial_checker(sol)
    else
        resaxisbox_pass_points = cf_box_check(sol, resaxis_pass_points[2], [0.0,0.1], resaxisboxcheckpoints)
    end
    if resaxisbox_pass_points == false
        return axial_checker(sol)
    else
        resiso_pass_points = cf_ressiocline_check(model, sol, resaxisbox_pass_points[2], res_Hopf_point, par)
    end
    if resiso_pass_points == false
        return axial_checker(sol)
    else
        rm_pass_points2 = cf_returnmap_check(sol, resiso_pass_points[2],  rm_point1, rm_point2, "second")
    end
    if rm_pass_points2 == false
        return axial_checker(sol)
    else
        return "canard"
    end
end #Function to run through each checker and returns either "canard", "axial", or "nothing" for a single simulation of the Rosenzweig-MacArthur or Yodzis-Innes models.


text:
using Base
using GraffSDK

# 1. Get a Synchrony configuration
# Assume that you're running in local directory
cd(joinpath(Pkg.dir("GraffSDK"),"examples"))
config = loadGraffConfig("synchronyConfig.json")
#Update the config (updates internal config too as it's by reference)
sessionId = "HexagonalDrive13"
config.sessionId = sessionId
println(config)

# 2a. Create session
# Creating a session allows us to ingest data and start refining the factor graph.
# We can also mine data across all our robots and our sessions.
# Using the defaults, we are going to create a session that is of type Pose2 (2D poses) and it auto initialized (it makes first pose and attaches a prior)
newSession = SessionDetailsRequest(sessionId, "A test dataset demonstrating data ingestion for a wheeled vehicle driving in a hexagon.")
if !isSessionExisting(sessionId)
    retSession = addSession(newSession)
    @show retSession
end
# Now we can get it as well if we want
retSession = getSession()

# 2b. Get all sessions
# Just for example's sake - let's retrieve all sessions associated with out robot
sessions = getSessions()
@show sessions

# Get all nodes for our session
nodes = getNodes()
@show nodes
# Get a specific node
nodeDetails = getNode(nodes.nodes[1].id)


text:
"""
    @argcheck cond text

Throw an `ArgumentError` if `cond` is `false`.
"""
macro argcheck(ex, msgs...)
    msg = isempty(msgs) ? ex : msgs[1]
    if isa(msg, AbstractString)
        msg = msg # pass-through
    elseif !isempty(msgs) && (isa(msg, Expr) || isa(msg, Symbol))
        # message is an expression needing evaluating
        msg = :(Main.Base.string($(esc(msg))))
    elseif isdefined(Main, :Base) && isdefined(Main.Base, :string) && applicable(Main.Base.string, msg)
        msg = Main.Base.string(msg)
    else
        # string() might not be defined during bootstrap
        msg = quote
            msg = $(Expr(:quote,msg))
            isdefined(Main, :Base) ? Main.Base.string(msg) :
                (Core.println(msg); "Error during bootstrap. See stdout.")
        end
    end
    return :($(esc(ex)) ? $(nothing) : throw(ArgumentError($msg)))
end