text:
"""
    frank_wolfe(f, grad!, lmo, x0; ...)

Simplest form of the Frank-Wolfe algorithm.
Returns a tuple `(x, v, primal, dual_gap, traj_data)` with:
- `x` final iterate
- `v` last vertex from the LMO
- `primal` primal value `f(x)`
- `dual_gap` final Frank-Wolfe gap
- `traj_data` vector of trajectory information.
"""
function frank_wolfe(
    f,
    grad!,
    lmo,
    x0;
    line_search::LineSearchMethod=Adaptive(),
    L=Inf,
    gamma0=0,
    step_lim=20,
    momentum=nothing,
    epsilon=1e-7,
    max_iteration=10000,
    print_iter=1000,
    trajectory=false,
    verbose=false,
    linesearch_tol=1e-7,
    emphasis::Emphasis=memory,
    gradient=nothing,
    callback=nothing,
    timeout=Inf,
    print_callback=print_callback,
)

    # format string for output of the algorithm
    format_string = "%6s %13s %14e %14e %14e %14e %14e\n"
    t = 0
    dual_gap = Inf
    primal = Inf
    v = []
    x = x0
    tt = regular
    traj_data = []
    if trajectory && callback === nothing
        callback = trajectory_callback(traj_data)
    end
    time_start = time_ns()

    if line_search isa Shortstep && !isfinite(L)
        println("FATAL: Lipschitz constant not set. Prepare to blow up spectacularly.")
    end

    if line_search isa FixedStep && gamma0 == 0
        println("FATAL: gamma0 not set. We are not going to move a single bit.")
    end

    if (!isnothing(momentum) && line_search isa Union{Shortstep,Adaptive,RationalShortstep})
        println(
            "WARNING: Momentum-averaged gradients should usually be used with agnostic stepsize rules.",
        )
    end

    if verbose
        println("\nVanilla Frank-Wolfe Algorithm.")
        NumType = eltype(x0)
        println(
            "EMPHASIS: $emphasis STEPSIZE: $line_search EPSILON: $epsilon MAXITERATION: $max_iteration TYPE: $NumType",
        )
        grad_type = typeof(gradient)
        println("MOMENTUM: $momentum GRADIENTTYPE: $grad_type")
        if emphasis === memory
            println("WARNING: In memory emphasis mode iterates are written back into x0!")
        end
        headers = ["Type", "Iteration", "Primal", "Dual", "Dual Gap", "Time", "It/sec"]
        print_callback(headers, format_string, print_header=true)
    end
    if emphasis == memory && !isa(x, Union{Array,SparseArrays.AbstractSparseArray})
        # if integer, convert element type to most appropriate float
        if eltype(x) <: Integer
            x = copyto!(similar(x, float(eltype(x))), x)
        else
            x = copyto!(similar(x), x)
        end
    end
    first_iter = true
    # instanciating container for gradient
    if gradient === nothing
        gradient = similar(x)
    end

    # container for direction
    d = similar(x)
    gtemp = if momentum === nothing
        nothing
    else
        similar(x)
    end

    while t <= max_iteration && dual_gap >= max(epsilon, eps())

        #####################
        # managing time and Ctrl-C
        #####################
        time_at_loop = time_ns()
        if t == 0
            time_start = time_at_loop
        end
        # time is measured at beginning of loop for consistency throughout all algorithms
        tot_time = (time_at_loop - time_start) / 1e9

        if timeout < Inf
            if tot_time ≥ timeout
                if verbose
                    @info "Time limit reached"
                end
                break
            end
        end

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

        if momentum === nothing || first_iter
            grad!(gradient, x)
            if momentum !== nothing
                gtemp .= gradient
            end
        else
            grad!(gtemp, x)
            @emphasis(emphasis, gradient = (momentum * gradient) + (1 - momentum) * gtemp)
        end
        first_iter = false
        v = if is_simple_lmo(lmo)
            compute_extreme_point(lmo, gradient)
        else
            compute_extreme_point(lmo, gradient; call_counter=call_counter)
        end
        # go easy on the memory - only compute if really needed
        if (
            (mod(t, print_iter) == 0 && verbose) ||
            callback !== nothing ||
            line_search isa Shortstep
        )
            primal = f(x)
            dual_gap = fast_dot(x, gradient) - fast_dot(v, gradient)
        end

        @emphasis(emphasis, d = x - v)

        gamma, L = line_search_wrapper(
            line_search,
            t,
            f,
            grad!,
            x,
            d,
            momentum === nothing ? gradient : gtemp, # use appropriate storage
            dual_gap,
            L,
            gamma0,
            linesearch_tol,
            step_lim,
            one(eltype(x)),
        )
        if callback !== nothing
            state = (
                t=t,
                primal=primal,
                dual=primal - dual_gap,
                dual_gap=dual_gap,
                time=tot_time,
                x=x,
                v=v,
                gamma=gamma,
                f=f,
                grad!=grad!,
                lmo=lmo,
            )
            callback(state)
        end

        @emphasis(emphasis, x = x - gamma * d)

        if (mod(t, print_iter) == 0 && verbose)
            tt = regular
            if t == 0
                tt = initial
            end

            rep = (
                st[Symbol(tt)],
                string(t),
                Float64(primal),
                Float64(primal - dual_gap),
                Float64(dual_gap),
                tot_time,
                t / tot_time,
            )
            print_callback(rep, format_string)

            flush(stdout)
        end
        t = t + 1
    end
    # recompute everything once for final verfication / do not record to trajectory though for now!
    # this is important as some variants do not recompute f(x) and the dual_gap regularly but only when reporting
    # hence the final computation.

    grad!(gradient, x)
    v = compute_extreme_point(lmo, gradient)
    primal = f(x)
    dual_gap = fast_dot(x, gradient) - fast_dot(v, gradient)
    if verbose
        tt = last
        tot_time = (time_ns() - time_start) / 1.0e9
        rep = (
            st[Symbol(tt)],
            string(t - 1),
            Float64(primal),
            Float64(primal - dual_gap),
            Float64(dual_gap),
            tot_time,
            t / tot_time,
        )
        print_callback(rep, format_string)
        print_callback(nothing, format_string, print_footer=true)
        flush(stdout)
    end
    return x, v, primal, dual_gap, traj_data
end

"""
    lazified_conditional_gradient

Similar to [`frank_wolfe`](@ref) but lazyfying the LMO:
each call is stored in a cache, which is looked up first for a good-enough direction.
The cache used is a [`FrankWolfe.MultiCacheLMO`](@ref) or a [`FrankWolfe.VectorCacheLMO`](@ref)
depending on whether the provided `cache_size` option is finite.
"""
function lazified_conditional_gradient(
    f,
    grad!,
    lmo_base,
    x0;
    line_search::LineSearchMethod=Adaptive(),
    L=Inf,
    gamma0=0,
    lazy_tolerance=2.0,
    cache_size=Inf,
    greedy_lazy=false,
    epsilon=1e-7,
    max_iteration=10000,
    print_iter=1000,
    trajectory=false,
    verbose=false,
    linesearch_tol=1e-7,
    step_lim=20,
    emphasis::Emphasis=memory,
    gradient=nothing,
    VType=typeof(x0),
    callback=nothing,
    timeout=Inf,
    print_callback=print_callback,
)

    # format string for output of the algorithm
    format_string = "%6s %13s %14e %14e %14e %14e %14e %14i\n"

    if isfinite(cache_size)
        lmo = MultiCacheLMO{cache_size,typeof(lmo_base),VType}(lmo_base)
    else
        lmo = VectorCacheLMO{typeof(lmo_base),VType}(lmo_base)
    end

    t = 0
    dual_gap = Inf
    primal = Inf
    v = []
    x = x0
    phi = Inf
    traj_data = []
    if trajectory && callback === nothing
        callback = trajectory_callback(traj_data)
    end
    tt = regular
    time_start = time_ns()

    if line_search isa Shortstep && !isfinite(L)
        println("FATAL: Lipschitz constant not set. Prepare to blow up spectacularly.")
    end

    if line_search isa Agnostic || line_search isa Nonconvex
        println("FATAL: Lazification is not known to converge with open-loop step size strategies.")
    end

    if verbose
        println("\nLazified Conditional Gradients (Frank-Wolfe + Lazification).")
        NumType = eltype(x0)
        println(
            "EMPHASIS: $emphasis STEPSIZE: $line_search EPSILON: $epsilon MAXITERATION: $max_iteration lazy_tolerance: $lazy_tolerance TYPE: $NumType",
        )
        grad_type = typeof(gradient)
        println("GRADIENTTYPE: $grad_type CACHESIZE $cache_size GREEDYCACHE: $greedy_lazy")
        if emphasis == memory
            println("WARNING: In memory emphasis mode iterates are written back into x0!")
        end
        headers =
            ["Type", "Iteration", "Primal", "Dual", "Dual Gap", "Time", "It/sec", "Cache Size"]
        print_callback(headers, format_string, print_header=true)
    end

    if emphasis == memory && !isa(x, Union{Array,SparseArrays.AbstractSparseArray})
        if eltype(x) <: Integer
            x = copyto!(similar(x, float(eltype(x))), x)
        else
            x = copyto!(similar(x), x)
        end
    end

    if gradient === nothing
        gradient = similar(x)
    end

    # container for direction
    d = similar(x)

    while t <= max_iteration && dual_gap >= max(epsilon, eps(float(eltype(x))))

        #####################
        # managing time and Ctrl-C
        #####################
        time_at_loop = time_ns()
        if t == 0
            time_start = time_at_loop
        end
        # time is measured at beginning of loop for consistency throughout all algorithms
        tot_time = (time_at_loop - time_start) / 1e9

        if timeout < Inf
            if tot_time ≥ timeout
                if verbose
                    @info "Time limit reached"
                end
                break
            end
        end

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

        grad!(gradient, x)

        threshold = fast_dot(x, gradient) - phi / lazy_tolerance

        # go easy on the memory - only compute if really needed
        if ((mod(t, print_iter) == 0 && verbose ) || callback !== nothing)
            primal = f(x)
        end

        v = compute_extreme_point(lmo, gradient, threshold=threshold, greedy=greedy_lazy)
        tt = lazy
        if fast_dot(v, gradient) > threshold
            tt = dualstep
            dual_gap = fast_dot(x, gradient) - fast_dot(v, gradient)
            phi = min(dual_gap, phi / 2)
        end

        @emphasis(emphasis, d = x - v)

        gamma, L = line_search_wrapper(
            line_search,
            t,
            f,
            grad!,
            x,
            d,
            gradient,
            dual_gap,
            L,
            gamma0,
            linesearch_tol,
            step_lim,
            1.0,
        )

        if callback !== nothing
            state = (
                t=t,
                primal=primal,
                dual=primal - dual_gap,
                dual_gap=dual_gap,
                time=tot_time,
                cache_size=length(lmo),
                x=x,
                v=v,
                gamma=gamma
            )
            callback(state)
        end

        @emphasis(emphasis, x = x - gamma * d)

        if verbose && (mod(t, print_iter) == 0 || tt == dualstep)
            if t == 0
                tt = initial
            end
            rep = (
                st[Symbol(tt)],
                string(t),
                Float64(primal),
                Float64(primal - dual_gap),
                Float64(dual_gap),
                tot_time,
                t / tot_time,
                length(lmo),
            )
            print_callback(rep, format_string)
            flush(stdout)
        end
        t += 1
    end

    # recompute everything once for final verfication / do not record to trajectory though for now!
    # this is important as some variants do not recompute f(x) and the dual_gap regularly but only when reporting
    # hence the final computation.
    grad!(gradient, x)
    v = compute_extreme_point(lmo, gradient)
    primal = f(x)
    dual_gap = fast_dot(x, gradient) - fast_dot(v, gradient)

    if verbose
        tt = last
        tot_time = (time_ns() - time_start) / 1.0e9
        rep = (
            st[Symbol(tt)],
            string(t - 1),
            Float64(primal),
            Float64(primal - dual_gap),
            Float64(dual_gap),
            tot_time,
            t / tot_time,
            length(lmo),
        )
        print_callback(rep, format_string)
        print_callback(nothing, format_string, print_footer=true)
        flush(stdout)
    end
    return x, v, primal, dual_gap, traj_data
end

"""
    stochastic_frank_wolfe(f::StochasticObjective, lmo, x0; ...)

Stochastic version of Frank-Wolfe, evaluates the objective and gradient stochastically,
implemented through the [FrankWolfe.StochasticObjective](@ref) interface.

Keyword arguments include `batch_size` to pass a fixed `batch_size`
or a `batch_iterator` implementing
`batch_size = FrankWolfe.batchsize_iterate(batch_iterator)` for algorithms like
Variance-reduced and projection-free stochastic optimization, E Hazan, H Luo, 2016.

Similarly, a constant `momentum` can be passed or replaced by a `momentum_iterator`
implementing `momentum = FrankWolfe.momentum_iterate(momentum_iterator)`.
"""
function stochastic_frank_wolfe(
    f::StochasticObjective,
    lmo,
    x0;
    line_search::LineSearchMethod=Nonconvex(),
    L=Inf,
    gamma0=0,
    momentum_iterator=nothing,
    momentum=nothing,
    epsilon=1e-7,
    max_iteration=10000,
    print_iter=1000,
    trajectory=false,
    verbose=false,
    linesearch_tol=1e-7,
    emphasis::Emphasis=memory,
    rng=Random.GLOBAL_RNG,
    batch_size=length(f.xs) ÷ 10 + 1,
    batch_iterator=nothing,
    full_evaluation=false,
    callback=nothing,
    timeout=Inf,
    print_callback=print_callback,
)

    # format string for output of the algorithm
    format_string = "%6s %13s %14e %14e %14e %14e %14e %6i\n"

    t = 0
    dual_gap = Inf
    primal = Inf
    v = []
    x = x0
    tt = regular
    traj_data = []
    if trajectory && callback === nothing
        callback = trajectory_callback(traj_data)
    end
    time_start = time_ns()

    if line_search == Shortstep && L == Inf
        println("FATAL: Lipschitz constant not set. Prepare to blow up spectacularly.")
    end

    if line_search == FixedStep && gamma0 == 0
        println("FATAL: gamma0 not set. We are not going to move a single bit.")
    end
    if momentum_iterator === nothing && momentum !== nothing
        momentum_iterator = ConstantMomentumIterator(momentum)
    end
    if batch_iterator === nothing
        batch_iterator = ConstantBatchIterator(batch_size)
    end

    if verbose
        println("\nStochastic Frank-Wolfe Algorithm.")
        NumType = eltype(x0)
        println(
            "EMPHASIS: $emphasis STEPSIZE: $line_search EPSILON: $epsilon max_iteration: $max_iteration TYPE: $NumType",
        )
        println("GRADIENTTYPE: $(typeof(f.storage)) MOMENTUM: $(momentum_iterator !== nothing) batch policy: $(typeof(batch_iterator)) ")
        if emphasis == memory
            println("WARNING: In memory emphasis mode iterates are written back into x0!")
        end
        headers = ("Type", "Iteration", "Primal", "Dual", "Dual Gap", "Time", "It/sec", "batch size")
        print_callback(headers, format_string, print_header=true)
    end

    if emphasis == memory && !isa(x, Union{Array, SparseArrays.AbstractSparseArray})
        if eltype(x) <: Integer
            x = copyto!(similar(x, float(eltype(x))), x)
        else
            x = copyto!(similar(x), x)
        end
    end
    first_iter = true
    gradient = 0
    while t <= max_iteration && dual_gap >= max(epsilon, eps())

        #####################
        # managing time and Ctrl-C
        #####################
        time_at_loop = time_ns()
        if t == 0
            time_start = time_at_loop
        end
        # time is measured at beginning of loop for consistency throughout all algorithms
        tot_time = (time_at_loop - time_start) / 1e9

        if timeout < Inf
            if tot_time ≥ timeout
                if verbose
                    @info "Time limit reached"
                end
                break
            end
        end

        #####################
        batch_size = batchsize_iterate(batch_iterator)

        if momentum_iterator === nothing
            gradient = compute_gradient(
                f,
                x,
                rng=rng,
                batch_size=batch_size,
                full_evaluation=full_evaluation,
            )
        elseif first_iter
            gradient = copy(compute_gradient(
                f,
                x,
                rng=rng,
                batch_size=batch_size,
                full_evaluation=full_evaluation,
            ))
        else
            momentum = momentum_iterate(momentum_iterator)
            compute_gradient(
                f,
                x,
                rng=rng,
                batch_size=batch_size,
                full_evaluation=full_evaluation,
            )
            # gradient = momentum * gradient + (1 - momentum) * f.storage
            LinearAlgebra.mul!(gradient, LinearAlgebra.I, f.storage, 1-momentum, momentum)
        end
        first_iter = false

        v = compute_extreme_point(lmo, gradient)

        # go easy on the memory - only compute if really needed
        if (mod(t, print_iter) == 0 && verbose) ||
           callback !== nothing ||
           !(line_search isa Agnostic || line_search isa Nonconvex || line_search isa FixedStep)
            primal = compute_value(f, x, full_evaluation=true)
            dual_gap = fast_dot(x, gradient) - fast_dot(v, gradient)
        end

        if line_search isa Agnostic
            gamma = 2 // (2 + t)
        elseif line_search isa Nonconvex
            gamma = 1 / sqrt(t + 1)
        elseif line_search isa Shortstep
            gamma = dual_gap / (L * norm(x - v)^2)
        elseif line_search isa RationalShortstep
            rat_dual_gap = sum((x - v) .* gradient)
            gamma = rat_dual_gap // (L * sum((x - v) .^ 2))
        elseif line_search isa FixedStep
            gamma = gamma0
        end

        if callback !== nothing
            state = (
                t=t,
                primal=primal,
                dual=primal - dual_gap,
                dual_gap=dual_gap,
                time=tot_time,
                x=x,
                v=v,
                gamma=gamma,
            )
            callback(state)
        end

        @emphasis(emphasis, x = (1 - gamma) * x + gamma * v)

        if mod(t, print_iter) == 0 && verbose
            tt = regular
            if t == 0
                tt = initial
            end
            rep = (
                st[Symbol(tt)],
                string(t),
                Float64(primal),
                Float64(primal - dual_gap),
                Float64(dual_gap),
                tot_time,
                t / tot_time,
                batch_size,
            )
            print_callback(rep, format_string)
            flush(stdout)
        end
        t += 1
    end
    # recompute everything once for final verfication / no additional callback call
    # this is important as some variants do not recompute f(x) and the dual_gap regularly but only when reporting
    # hence the final computation.
    # last computation done with full evaluation for exact gradient

    (primal, gradient) = compute_value_gradient(f, x, full_evaluation=true)
    v = compute_extreme_point(lmo, gradient)
    # @show (gradient, primal)
    dual_gap = fast_dot(x, gradient) - fast_dot(v, gradient)
    if verbose
        tt = last
        tot_time = (time_ns() - time_start) / 1.0e9
        rep = (
            st[Symbol(tt)],
            string(t - 1),
            Float64(primal),
            Float64(primal - dual_gap),
            Float64(dual_gap),
            tot_time,
            t / tot_time,
            batch_size
        )
        print_callback(rep, format_string)
        print_callback(nothing, format_string, print_footer=true)
        flush(stdout)
    end
    return x, v, primal, dual_gap, traj_data
end

ext: jl


text:
using NSDE
using Test

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


text:
include("fpsquareOmega.jl")

# A projective point on the curve.
# The curve itself is not explicitly defined
struct TwiHesPoint{T <: Integer}
    X::Fp2Elem{T}
    Y::Fp2Elem{T}
    Z::Fp2Elem{T}
end

# Element-Constructor
# NB! Assumes point is valid! Checking everytime you construct a new point is very slow....
function EC(X::Fp2Elem{T}, Y::Fp2Elem{T}, Z::Fp2Elem{T}) where {T <: Integer}
    return TwiHesPoint{T}(X, Y, Z)
end

############################
#
# Printing
#
############################

function show(io::IO, P::TwiHesPoint{T}) where {T <: Integer}
    print(io, "($(P.X) : $(P.Y) : $(P.Z))")
end

############################
#
# Some helper functions
#
############################

function isValid(P::TwiHesPoint{T}, a::Fp2Elem{T}, d::Fp2Elem{T}) where {T <: Integer}
    return a*P.X^3 + P.Y^3 + P.Z^3 == d*P.X*P.Y*P.Z && !isZeroCoordinate(P.X,P.Y,P.Z)
end

function isZeroCoordinate(X::Fp2Elem{T}, Y::Fp2Elem{T}, Z::Fp2Elem{T}) where {T <: Integer}
    return (X,Y,Z)==(zero(X.field),zero(X.field),zero(X.field))
end

function identity(F::Fp2{T}) where {T <: Integer}
    return EC(zero(F),-one(F),one(F))
end

function isIdentity(P::TwiHesPoint{T}) where {T <: Integer}
    return P == identity(P.X.field)
end

function ==(P::TwiHesPoint{T}, Q::TwiHesPoint{T}) where {T <: Integer}
    P1 = normalized(P)
    Q1 = normalized(Q)
    return P1.X == Q1.X && P1.Y == Q1.Y && P1.Z == Q1.Z
end

# Prefers Z = 1, then Y = 1
function normalized(P::TwiHesPoint{T}) where {T <: Integer}
    if !iszero(P.Z)
        Zinv = inv(P.Z)
        return EC(P.X*Zinv, P.Y*Zinv, one(P.X.field))
    else
        return EC(P.X*inv(P.Y), one(P.X.field), zero(P.X.field))
    end
end

# Recover curve coefficient d from a and a point P, with [3]P != O
function recover_d(P::TwiHesPoint{T}, a::Fp2Elem{T}) where {T <: Integer}
    return (a*P.X*P.X*P.X + P.Y*P.Y*P.Y + P.Z*P.Z*P.Z)*inv(P.X*P.Y*P.Z)
end

############################
#
# Basic operators
#
############################

# Inverse point
# Input: A point P
# Output: -P
function -(P::TwiHesPoint{T}) where {T <: Integer}
    return EC(P.X, P.Z, P.Y)
end

# Addition based on the rotated addition law.
# Input: Points P, Q on curve E(a,d), and curve-parameter a
# Output: P+Q
function add(P::TwiHesPoint{T}, Q::TwiHesPoint{T}, a::Fp2Elem{T}) where {T <: Integer}
    A = P.X*Q.Z
    B = P.Z*Q.Z
    C = P.Y*Q.X
    D = P.Y*Q.Y
    E = P.Z*Q.Y
    F = a*P.X*Q.X
    G = (D + B)*(A - C)
    H = (D - B)*(A + C)
    J = (D + F)*(A - E)
    K = (D - F)*(A + E)
    X = G - H
    Y = K - J
    Z = J + K - G - H - 2*(B - F)*(C + E)
    return EC(X, Y, Z)
end

# This function uses the simple double and add.
# Input: integer n and point P
# Output: [n]P
function mul(n::Integer, P::TwiHesPoint{T}, a::Fp2Elem{T}) where {T <: Integer}
    if n == 1
        return P
    end
    if n == 0
        return identity(a.field)
    end
    if n < 0
        P = -P
        n = -n
    end
    bit = T(1) << (ndigits(n, base = 2) - 1)
    Q = P
    bit >>= 1
    while bit != 0
        Q = double(Q, a)
        if (bit & n) != 0
            Q = add(Q, P, a)
        end
        bit >>=1
    end
    return Q
end

# Doubling of a point
# Costs 6M + 3S + 1Ma
# Can be improved, but we wish to avoid the need for d
# Input: A point P on E(a,d), and curve parameter a
# Output: [2]P
function double(P::TwiHesPoint{T}, a::Fp2Elem{T}) where {T <: Integer}
    X, Y, Z = P.X, P.Y, P.Z
    #X = Z^3*X - Y^3*X
    #Y = Y^3*Z - a*X^3*Z
    #Z = a*X^3*Y - Z^3*Y
    X3 = P.X*P.X*P.X
    Y3 = P.Y*P.Y*P.Y
    Z3 = P.Z*P.Z*P.Z
    aX3 = a*X3
    X_new = (Z3 - Y3)*P.X
    Y_new = (Y3 - aX3)*P.Z
    Z_new = (aX3 - Z3)*P.Y
    return EC(X_new, Y_new, Z_new)
end

#Repeated doubling
#Input point P, integer e
#Output [2^e]P
function doublePow(P::TwiHesPoint{T}, a::Fp2Elem{T}, e::Integer) where {T <: Integer}
    Q = P
    for i in 1:e
        Q = double(Q, a)
    end
    return Q
end

############################
#
# Isogeny-related
#
############################

#Input: Curve a and d from E(a,d) in twisted Hessian form
#Output: j(E)
function jInvariant(a::Fp2Elem{T}, d::Fp2Elem{T}) where {T <: Integer}
    d3 = d*d*d
    o = (216*a+d3)*d
    u = (d3-27*a)
    over = o*o*o
    under = a*u*u*u
    return over*inv(under)
end

#input: A point P, and precomputed values a2 = a^2, b2 = b^2 ab = a*b from the kernel (a,b,b)
#Output: Phi(P)
function iso2eval(P::TwiHesPoint{T}, a2::Fp2Elem{T}, b2::Fp2Elem{T}, ab::Fp2Elem{T}) where {T <: Integer}
    X2 = P.X*P.X
    Y2 = P.Y*P.Y
    Z2 = P.Z*P.Z
    XY = P.X*P.Y
    XZ = P.X*P.Z
    YZ = P.Y*P.Z
    X_new = P.X*(X2*b2 - a2*YZ)
    Y_new = P.Y*(Z2*ab - b2*XY)
    Z_new = P.Z*(Y2*ab - b2*XZ)
    return EC(X_new, Y_new, Z_new)
end

# Some precomputed values for faster 2-isogenies
# Input: Point G of order 2
# Output: Precomputed values for evaluating isogeny generated by G
function iso2precomp(G::TwiHesPoint{T}) where {T <: Integer}
    return G.X*G.X, G.Y*G.Y, G.X*G.Y
end

# 3-isogeny evaluation when the kernel is <1 : -c : 0)>
# Input: Point P, and cb and b such that c^3 = a
# Output: Phi(P)
function iso3eval_c2(P::TwiHesPoint{T}, c::Fp2Elem{T}) where {T <: Integer}
    YZ = P.Y*P.Z
    cX = c*P.X
    c2X2 = cX*cX
    X_new = P.X*YZ
    Y_new = c2X2*P.Z + cX*P.Y*P.Y + YZ*P.Z
    Z_new = c2X2*P.Y + cX*P.Z*P.Z + YZ*P.Y
    return EC(X_new, Y_new, Z_new)
end

# Evaluating the isomorphism \varphi: H(a, d) -> H(1, d/a^(1/3))
# Input A point P, and a cuberoot of A
# Output \varphi(P)
function isoHessianEval(P::TwiHesPoint{T}, cuberoota::Fp2Elem{T}) where {T <: Integer}
    X, Y, Z = P.X, P.Y, P.Z
    return EC(cuberoota*X, Y, Z)
end

###############################
#
# Optimal isogeny strategies
#
###############################

#Input: A generator point G of order 2^e, and two other points P, Q, curve parameter a, isogeny degree e, and an optimal strategy S
#Output: E/<G>, and phi(P), phi(Q)
function iso2pow_opt(G::TwiHesPoint{T}, P::TwiHesPoint{T}, Q::TwiHesPoint{T}, a::Fp2Elem{T}, e::Integer, Strat::Vector{Int64}) where {T <: Integer}
    Pi = P
    Qi = Q
    ai = a
    Si = Vector{Tuple{Integer,TwiHesPoint{T}}}()
    push!(Si, (e, G))
    i = 1
    while !isempty(Si)
        (h, Ki) = pop!(Si)
        if h == 1
            ai = ai*ai
            a2, b2, ab = iso2precomp(Ki)
            S_new = Vector{Tuple{Integer,TwiHesPoint{T}}}()
            while !isempty(Si)
                (h, Xi) = popfirst!(Si)
                push!(S_new, (h-1, iso2eval(Xi, a2, b2, ab)))
            end
            Si = S_new
            Pi = iso2eval(Pi, a2, b2, ab)
            Qi = iso2eval(Qi, a2, b2, ab)
        else
            strat_i = Strat[i]
            if 0 < strat_i < h
                push!(Si, (h, Ki))
                push!(Si, (h-strat_i, doublePow(Ki, ai, strat_i)))
                i += 1
            else
                error("Strategy seems to be invalid")
            end
        end
    end
    return Pi, Qi, ai
end

#Second part, no points to carry
#Input: A generator point G of order 2^e, curve parameter a, isogeny degree e, and an optimal strategy S
#Output: parameters a' and d' of ending curve
function iso2pow_opt(G::TwiHesPoint{T}, a::Fp2Elem{T}, e::Integer, Strat::Vector{Int64}) where {T <: Integer}
    ai = a
    Si = Vector{Tuple{Integer,TwiHesPoint{T}}}()
    push!(Si, (e, G))
    i = 1
    while !isempty(Si)
        (h, Ki) = pop!(Si)
        if h == 1
            if isempty(Si)
                #Second to last curve, must recover d here, and use it to find d on the final curve
                di = recover_d(Ki, ai)
                a_final = ai*ai
                d_final = Ki.Y*(6*Ki.Y-Ki.X*di)*inv(Ki.X*Ki.X)
                return a_final, d_final
            end
            ai = ai*ai
            a2, b2, ab = iso2precomp(Ki)
            S_new = Vector{Tuple{Integer,TwiHesPoint{T}}}()
            while !isempty(Si)
                (h, Xi) = popfirst!(Si)
                push!(S_new, (h-1, iso2eval(Xi, a2, b2, ab)))
            end
            Si = S_new
        else
            strat_i = Strat[i]
            if 0 < strat_i < h
                push!(Si, (h, Ki))
                push!(Si, (h-strat_i, doublePow(Ki, ai, strat_i)))
                i += 1
            else
                error("Strategy seems to be invalid")
            end
        end
    end
end

function Xnormalized(P::TwiHesPoint{T}) where {T <: Integer}
    X, Y, Z = P.X, P.Y, P.Z
    F = X.field
    if !iszero(X)
        Xinv = inv(X)
        return EC(one(F), Y*Xinv, Z*Xinv)
    else
        Zinv = inv(Z)
        return EC(zero(F), Y*Zinv, one(F))
    end
end


text:
"""
    struct Crane

Struct representing a rail-mounted gantry crane as modeled in the GCSPET. Straight forward model of crane-movement at constant `speed` (`default = 1`), with a zone of `safety` (`default = 1`) wide that must be kept free on either side of the crane.

# Functionality
- id(q): identifier and position in the order of cranes along the non-crossing axis
- speed(q): movement speed of the crane
- starting_position(q): initial position of the crane along the non-crossing axis
- zone(q): minimum and maximum position of the crane along the non-crossing axis
- safety(q): minimal distance that must be kept free on either side of the crane
"""
struct Crane
    id::Int
    l0::Int
    speed::Int
    zone::UnitRange{Int64}
    safety::Int
end

id(c::Crane) = c.id
starting_position(c::Crane) = c.l0
speed(c::Crane) = c.speed
zone(c::Crane) = c.zone
safety(c::Crane) = c.safety


text:
function constant!(e::EGraph, x::Symbol)
    t = Term(x, [])
    push!(e, t)
end

term!(e::EGraph, f::Symbol) = (args...) -> begin
  t = Term(f, collect(args))
  push!(e,  t)
end

using LightGraphs
using GraphPlot
using Colors

function graph(e::EGraph)
    nverts = length(e.memo)
    g = SimpleDiGraph(nverts)
    vertmap = Dict([ t => n for (n, (t, id)) in enumerate(e.memo)])
    #=for (id, cls) in e.classes
        for t1 in cls.nodes
            for (t2,_) in cls.parents
                add_edge!(g, vertmap[t2], vertmap[t1])
            end
        end
    end =#
    for (n,(t,id)) in enumerate(e.memo)
        for n2 in t.args
            for t2 in e.classes[n2].nodes
                add_edge!(g, n, vertmap[t2])
            end
        end
    end
    nodelabel = [ t.head for (t, id) in e.memo]
    classmap = Dict([ (id,n)  for (n,id) in enumerate(Set([ id.id for (t, id) in e.memo]))])
    nodecolor = [classmap[id.id] for (t, id) in e.memo]

    return g, nodelabel, nodecolor
end

function graphplot(e::EGraph)
    g, nodelabel,nodecolor = graph(e)

    # Generate n maximally distinguishable colors in LCHab space.
    nodefillc = distinguishable_colors(maximum(nodecolor), colorant"blue")
    gplot(g, nodelabel=nodelabel, nodefillc=nodefillc[nodecolor])
end


text:
function get_affected_values(ff::AbstractAffectFeedforward)
    return ff.affected_values
end

function attach_feedforward!(model, ff::T) where {T <: AbstractAffectFeedforward}
    if !isempty(model.feedforwards)
        ff_k = [get_optimization_container_key(v) for v in model.feedforwards if isa(v, T)]
        if get_optimization_container_key(ff) ∈ ff_k
            return
        end
    end
    push!(model.feedforwards, ff)
    return
end

function get_component_type(ff::AbstractAffectFeedforward)
    return get_component_type(get_optimization_container_key(ff))
end

"""
Adds an upper bound constraint to a variable.
"""
struct UpperBoundFeedforward <: AbstractAffectFeedforward
    optimization_container_key::OptimizationContainerKey
    affected_values::Vector
    function UpperBoundFeedforward(;
        component_type::Type{<:PSY.Component},
        source::Type{T},
        affected_values::Vector{DataType},
        meta=CONTAINER_KEY_EMPTY_META,
    ) where {T}
        values_vector = Vector(undef, length(affected_values))
        for (ix, v) in enumerate(affected_values)
            if v <: VariableType
                values_vector[ix] =
                    get_optimization_container_key(v(), component_type, meta)
            else
                error(
                    "UpperBoundFeedforward is only compatible with VariableType affected values",
                )
            end
        end
        new(get_optimization_container_key(T(), component_type, meta), values_vector)
    end
end

get_default_parameter_type(::UpperBoundFeedforward, _) = UpperBoundValueParameter()
get_optimization_container_key(ff::UpperBoundFeedforward) = ff.optimization_container_key

"""
Adds a lower bound constraint to a variable.
"""
struct LowerBoundFeedforward <: AbstractAffectFeedforward
    optimization_container_key::OptimizationContainerKey
    affected_values::Vector{<:OptimizationContainerKey}
    function LowerBoundFeedforward(;
        component_type::Type{<:PSY.Component},
        source::Type{T},
        affected_values::Vector{DataType},
        meta=CONTAINER_KEY_EMPTY_META,
    ) where {T}
        values_vector = Vector{VariableKey}(undef, length(affected_values))
        for (ix, v) in enumerate(affected_values)
            if v <: VariableType
                values_vector[ix] =
                    get_optimization_container_key(v(), component_type, meta)
            else
                error(
                    "LowerBoundFeedforward is only compatible with VariableType affected values",
                )
            end
        end
        new(get_optimization_container_key(T(), component_type, meta), values_vector)
    end
end

get_default_parameter_type(::LowerBoundFeedforward, _) = LowerBoundValueParameter()
get_optimization_container_key(ff::LowerBoundFeedforward) = ff.optimization_container_key

"""
Adds a constraint to make the bounds of a variable 0.0. Effectively allows to "turn off" a value.
"""
struct SemiContinuousFeedforward <: AbstractAffectFeedforward
    optimization_container_key::OptimizationContainerKey
    affected_values::Vector{<:OptimizationContainerKey}
    function SemiContinuousFeedforward(;
        component_type::Type{<:PSY.Component},
        source::Type{T},
        affected_values::Vector{DataType},
        meta=CONTAINER_KEY_EMPTY_META,
    ) where {T}
        values_vector = Vector{VariableKey}(undef, length(affected_values))
        for (ix, v) in enumerate(affected_values)
            if v <: VariableType
                values_vector[ix] =
                    get_optimization_container_key(v(), component_type, meta)
            else
                error(
                    "SemiContinuousFeedforward is only compatible with VariableType affected values",
                )
            end
        end
        new(get_optimization_container_key(T(), component_type, meta), values_vector)
    end
end

get_default_parameter_type(::SemiContinuousFeedforward, _) = OnStatusParameter()
get_optimization_container_key(f::SemiContinuousFeedforward) = f.optimization_container_key

function has_semicontinuous_feedforward(
    model::DeviceModel,
    ::Type{T},
)::Bool where {T <: Union{VariableType, ExpressionType}}
    if isempty(model.feedforwards)
        return false
    end
    sc_feedforwards = [x for x in model.feedforwards if isa(x, SemiContinuousFeedforward)]
    if isempty(sc_feedforwards)
        return false
    end

    keys = get_affected_values(sc_feedforwards[1])

    return T ∈ get_entry_type.(keys)
end

function has_semicontinuous_feedforward(
    model::DeviceModel,
    ::Type{T},
)::Bool where {T <: Union{ActivePowerRangeExpressionUB, ActivePowerRangeExpressionLB}}
    return has_semicontinuous_feedforward(model, ActivePowerVariable)
end

"""
Adds a constraint to limit the sum of a variable over the number of periods to the source value
"""
struct EnergyLimitFeedforward <: AbstractAffectFeedforward
    optimization_container_key::OptimizationContainerKey
    affected_values::Vector{<:OptimizationContainerKey}
    number_of_periods::Int
    function EnergyLimitFeedforward(;
        component_type::Type{<:PSY.Component},
        source::Type{T},
        affected_values::Vector{DataType},
        number_of_periods::Int,
        meta=CONTAINER_KEY_EMPTY_META,
    ) where {T}
        values_vector = Vector{VariableKey}(undef, length(affected_values))
        for (ix, v) in enumerate(affected_values)
            if v <: VariableType
                values_vector[ix] =
                    get_optimization_container_key(v(), component_type, meta)
            else
                error(
                    "EnergyLimitFeedforward is only compatible with VariableType or ParamterType affected values",
                )
            end
        end
        new(
            get_optimization_container_key(T(), component_type, meta),
            values_vector,
            number_of_periods,
        )
    end
end

get_default_parameter_type(::EnergyLimitFeedforward, _) = EnergyLimitParameter()
get_optimization_container_key(ff) = ff.optimization_container_key
get_number_of_periods(ff) = ff.number_of_periods

"""
Fixes a Variable or Parameter Value in the model. Is the only Feed Forward that can be used
with a Parameter or a Variable as the affected value.
"""
struct FixValueFeedforward <: AbstractAffectFeedforward
    optimization_container_key::OptimizationContainerKey
    affected_values::Vector
    function FixValueFeedforward(;
        component_type::Type{<:PSY.Component},
        source::Type{T},
        affected_values::Vector{DataType},
        meta=CONTAINER_KEY_EMPTY_META,
    ) where {T}
        values_vector = Vector(undef, length(affected_values))
        for (ix, v) in enumerate(affected_values)
            if v <: VariableType || v <: ParameterType
                values_vector[ix] =
                    get_optimization_container_key(v(), component_type, meta)
            else
                error(
                    "UpperBoundFeedforward is only compatible with VariableType affected values",
                )
            end
        end
        new(get_optimization_container_key(T(), component_type, meta), values_vector)
    end
end

get_default_parameter_type(::FixValueFeedforward, _) = FixValueParameter()
get_optimization_container_key(ff::FixValueFeedforward) = ff.optimization_container_key

"""
Adds a constraint to enforce a minimum energy level target with a slack variable associated witha penalty term.
"""
struct EnergyTargetFeedforward <: AbstractAffectFeedforward
    optimization_container_key::OptimizationContainerKey
    affected_values::Vector{<:OptimizationContainerKey}
    target_period::Int
    penalty_cost::Float64
    function EnergyTargetFeedforward(;
        component_type::Type{<:PSY.Component},
        source::Type{T},
        affected_values::Vector{DataType},
        target_period::Int,
        penalty_cost::Float64,
        meta=CONTAINER_KEY_EMPTY_META,
    ) where {T}
        values_vector = Vector{VariableKey}(undef, length(affected_values))
        for (ix, v) in enumerate(affected_values)
            if v <: VariableType
                values_vector[ix] =
                    get_optimization_container_key(v(), component_type, meta)
            else
                error(
                    "EnergyTargetFeedforward is only compatible with VariableType or ParamterType affected values",
                )
            end
        end
        new(
            get_optimization_container_key(T(), component_type, meta),
            values_vector,
            target_period,
            penalty_cost,
        )
    end
end

get_default_parameter_type(::EnergyTargetFeedforward, _) = EnergyTargetParameter()
get_optimization_container_key(ff::EnergyTargetFeedforward) = ff.optimization_container_key


text:
function boundary_traction_operator(
    basis,
    quad,
    normal,
    stiffness,
    dim,
    facedetjac,
    jac,
    vectosymmconverter,
)

    numqp = length(quad)
    normals = repeat(normal, inner = (1, numqp))
    scalearea = repeat([facedetjac], numqp)
    return -1.0 * surface_traction_operator(
        basis,
        quad,
        quad,
        normals,
        stiffness,
        dim,
        scalearea,
        jac,
        vectosymmconverter,
    )
end

function assemble_boundary_traction_operator!(
    sysmatrix,
    basis,
    facequads,
    stiffness,
    mesh,
    onboundary,
)

    facemidpoints = reference_face_midpoints()
    nfaces = number_of_faces_per_cell(facequads)
    facedetjac = face_determinant_jacobian(mesh)
    jac = jacobian(mesh)
    dim = dimension(mesh)
    vectosymmconverter = vector_to_symmetric_matrix_converter()
    normals = reference_face_normals()
    uniformquads = uniform_face_quadratures(facequads)

    faceids = 1:nfaces

    uniformtop1 = [
        vec(
            boundary_traction_operator(
                basis,
                q,
                n,
                stiffness[+1],
                dim,
                d,
                jac,
                vectosymmconverter,
            ),
        ) for (q, n, d) in zip(uniformquads, normals, facedetjac)
    ]
    uniformtop2 = [
        vec(
            boundary_traction_operator(
                basis,
                q,
                n,
                stiffness[-1],
                dim,
                d,
                jac,
                vectosymmconverter,
            ),
        ) for (q, n, d) in zip(uniformquads, normals, facedetjac)
    ]
    uniformtop = [uniformtop1, uniformtop2]

    ncells = number_of_cells(mesh)

    for cellid = 1:ncells
        cellsign = cell_sign(mesh, cellid)
        cellmap1 = cell_map(mesh, +1, cellid)
        cellmap2 = cell_map(mesh, -1, cellid)

        check_cellsign(cellsign)
        for faceid in faceids
            if cell_connectivity(mesh, faceid, cellid) == 0 && (
                onboundary(cellmap1(facemidpoints[faceid])) ||
                onboundary(cellmap2(facemidpoints[faceid]))
            )

                if cellsign == +1 || cellsign == -1
                    row = cell_sign_to_row(cellsign)
                    nodeids = nodal_connectivity(mesh, cellsign, cellid)
                    assemble_cell_matrix!(
                        sysmatrix,
                        nodeids,
                        dim,
                        uniformtop[row][faceid],
                    )
                else
                    assemble_cut_cell_boundary_traction_operator!(
                        sysmatrix,
                        basis,
                        facequads,
                        normals,
                        stiffness,
                        mesh,
                        +1,
                        faceid,
                        cellid,
                        facedetjac,
                        jac,
                        vectosymmconverter,
                    )
                    assemble_cut_cell_boundary_traction_operator!(
                        sysmatrix,
                        basis,
                        facequads,
                        normals,
                        stiffness,
                        mesh,
                        -1,
                        faceid,
                        cellid,
                        facedetjac,
                        jac,
                        vectosymmconverter,
                    )
                end

            end
        end
    end
end

function assemble_cut_cell_boundary_traction_operator!(
    sysmatrix,
    basis,
    facequads,
    normals,
    stiffness,
    mesh,
    cellsign,
    faceid,
    cellid,
    facedetjac,
    jac,
    vectosymmconverter,
)

    dim = dimension(mesh)
    operator = vec(
        boundary_traction_operator(
            basis,
            facequads[cellsign, faceid, cellid],
            normals[faceid],
            stiffness[cellsign],
            dim,
            facedetjac[faceid],
            jac,
            vectosymmconverter,
        ),
    )
    nodeids = nodal_connectivity(mesh, cellsign, cellid)
    assemble_cell_matrix!(sysmatrix, nodeids, dim, operator)
end

function boundary_mass_operator(basis, quad, dim, scale)
    numqp = length(quad)
    scalearea = repeat([scale], numqp)
    return mass_matrix(basis, quad, quad, dim, scalearea)
end

function assemble_boundary_mass_operator!(
    sysmatrix,
    basis,
    facequads,
    mesh,
    onboundary,
    penalty,
)

    uniformquads = uniform_face_quadratures(facequads)
    facedetjac = face_determinant_jacobian(mesh)
    nfaces = number_of_faces_per_cell(facequads)
    dim = dimension(mesh)
    facemidpoints = reference_face_midpoints()

    faceids = 1:nfaces

    uniformmassops = [
        vec(boundary_mass_operator(basis, q, dim, penalty * d)) for
        (q, d) in zip(uniformquads, facedetjac)
    ]

    ncells = number_of_cells(mesh)

    for cellid = 1:ncells
        cellsign = cell_sign(mesh, cellid)
        cellmap1 = cell_map(mesh, +1, cellid)
        cellmap2 = cell_map(mesh, -1, cellid)

        check_cellsign(cellsign)
        for faceid in faceids
            if cell_connectivity(mesh, faceid, cellid) == 0 && (
                onboundary(cellmap1(facemidpoints[faceid])) ||
                onboundary(cellmap2(facemidpoints[faceid]))
            )

                if cellsign == +1 || cellsign == -1
                    nodeids = nodal_connectivity(mesh, cellsign, cellid)
                    assemble_cell_matrix!(
                        sysmatrix,
                        nodeids,
                        dim,
                        uniformmassops[faceid],
                    )
                else
                    assemble_cut_cell_boundary_mass_operator!(
                        sysmatrix,
                        basis,
                        facequads,
                        mesh,
                        +1,
                        faceid,
                        cellid,
                        facedetjac,
                        penalty,
                    )
                    assemble_cut_cell_boundary_mass_operator!(
                        sysmatrix,
                        basis,
                        facequads,
                        mesh,
                        -1,
                        faceid,
                        cellid,
                        facedetjac,
                        penalty,
                    )
                end
            end
        end
    end
end

function assemble_cut_cell_boundary_mass_operator!(
    sysmatrix,
    basis,
    facequads,
    mesh,
    cellsign,
    faceid,
    cellid,
    facedetjac,
    penalty,
)

    dim = dimension(mesh)
    operator = vec(
        boundary_mass_operator(
            basis,
            facequads[cellsign, faceid, cellid],
            dim,
            penalty * facedetjac[faceid],
        ),
    )
    nodeids = nodal_connectivity(mesh, cellsign, cellid)
    assemble_cell_matrix!(sysmatrix, nodeids, dim, operator)
end

function assemble_boundary_displacement_rhs!(
    sysrhs,
    bcfunc,
    basis,
    facequads,
    mesh,
    onboundary,
    penalty,
)

    facemidpoints = reference_face_midpoints()
    nfaces = number_of_faces_per_cell(facequads)
    facedetjac = face_determinant_jacobian(mesh)
    ncells = number_of_cells(mesh)
    faceids = 1:nfaces

    for cellid = 1:ncells
        cellsign = cell_sign(mesh, cellid)
        cellmap1 = cell_map(mesh, +1, cellid)
        cellmap2 = cell_map(mesh, -1, cellid)

        check_cellsign(cellsign)
        for faceid in faceids
            if cell_connectivity(mesh, faceid, cellid) == 0 && (
                onboundary(cellmap1(facemidpoints[faceid])) ||
                onboundary(cellmap2(facemidpoints[faceid]))
            )

                if cellsign == +1 || cellsign == 0
                    assemble_face_boundary_displacement_rhs!(
                        sysrhs,
                        bcfunc,
                        basis,
                        facequads,
                        mesh,
                        +1,
                        faceid,
                        cellid,
                        facedetjac,
                        penalty,
                    )
                end
                if cellsign == -1 || cellsign == 0
                    assemble_face_boundary_displacement_rhs!(
                        sysrhs,
                        bcfunc,
                        basis,
                        facequads,
                        mesh,
                        -1,
                        faceid,
                        cellid,
                        facedetjac,
                        penalty,
                    )
                end
            end
        end
    end
end

function assemble_face_boundary_displacement_rhs!(
    sysrhs,
    bcfunc,
    basis,
    facequads,
    mesh,
    cellsign,
    faceid,
    cellid,
    facedetjac,
    penalty,
)

    dim = dimension(mesh)
    cellmap = cell_map(mesh, cellsign, cellid)
    rhsvals =
        penalty * linear_form(
            bcfunc,
            basis,
            facequads[cellsign, faceid, cellid],
            cellmap,
            dim,
            facedetjac[faceid],
        )
    nodeids = nodal_connectivity(mesh, cellsign, cellid)
    assemble_cell_rhs!(sysrhs, nodeids, dim, rhsvals)
end

function assemble_penalty_displacement_bc!(
    sysmatrix,
    sysrhs,
    bcfunc,
    basis,
    facequads,
    stiffness,
    mesh,
    onboundary,
    penalty,
)

    assemble_boundary_traction_operator!(
        sysmatrix,
        basis,
        facequads,
        stiffness,
        mesh,
        onboundary,
    )
    assemble_boundary_mass_operator!(
        sysmatrix,
        basis,
        facequads,
        mesh,
        onboundary,
        penalty,
    )
    assemble_boundary_displacement_rhs!(
        sysrhs,
        bcfunc,
        basis,
        facequads,
        mesh,
        onboundary,
        penalty,
    )
end


text:
###############################################################################
#                                  GRAD CONTEXT                               #
###############################################################################

struct GradCtx
    # map from primal var to its pullback var
    pullbacks::Dict{Variable,Variable}
    # map from primal var to its derivative var
    derivs::Dict{Variable,Variable}
end

GradCtx() = GradCtx(Dict(), Dict())

function rebind_context!(tape::Tape{GradCtx}, st::Dict)
    for (v, dv) in tape.c.derivs
        rebind!(tape, v, st)
        rebind!(tape, dv, st)
    end
    for (v, dv) in tape.c.pullbacks
        rebind!(tape, v, st)
        rebind!(tape, dv, st)
    end
end

const BASE_CTX = BaseCtx()
const CR_CTX = ChainRulesCtx()
const YOTA_RULE_CONFIG = YotaRuleConfig()

function isprimitive(::GradCtx, f, args...)
    return (isprimitive(BASE_CTX, f, args...) || isprimitive(CR_CTX, f, args...))
end

"""
    record_primitive!(tape::Tape{GradCtx}, v_fargs...)

Replace ChainRules primitives `f(args...)` with a sequence:

    rr = push!(tape, mkcall(rrule, f, args...))   # i.e. rrule(f, args...)
    val = push!(tape, mkcall(getfield, rr, 1)     # extract value
    pb = push!(tape, mkcall(getfield, rr, 2)      # extract pullback
"""
function record_primitive!(tape::Tape{GradCtx}, v_fargs...)
    v_f, v_args... = v_fargs
    f, args... = [v isa V ? tape[v].val : v for v in v_fargs]
    if isprimitive(CR_CTX, f, args...)
        rr_op = (is_kwfunc(f) ?
                    mkcall(Core.kwfunc(rrule), v_args[1], rrule, YOTA_RULE_CONFIG, v_args[2:end]...) :
                    mkcall(rrule, YOTA_RULE_CONFIG, v_f, v_args...))
        @assert rr_op.val !== nothing "rrule($f, ...) returned nothing"
        v_rr = push!(tape, rr_op)
        v_val = push!(tape, mkcall(_getfield, v_rr, 1))
        v_pb = push!(tape, mkcall(_getfield, v_rr, 2))
        tape.c.pullbacks[v_val] = v_pb
        return v_val
    else
        return push!(tape, mkcall(v_fargs...))
    end
end

###############################################################################
#                            BCAST GRAD CONTEXT                               #
###############################################################################

struct BcastGradCtx
    inner
end

# get_static_params is broken for BcastGradCtx, so turning off
# this feature for now
Umlaut.get_static_params(::Tracer{BcastGradCtx}, v_fargs) = Core.svec([])

function record_or_recurse!(t::Tracer{BcastGradCtx}, v_fargs...)
    fargs = [v isa V ? t.tape[v].val : v for v in v_fargs]
    # global STATE = (t, v_fargs)
    v_f, v_args... = v_fargs
    f, args... = [v isa V ? t.tape[v].val : v for v in v_fargs]
    return if isprimitive(t.tape.c.inner, fargs...)
        # push!(t.tape, mkcall(broadcasted, vs...))
        rr_op = (is_kwfunc(f) ?
                    mkcall(Core.kwfunc(bcast_rrule), v_args[1], bcast_rrule, YOTA_RULE_CONFIG, broadcasted, v_args[2:end]...) :
                    mkcall(bcast_rrule, YOTA_RULE_CONFIG, broadcasted, v_f, v_args...))
        v_rr = push!(t.tape, rr_op)
        v_val = push!(t.tape, mkcall(_getfield, v_rr, 1))
        v_pb = push!(t.tape, mkcall(_getfield, v_rr, 2))
        t.tape.c.inner.pullbacks[v_val] = v_pb
        return v_val
    else
        types = map(eltype, fargs[2:end])
        trace!(t, getcode(fargs[1], types), v_fargs...)
    end
end

###############################################################################
#                                   GRAD                                      #
###############################################################################

getderiv(tape::Tape, v::Variable) = get(tape.c.derivs, bound(tape, v), nothing)
setderiv!(tape::Tape, x::Variable, dx::Variable) = (
    tape.c.derivs[bound(tape, x)] = bound(tape, dx)
)
hasderiv(tape::Tape, v::Variable) = getderiv(tape, v) !== nothing

function set_or_add_deriv!(tape::Tape, x::Variable, dx::Variable)
    if !hasderiv(tape, x)
        setderiv!(tape, x, dx)
    else
        old_dx = getderiv(tape, x)
        if tape[dx].val isa Tangent || tape[old_dx].val isa Tangent
            new_dx = push!(tape, mkcall(+, dx, old_dx))
        else
            new_dx = push!(tape, mkcall(broadcast, +, dx, old_dx))
        end
        setderiv!(tape, x, new_dx)
    end
end

"""
Collect variables that we need to step through during the reverse pass.
The returned vector is already deduplicated and reverse-sorted
"""
function todo_list(tape::Tape{GradCtx}, y=tape.result)
    y_orig = y
    y_fargs = [tape[y].fn; tape[y].args...]
    is_rrule_based = haskey(tape.c.pullbacks, y)
    if is_rrule_based
        # use rrule instead
        y = tape[y].args[1]
        y_fargs = is_kwfunc(y._op.fn) ? tape[y].args[3:end] : tape[y].args
    end
    y_todo = [x for x in y_fargs if x isa V && tape[x] isa Call]
    x_todos = [todo_list(tape, x) for x in y_todo]
    # include y itself (original), its parents and their parents recursively
    todo = [[y_orig]; y_todo; vcat(x_todos...)]
    # deduplicate
    todo = collect(Set([bound(tape, v) for v in todo]))
    todo = sort(todo, by=v->v.id, rev=true)
    return todo
end

"""
Make a single step of backpropagation.
"""
function step_back!(tape::Tape, y::Variable)
    @debug "step_back!() for $(tape[y])"
    if !hasderiv(tape, y)
        @debug "No derivative found for y = $y, stop propagating the gradient in this path"
        return
    end
    dy = tape.c.derivs[y]
    if haskey(tape.c.pullbacks, y)
        # ChainRules
        pb = tape.c.pullbacks[y]
        @debug "Found pullback: $pb"
        dxs = push!(tape, mkcall(pb, dy))
        # propage derivs to rrule variable
        rr = tape[y].args[1]
        y_fargs = is_kwfunc(rr._op.fn) ? tape[rr].args[4:end] : tape[rr].args[2:end]
    else
        sig_str = join(["::$T" for T in Umlaut.call_signature(tape, tape[y]).parameters], ", ")
        error("No deriative rule found for op $(tape[y]), " *
              "try defining it using \n\n\tChainRulesCore.rrule($sig_str) = ...\n")
    end
    for (i, x) in enumerate(y_fargs)
        if x isa V
            dx = push!(tape, mkcall(getfield, dxs, i))
            @debug "Updating derivative: $x -> $dx"
            set_or_add_deriv!(tape, x, dx)
        end
    end
end

const BACKPROP_STATE = Ref{Tuple}()

"""
Backpropagate through the tape, record derivatives as new operations.
"""
function back!(tape::Tape; seed=1)
    # z - final variable (usually a loss)
    # y - resulting variable of current op
    # x - dependencies of y
    # dy - derivative of z w.r.t. y
    z = tape.result
    if (seed == 1) && (ndims(tape[z].val) > 0)
        error("Gradient of a vector-valued function requires a seed")
    elseif seed == :auto
        zval = tape[z].val
        @assert zval isa Number || zval isa AbstractArray
        seed = zval isa Number ? one(zval) : ones(eltype(zval), size(zval))
    end
    dy = push!(tape, Constant(seed))
    # save seed var to use in compilation later
    tape.meta[:seed] = dy
    # set initial derivative value
    tape.c.derivs[z] = dy
    # queue of variables to calculate derivatives for
    deriv_todo = todo_list(tape)
    for y in deriv_todo
        try
            step_back!(tape, y)
        catch e
            BACKPROP_STATE[] = (tape, y, deriv_todo)
            rethrow(e)
        end
    end
end

function gradtape!(tape::Tape; seed=1)
    # backpropagate gradients
    back!(tape; seed=seed)
    # add a tuple of (val, (gradients...))
    deriv_vars = [hasderiv(tape, v) ? getderiv(tape, v) : ZeroTangent() for v in inputs(tape)]
    deriv_tuple = push!(tape, mkcall(tuple, deriv_vars...))
    # unthunk results
    deriv_tuple_unthunked = push!(tape, mkcall(map, ChainRules.unthunk, deriv_tuple))
    new_result = push!(tape, mkcall(tuple, tape.result, deriv_tuple_unthunked))
    # set result
    tape.result = new_result
    return tape
end

"""
    gradtape(f::Union{Function, DataType}, args...; ctx=GradCtx(), seed=1)
    gradtape!(tape::Tape; seed=1)

Calculate and record to the tape gradients of `tape[tape.resultid]` w.r.t. `Input` nodes.
See grad() for more high-level API.
"""
function gradtape(f::Union{Function,DataType}, args...; ctx=GradCtx(), seed=1)
    _, tape = trace(f, args...; ctx=ctx)
    tape = gradtape!(tape; seed=seed)
    return tape
end

"Like Umlaut.compile, but adds Yota specific ops"
function grad_compile(tape::Tape)
    ex = Umlaut.to_expr(tape)
    seed_var = tape.meta[:seed]
    seed_default_val = tape[seed_var].val
    insert!(ex.args[1].args, 2, Expr(:parameters, Expr(:kw, :seed, seed_default_val)))
    seed_var_name = Umlaut.make_name(seed_var.id)
    op_exprs = ex.args[2].args
    for op_ex in op_exprs
        if op_ex.args[1] == seed_var_name
            op_ex.args[2] = :seed
            break
        end
    end
    return Base.eval(@__MODULE__, ex)
end

const GRAD_CACHE = Dict{Any,Any}()

reset!() = empty!(GRAD_CACHE)

"""
    grad(f, args...; seed=1)

Find gradient of a callable `f` w.r.t. its arguments.

`grad()` returns two things: value of `f(args...)` and a tuple of
grafients w.r.t. to its inputs (including the callable itself).

```jldoctest
using Yota   # hide

val, g = grad(x -> sum(x .+ 1), [1.0, 2.0, 3.0])

# output
(9.0, (ChainRulesCore.ZeroTangent(), [1.0, 1.0, 1.0]))
```

By default, `grad()` expects the callable to return a scalar.
Vector-valued functions can be differentiated if a seed (starting value)
is provided. Seed is equivalent to the vector in VJP notation.

```jldoctest
using Yota   # hide

val, g = grad(x -> 2x, [1.0, 2.0, 3.0]; seed=ones(3))

# output
([2.0, 4.0, 6.0], (ChainRulesCore.ZeroTangent(), [2.0, 2.0, 2.0]))
```

All gradients can be applied to original variables using `update!()` function.

See also: [gradtape](@ref)
"""
function grad(f::Union{Function,DataType}, args...; seed=1)
    # key consists of function type and type of argument (for structs) or its size
    # cache_key = (f, ([(isstruct(arg) || args isa Function) ? typeof(arg) : size(arg) for arg in args]...,))
    cache_key = map(typeof, (f, args...))
    if haskey(GRAD_CACHE, cache_key)
        gf = GRAD_CACHE[cache_key]
        return gf(f, args...; seed=seed)
    else
        tape = gradtape(f, args...; seed=seed)
        gf = grad_compile(tape)
        GRAD_CACHE[cache_key] = gf
        return tape.retval
    end
end


text:
###########################
### Auxiliary Functions ###
###########################

function haslanes(w::OpenStreetMapX.Way)
    v = get(w.tags, "lanes", "")
    length(v)==1 && '1' <= v[1] <= '9'
end

getlanes(w::OpenStreetMapX.Way) = parse(Int, w.tags["lanes"])

visible(obj::T) where T <: OSMElement = (get(obj.tags, "visible", "") != "false")

services(w::OpenStreetMapX.Way) = (get(w.tags,"highway", "") == "services")

########################
### Extract Highways ###
########################

extract_highways(ways::Vector{OpenStreetMapX.Way}) = [way for way in ways if isdefined(way,:tags) && haskey(way.tags, "highway")]

filter_highways(ways::Vector{OpenStreetMapX.Way}) = [way for way in ways if OpenStreetMapX.visible(way) && !OpenStreetMapX.services(way)]

##############################################
### Filter and Classify Highways for Cars ###
##############################################

filter_roadways(ways::Vector{OpenStreetMapX.Way}, classes::Dict{String, Int} = OpenStreetMapX.ROAD_CLASSES; levels::Set{Int} = Set(1:length(OpenStreetMapX.ROAD_CLASSES))) = [way for way in ways if way.tags["highway"] in keys(classes) && classes[way.tags["highway"]] in levels]

classify_roadways(ways::Vector{OpenStreetMapX.Way}, classes::Dict{String, Int} = OpenStreetMapX.ROAD_CLASSES) = Dict{Int,Int}(way.id => classes[way.tags["highway"]] for  way in ways if haskey(classes, way.tags["highway"]))

####################################################
### Filter and Classify Highways for Pedestrians ###
####################################################

function filter_walkways(ways::Vector{OpenStreetMapX.Way},classes::Dict{String, Int} = OpenStreetMapX.PED_CLASSES; levels::Set{Int} = Set(1:length(OpenStreetMapX.PED_CLASSES)))
    walkways = OpenStreetMapX.Way[]
    for way in ways
        sidewalk = get(way.tags, "sidewalk", "")
        if sidewalk != "no"
            if haskey(classes, "sidewalk:$(sidewalk)") && classes["sidewalk:$(sidewalk)"] in levels
                push!(walkways,way)
            elseif haskey(classes, way.tags["highway"]) && classes[way.tags["highway"]] in levels
                push!(walkways,way)
            end
        end
    end
    return walkways
end

function classify_walkways(ways::Vector{OpenStreetMapX.Way},classes::Dict{String, Int} = OpenStreetMapX.PED_CLASSES)
    walkways = Dict{Int,Int}()
    for way in ways
        sidewalk = get(way.tags, "sidewalk", "")
        if sidewalk != "no"
            if haskey(classes, "sidewalk:$(sidewalk)")
                walkways[way.id] = classes["sidewalk:$(sidewalk)"]
            elseif haskey(classes, way.tags["highway"])
                walkways[way.id] = classes[way.tags["highway"]]
            end
        end
    end
    return walkways
end

###############################################
### Filter and Classify Highways for Cycles ###
###############################################

function filter_cycleways(ways::Vector{OpenStreetMapX.Way}, classes::Dict{String, Int} = OpenStreetMapX.CYCLE_CLASSES; levels::Set{Int} = Set(1:length(OpenStreetMapX.CYCLE_CLASSES)))
    cycleways = OpenStreetMapX.Way[]
    for way in ways
        bicycle = get(way.tags, "bicycle", "")
        cycleway = get(way.tags, "cycleway", "")
        highway = get(way.tags, "highway", "")

        cycleclass = "cycleway:$(cycleway)"
        bikeclass = "bicycle:$(bicycle)"

        if bicycle != "no"
            if haskey(classes, cycleclass) && classes[cycleclass] in levels
                push!(cycleways, way)
            elseif haskey(classes, bikeclass) && classes[bikeclass] in levels
                push!(cycleways, way)
            elseif haskey(classes, highway) && classes[highway] in levels
                push!(cycleways, way)
            end
        end
    end
    return cycleways
end

function classify_cycleways(ways::Vector{OpenStreetMapX.Way}, classes::Dict{String, Int} = OpenStreetMapX.CYCLE_CLASSES)
    cycleways = Dict{Int,Int}()
    for way in ways
        bicycle = get(way.tags, "bicycle", "")
        cycleway = get(way.tags, "cycleway", "")
        highway = get(way.tags, "highway", "")

        cycleclass = "cycleway:$(cycleway)"
        bikeclass = "bicycle:$(bicycle)"

        if bicycle != "no"
            if haskey(classes, cycleclass)
                cycleways[way.id] = classes[cycleclass]
            elseif haskey(classes, bikeclass)
                cycleways[way.id] = classes[bikeclass]
            elseif haskey(classes, highway)
                cycleways[way.id] = classes[highway]
            end
        end
    end
    return cycleways
end

##############################################
### Extract, Filter and Classify Buildings ###
##############################################

extract_buildings(ways::Vector{OpenStreetMapX.Way}) = [way for way in ways if isdefined(way,:tags) && haskey(way.tags, "building")]

filter_buildings(buildings::Vector{OpenStreetMapX.Way}, classes::Dict{String, Int} = OpenStreetMapX.BUILDING_CLASSES; levels::Set{Int} = Set(1:length(OpenStreetMapX.BUILDING_CLASSES))) = [building for building in buildings if (building.tags["building"] in keys(classes)) && (classes[building.tags["building"]] in levels) && OpenStreetMapX.visible(building)]

classify_buildings(buildings::Vector{OpenStreetMapX.Way}, classes::Dict{String, Int} = OpenStreetMapX.BUILDING_CLASSES) = Dict{Int,Int}(building.id => classes[building.tags["building"]] for  building in buildings if haskey(classes, building.tags["building"]))

#############################################
### Extract, Filter and Classify Features ###
#############################################

filter_features(features::Dict{Int,Tuple{String,String}}, classes::Dict{String, Int} = OpenStreetMapX.FEATURE_CLASSES; levels::Set{Int} = Set(1:length(OpenStreetMapX.FEATURE_CLASSES))) = Dict{Int,Tuple{String,String}}(key => feature for (key,feature) in features if classes[feature[1]] in levels)

function filter_features!(osmdata::OpenStreetMapX.OSMData, classes::Dict{String, Int} = OpenStreetMapX.FEATURE_CLASSES; levels::Set{Int} = Set(1:length(OpenStreetMapX.FEATURE_CLASSES)))
	osmdata.features = filter_features(osmdata.features, classes, levels = levels)
end

classify_features(features::Dict{Int,Tuple{String,String}}, classes::Dict{String, Int} = OpenStreetMapX.FEATURE_CLASSES) = Dict{Int,Int}(key =>  classes[feature[1]] for (key, feature) in features if haskey(classes, feature[1]))

### filter features in graph ###

function filter_graph_features(features::Dict{Int,Tuple{String,String}}, graphFeatures::Dict{Int,Int},classes::Dict{String,Int},class::String)
    if !haskey(classes,class)
        error("class not in classes")
    end
    level = classes[class]
    Dict{Int,Int}(key => node for (key,node) in graphFeatures if classes[features[key][1]] == level)
end


text:
using PyCall
@pyimport matplotlib.pyplot as plt
@pyimport matplotlib.lines as lines
export SeisPlot,
SeisPlotCoordinates
include("SeisPlot.jl")
include("SeisPlotCoordinates.jl")

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
"""
    AbstractDataArray{T, N}

An `N`-dimensional `AbstractArray` whose entries can take on values of type
`T` or the value `missing`.
"""
abstract type AbstractDataArray{T, N} <: AbstractArray{Union{T,Missing}, N} end

"""
    AbstractDataVector{T}

A 1-dimensional [`AbstractDataArray`](@ref) with element type `T`.
"""
const AbstractDataVector{T} = AbstractDataArray{T, 1}

"""
    AbstractDataMatrix{T}

A 2-dimensional [`AbstractDataArray`](@ref) with element type `T`.
"""
const AbstractDataMatrix{T} = AbstractDataArray{T, 2}

Base.eltype(d::AbstractDataArray{T, N}) where {T, N} = Union{T,Missing}

# Generic iteration over AbstractDataArray's

Base.start(x::AbstractDataArray) = 1
Base.next(x::AbstractDataArray, state::Integer) = (x[state], state + 1)
Base.done(x::AbstractDataArray, state::Integer) = state > length(x)

# FIXME: type piracy
"""
    ismissing(a::AbstractArray, i) -> Bool

Determine whether the element of `a` at index `i` is missing, i.e. `missing`.

# Examples

```jldoctest
julia> X = @data [1, 2, missing];

julia> ismissing(X, 2)
false

julia> ismissing(X, 3)
true
```
"""
Missings.ismissing(a::AbstractArray{T}, i::Real) where {T} = Missing <: T ? isa(a[i], Missing) : false # -> Bool


text:
using CSV, Dates, BenchmarkTools, InMemoryDatasets
using Revise
using AbnormalReturns

##
ds_firm = CSV.File(joinpath("data", "firm_ret.csv")) |> Dataset
ds_mkt = CSV.File(joinpath("data", "mkt_ret.csv")) |> Dataset
ds_events = CSV.File(joinpath("data", "event_dates.csv")) |> Dataset |> unique

##
@time ds_all = innerjoin(
    ds_firm,
    ds_mkt,
    on=[:date]
)
# 8.563892 seconds (4.58 k allocations: 7.422 GiB, 19.65% gc time, 0.05% compilation time)
##

@time ds_event_joined = innerjoin(
    ds_all,
    ds_events[:, [:firm_id, :event_date, :est_window_start, :est_window_end]],
    on=[:firm_id => :firm_id, :date => (:est_window_start, :est_window_end)]
)
# 52.283242 seconds (11.54 k allocations: 23.178 GiB, 7.14% gc time)
##
@time groupby!(ds_event_joined, [:firm_id, :event_date])
# 45.565012 seconds (1.01 M allocations: 26.786 GiB, 16.40% gc time)
##
function simple_reg(xs...)
    pred = disallowmissing(hcat(xs[2:end]...))
    resp = disallowmissing(xs[1])
    [cholesky!(Symmetric(pred' * pred)) \ (pred' * resp)]
end

@time InMemoryDatasets.combine(ds_event_joined, (:ret, :mkt, :hml, :umd, :smb) => simple_reg)
# 15.090089 seconds (35.68 M allocations: 19.679 GiB, 56.77% gc time, 8.81% compilation time)

##

@time ds_event_joined = innerjoin(
    ds_all,
    ds_events[:, [:firm_id, :event_date, :event_window_start, :event_window_end]],
    on=[:firm_id => :firm_id, :date => (:event_window_start, :event_window_end)]
)

##

@time groupby!(ds_event_joined, [:firm_id, :event_date])


text:
read_probabilities(str::AbstractString, T::Type=Int) = open(io -> read_probabilities(io, T), str)
read_probabilities(io::IOStream, T::Type=Int) = AliasTable(parse.(T, readlines(io)))

write_graph(str::AbstractString, edge_list::Vector{Vector{T}}) where T = open(io -> write_graph(io, edge_list), str, "w")
function write_graph(io::IOStream, edge_list::Vector{Vector{T}}) where T
    for edge in edge_list
        join(io, edge, ' ')
        println(io)
    end
end

max_stars_repo_name: LilithHafner/FBD

max_issues_repo_name: LilithHafner/FBD

max_forks_repo_name: LilithHafner/FBD


text:
using ChebyshevQuantum
using Test

tol_val = 1e-11

deriv=:forwarddiff

@testset "find edge function " begin

    #step function
    function f(x)
        if x < 0.2
            return 0.0
        else
            return 1.0
        end
    end

    c, stddev = ChebyshevQuantum.EdgeDetect.find_edges(f, deriv=deriv)

    @test isapprox(c[1,1], 0.2,  atol=1e-4)
    @test isapprox(c[1,2], 0.2,  atol=1e-4)

end

@testset "find edge f, f' " begin

    #two edges
    function f(x)
        if x < 0.3
            return 1.0 * abs(x-0.1) + x^3
        else
            return 1.0 * abs(x-0.1) + x^3 - 1
        end
    end
    #f(0.1)
    #f.([0.1,0.2])
    c, contin = ChebyshevQuantum.EdgeDetect.find(f, deriv=deriv)

    @test isapprox(c[1,1], 0.1,  atol=1e-4)
    @test isapprox(c[1,2], 0.1,  atol=1e-4)

    @test isapprox(c[2,1], 0.3,  atol=1e-4)
    @test isapprox(c[2,2], 0.3,  atol=1e-4)

    @test contin[1] == true
    @test contin[2] == false

end

@testset "find edge singularity " begin

    #singularity
    function f(x)
        return 2.0 * sqrt(abs(x-0.1))
    end

    c, contin = ChebyshevQuantum.EdgeDetect.find(f)

    @test isapprox(c[1,1], 0.1,  atol=1e-4)
    @test isapprox(c[1,2], 0.1,  atol=1e-4)

    @test contin[1] == true

    #singularity
    function f2(x)
        return -max(1.0 / (x-0.1), -1e10)
    end

    c2, contin2 = ChebyshevQuantum.EdgeDetect.find(f2)

    @test isapprox(c2[1,1], 0.1,  atol=1e-4)
    @test isapprox(c2[1,2], 0.1,  atol=1e-4)

    @test contin2[1] == false

end

@testset "don't find  " begin

    #no edges
    function f(x)
        return sin(30*x)
    end

    c, contin = ChebyshevQuantum.EdgeDetect.find(f)

    @test size(c)[1] == 0

    #no edges
    function f2(x)
        return exp(-1/x^2)
    end

    c2, contin2 = ChebyshevQuantum.EdgeDetect.find(f2)

    @test size(c2)[1] == 0

    #edge at end, don't find
    function f3(x)
        return max(-1/x, -1e10)
    end

    c3, contin3 = ChebyshevQuantum.EdgeDetect.find(f3, a=0, b =1)

    @test size(c3)[1] == 0

end


text:
module Displaz
using StaticArrays
using Colors

export plot3d, plot3d!, plotimage, plotimage!, clearplot, viewplot
export KeyEvent, CursorPosition, event_loop

"""
    set_displaz_cmd(cmd)

Set name or full path for where the displaz binary will be found to `cmd`.
Defaults to the environment variable `DISPLAZ_CMD`, or the string `"displaz"`
if that variable is not found.
"""
set_displaz_cmd(cmd) = global _displaz_cmd = cmd

function __init__()
    set_displaz_cmd(get(ENV, "DISPLAZ_CMD", "displaz"))
end

# Convert julia array into a type name and type appropriate for putting in the
# ply header
ply_type_convert(a::AbstractArray{UInt8})    = ("uint8",   a)
ply_type_convert(a::AbstractArray{UInt16})   = ("uint16",  a)
ply_type_convert(a::AbstractArray{UInt32})   = ("uint32",  a)
ply_type_convert(a::AbstractArray{Int8})     = ("int8",    a)
ply_type_convert(a::AbstractArray{Int16})    = ("int16",   a)
ply_type_convert(a::AbstractArray{Int32})    = ("int32",   a)
ply_type_convert(a::AbstractArray{Float32})  = ("float32", a)
ply_type_convert(a::AbstractArray{Float64})  = ("float64", a)
# Generic cases - actually do a conversion
ply_type_convert(a::AbstractArray{T}) where T <: Unsigned = ("uint32",  map(UInt32,a))
ply_type_convert(a::AbstractArray{T}) where T <: Integer = ("int32",   map(Int32,a))
ply_type_convert(a::AbstractArray{T}) where T <: Real = ("float64", map(Float64,a))

const array_semantic = 0
const vector_semantic = 1
const color_semantic = 2

# get ply header property name for given field index and semantic
function ply_property_name(semantic, idx)
    if semantic == array_semantic
        string(idx-1)
    elseif semantic == vector_semantic
        ("x", "y", "z", "w")[idx]
    elseif semantic == color_semantic
        ("r", "g", "b")[idx]
    end
end

# Write a set of points to displaz-native ply format
function write_ply_points(filename, nvertices, fields)
    converted_fields = [ply_type_convert(value) for (_,__,value) in fields]
    open(filename, "w") do fid
        write(fid, "ply\n")
        write(fid, "format binary_little_endian 1.0\n")
        write(fid, "comment Displaz native\n")
        for ((name,semantic,_), (typename, value)) in zip(fields, converted_fields)
            n = size(value,2)
            @assert(n == nvertices || n == 1)
            write(fid, "element vertex_$name $n\n")
            for i = 1:size(value,1)
                propname = ply_property_name(semantic, i)
                write(fid, "property $typename $propname\n")
            end
        end
        write(fid, "end_header\n")
        for (_,value) in converted_fields
            write(fid, value)
        end
    end
end

function write_ply_lines(filename, position, color, linebreak)
    nvalidvertices = size(position,2)

    # Create and write to ply file
    fid = open(filename, "w")
    write(fid,
        """
        ply
        format binary_little_endian 1.0
        element vertex $nvalidvertices
        property double x
        property double y
        property double z
        element color $nvalidvertices
        property float r
        property float g
        property float b
        element edge $(length(linebreak))
        property list int int vertex_index
        end_header
        """
    )

    write(fid,convert(Array{Float64,2},position))
    write(fid,color)

    realstart = 0
    linelen = []
    range = []
    for i = 1:size(linebreak,1)
        if i != size(linebreak,1)
            linelen = linebreak[i+1] - linebreak[i]
            range = realstart:realstart + linelen-1
        else
            linelen = size(position,2) - linebreak[i] + 1
            range = realstart:realstart + linelen - 1
        end
        write(fid,Int32(linelen))
        write(fid,UnitRange{Int32}(range))
        realstart = realstart + linelen
    end
    close(fid)
end

#const standard_elements = [:position  => (vector_semantic,3),
#                           :color     => (color_semantic,3),
#                           :marksize  => (array_semantic,1),
#                           :markshape => (array_semantic,1)]

const _color_names = Dict('r' => [1.0, 0,   0],
                          'g' => [0.0, 0.8, 0],
                          'b' => [0.0, 0,   0.8],
                          'c' => [0.0, 1,   1],
                          'm' => [1.0, 0,   1],
                          'y' => [1.0, 1,   0],
                          'k' => [0.0, 0,   0],
                          'w' => [1.0, 1,   1])

const _shape_ids = Dict('.' => 0,
                        's' => 1,
                        'o' => 2,
                        'x' => 3,
                        '+' => 4)

interpret_color(color) = color
interpret_color(s::AbstractString) = length(s) == 1 ? interpret_color(s[1]) : error("Unknown color abbreviation $s")
interpret_color(c::Char) = _color_names[c]
interpret_color(c::AbstractRGB) = [red(c), green(c), blue(c)]
interpret_color(c::Color) = interpret_color(RGB(c))
function interpret_color(cs::AbstractVector{T}) where T <: Color
    a = zeros(eltype(T), (3,length(cs)))
    for i=1:length(cs)
        c = RGB(cs[i])
        a[1,i] = red(c)
        a[2,i] = green(c)
        a[3,i] = blue(c)
    end
    a
end

interpret_shape(markershape) = markershape
interpret_shape(c::Char) = [_shape_ids[c]]
interpret_shape(s::Vector{Char}) = Int[_shape_ids[c] for c in s]
interpret_shape(s::AbstractString) = Int[_shape_ids[c] for c in s]

interpret_linebreak(nvertices, linebreak) = linebreak
interpret_linebreak(nvertices, i::Integer) = i == 1 ? [1] : 1:i:nvertices

interpret_position(pos::AbstractMatrix) = pos
function interpret_position(pos::AbstractVector{V}) where V <: StaticVector
    size(eltype(pos)) == (3,) || error("position should be a 3-vector")
    T = eltype(V)
    isbitstype(T) || error("Can't reinterpret position with elements $T")
    nvertices = length(pos)
    return reshape(reinterpret(T, pos), (3, nvertices))
end

function interpret_position(pos::AbstractVector{T}) where T <: Real
    size(pos) == (3,) || error("position should be a 3-vector")
    return pos
end

# Multiple figure window support
# TODO: Consider how the API below relates to Plots.jl and its tendency to
# create a lot of new figure windows rather than clearing existing ones.
mutable struct DisplazWindow
    name::AbstractString
end

_current_figure = DisplazWindow("default")
"Get handle to current figure window"
function current()
    _current_figure
end

"Get figure window by name may be new or already existing"
function figure(name::AbstractString)
    global _current_figure
    _current_figure = DisplazWindow(name)
end
"Get figure window with given id"
figure(id::Integer) = figure("Figure $id")

_figure_id = 1
"Create next incrementally named figure window, counting automatically from \"Figure 1\""
function newfigure()
    global _figure_id
    id = _figure_id
    _figure_id += 1
    figure(id)
end

"""
Add 3D points or lines to the current plot.

```
  plot3d([plotobj,] position; attr1=value1, ...)
```

The `position` array should be a set of N vertex positions, specified as 3xN
array or a `Vector` of `FixedVector{3}`.  The `plotobj` argument is optional
and determines which plot window to send the data to.  If it's not used the
data will be sent to the plot window returned by `current()`.

TODO: Revisit the nasty decision of the shape of position again - the above
choice is somewhat inconsistent with supplying markersize / markershape as a
column vector :-(  Can we have a set of consistent broadcasting rules for this?
It seems like the case of a 3x3 matrix will always be ambiguous if we try
to guess what the user wants.

### Data set attributes

The following attributes can be attached to a dataset on each call to `plot3d`:

  * `label` - A string labeling the data set

### Vertex attributes

Each point may have a set of vertex attributes attached to control the visual
representation and tag the point for inspection. You can pass any `Vector` of
`Float32` values for any custom information you like, but the following are
supported by the default shader:

  * `color`       - A color or vector of colors for each point; see below for
                    ways to specify these.
  * `intensity`   - A vector of the intensity of each point (between 0 and 1)
  * `markersize`  - Vertex marker size
  * `markershape` - Vector of vertex marker shapes.  Shape can be represented
                    either by a Char or a numeric identifier between 0 and 4:

```
                    sphere - '.' 0    square - 's' 1
                    circle - 'o' 2    times  - 'x' 3
                    cross  - '+' 4
```

#### Specifying colors

Colors may be provided in any of three ways:

* As instances of types from the ColorTypes package, for example, HSV(180,1,1).
  These are converted to RGB using the RGB constructor.
* As a `Vector` of three elements, red, green and blue, between 0.0 and 1.0.
* Using a matlab-like single color letter name string or `Char`.  Supported are
  red, green, blue, cyan, magenta, yellow, black and white; all are abbreviated
  with the first letter of the color name except black for which 'k' is used.

A color per point may be supplied as a `Vector` of `Color` subtypes or a 3xN
matrix with red, green and blue in the rows.

### Plotting points

To plot 10000 random points with distance from the origin determining the
color, and random marker shapes:
```
  P = randn(3,10000)
  c = 0.5./sumabs2(P,1) .* [1,0,0]
  plot3d(P, color=c, markershape=rand(1:4,10000))
```

### Plotting lines

To plot a piecewise linear curve between 10000 random vertices
```
  plot3d(randn(3,10000), markershape="-")
```

When plotting lines, the `linebreak` keyword argument can be used to break the
position array into multiple line segments.  Each index in the line break array
is the initial index of a line segment.
"""
function plot3d(plotobj::DisplazWindow, position; color=[1,1,1], markersize=[0.1], markershape=[0],
                label=nothing, linebreak=[1], _overwrite_label=false, shader="generic_points.glsl", kwargs...)
    position = interpret_position(position)
    nvertices = size(position, 2)
    color = interpret_color(color)
    linebreak = interpret_linebreak(nvertices, linebreak)
    size(position, 1) == 3 || error("position must be a 3xN array")
    size(color, 1)    == 3 || error("color must be a 3xN array")
    # FIXME in displaz itself.  No repeat waste should be required.
    if size(color,2) == 1
        color = repeat(color, 1, nvertices)
    end

    # works for vectors only at this stage...
    extra_fields = ((x -> (size(x[2]) == (nvertices,) || error("extra fields must be vectors of same length as number of points in position array") ;
        (x[1], array_semantic, transpose(Float32.(x[2])))) for x = pairs(kwargs)))

    # Ensure all fields are floats for now, to avoid surprising scaling in the
    # shader
    color = map(Float32,color)
    markersize = map(Float32,markersize)
    size(color,2) == nvertices || error("color must have same number of rows as position array")
    filename = tempname()*".ply"
    seriestype = "Points"
    if markershape == "-" || markershape == '-'
        # Plot lines
        # FIXME: The way this is detected is a bit of a mess - lines vs points
        # should be plotted using separate functions.
        seriestype = "Line"
        write_ply_lines(filename, position, color, linebreak)
    else # Plot points
        if length(markersize) == 1
            markersize = repeat(markersize, 1, nvertices)
        end
        markershape = interpret_shape(markershape)
        if length(markershape) == 1
            markershape = repeat(markershape, 1, nvertices)
        end
        write_ply_points(filename, nvertices, (
                         (:position, vector_semantic, position),
                         (:color, color_semantic, color),
                         # FIXME: shape of markersize??
                         (:markersize, array_semantic, vec(markersize)'),
                         (:markershape, array_semantic, vec(markershape)'),
                         extra_fields...
                         ))
    end
    if label === nothing
        label = "$seriestype [$nvertices vertices]"
    end
    addopt = _overwrite_label ? [] : "-add"
    run(`$_displaz_cmd -script $addopt -server $(plotobj.name) -label $label -shader $shader -rmtemp $filename`)
    nothing
end

"""
    mutate!([plotobj,] label, index; attr1=value1, ...))

Mutate the data in an existing displaz data set `label`, for instance to change the
position or other attribute of a subset of points (with the advantage of
reducing the amount of communication between Julia and displaz, and therefore
increasing speed).

`index` is vector of indices with reference to the original plot. The attribute
`label` is used to match with the correct data set within displaz. The
`position` attribute controls the vertex positions, and the remainder match
the original plotting command.
"""
function mutate!(label::AbstractString, index::AbstractVector{I}; kwargs...) where I <: Integer
    plotobj = _current_figure
    mutate!(plotobj, label, index; kwargs...)
end

function mutate!(plotobj::DisplazWindow, label::AbstractString, index::AbstractVector{I}; kwargs...) where I <: Integer
    nvertices = length(index)

    fields = Vector{Any}()

    push!(fields, (:index, array_semantic, map(UInt32, (index.-1)'))) # It turns out the -1 is rather important (0-based indexing)... :)

    for (fieldname, fielddata) ∈ kwargs
        if fieldname == :position
            fielddaata = interpret_position(fielddata)
            size(fielddata) == (3,nvertices) || error("position must be a 3x$nvertices array")

            push!(fields, (:position, vector_semantic, fielddata))
        elseif fieldname == :color
            fielddata = interpret_color(fielddata)
            size(fielddata, 1) == 3 || error("color must be a 3xN array")
            if size(fielddata,2) == 1
                fielddata = repeat(fielddata, 1, nvertices)
            end
            size(fielddata) == (3,nvertices,) || error("wrong number of color points")
            fielddata = map(Float32, fielddata)

            push!(fields, (:color, color_semantic, fielddata))
        elseif fieldname == :markershape
            if length(fielddata) == 1
                fielddata = repeat(fielddata, 1, nvertices)
            end
            size(fielddata) == (nvertices,) || error("wrong number of markershape points")
            fielddata = interpret_shape(fielddata)

            push!(fields, (:markershape, array_semantic, vec(fielddata)'))
        elseif fieldname == :linebreak
            if length(fielddata) == 1
                fielddata = repeat(fielddata, 1, nvertices)
            end
            fielddata = interpret_linebreak(fielddata)

            push!(fields, (:linebreak, array_semantic, vec(fielddata)'))
        else
            if length(fielddata) == 1
                fielddata = repeat(fielddata, nvertices)
            end
            size(fielddata) == (nvertices,) || error("extra fields must be vectors of same length as index array")
            fielddata = map(Float32, fielddata)

            push!(fields, (fieldname, array_semantic, vec(fielddata)'))
        end
    end

    filename = tempname()*".ply"

    write_ply_points(filename, nvertices, fields)
    run(`$_displaz_cmd -script -modify -server $(plotobj.name) -label $label -shader generic_points.glsl $filename`)
    nothing

end

"""
Overwrite points or lines with the same label on the 3D plot

See plot3d for documentation
"""
function plot3d!(plotobj::DisplazWindow, position; kwargs...)
    plot3d(plotobj, position; _overwrite_label=true, kwargs...)
end

# Plot to current window
plot3d!(position; kwargs...) = plot3d!(current(), position; kwargs...)
plot3d(position; kwargs...)  = plot3d(current(), position; kwargs...)

#-------------------------------------------------------------------------------
"""
    clearplot([plotobj::DisplazWindow=current()], [pattern1, ...])

Clear all or a subset of datasets in a plot window.

If not specified, `plotobj` is the current plot window.  If no patterns are
supplied, clears all data sets.  If one or more patterns are given, the dataset
labels matching those patterns will be removed.  Patterns shoudl be specified
using unix shell glob pattern syntax.
"""
function clearplot(plotobj::DisplazWindow, patterns...)
    unload_args = isempty(patterns) ? ["-clear"] : ["-unload", patterns...]
    run(`$_displaz_cmd -script -server $(plotobj.name) $unload_args`)
    nothing
end
clearplot(patterns...) = clearplot(current(), patterns...)

#-------------------------------------------------------------------------------
# Texture file

# Write image name (including absolute path) and vertices of the four corner to ply file
function write_ply_texture(texturefile::String, filename::String, vertices::AbstractArray)

    open(filename, "w") do f
        write(f,
            """
            ply
            format ascii 1.0
            comment TextureFile $texturefile
            element vertex 4
            property double x
            property double y
            property double z
            property double u
            property double v
            element face 1
            property list int int vertex_index
            end_header
            $(vertices[1,1]) $(vertices[1,2]) $(vertices[1,3]) 0 1
            $(vertices[2,1]) $(vertices[2,2]) $(vertices[2,3]) 1 1
            $(vertices[3,1]) $(vertices[3,2]) $(vertices[3,3]) 1 0
            $(vertices[4,1]) $(vertices[4,2]) $(vertices[4,3]) 0 0
            4 0 1 2 3
            """
        )
    end
    nothing
end

"""
Add images to the current plot.

```
  plotimage([plotobj,] texturefile, vertices; label=nothing, _overwrite_label=false)
```

The `texturefile` string should include the path (relative or absolute) to the
image to be loaded (if not in the current folder).
The `vertices` array should be a set of vertex positions specifying the corners
of the image to plot. The order is anticlockwise starting from bottom left (i.e.
(0,1), (1,1), (1,0), (0,0) in texture coordinates), and the vertices
should be specified as a 3x4 array.
The `plotobj` argument is optional and determines which plot window
to send the data to.  If it's not used the data will be sent to the plot window
returned by `current()`.
"""
function plotimage(plotobj::DisplazWindow, texturefile::String, vertices::AbstractArray; label=nothing, _overwrite_label=false)

    filename = tempname()*".ply"

    write_ply_texture(abspath(texturefile), filename, vertices)

    if label === nothing
        label = "$texturefile"
    end
    addopt = _overwrite_label ? [] : "-add"
    run(`$_displaz_cmd -script $addopt -server $(plotobj.name) -label $label -rmtemp $filename`)
    nothing
end

"""
Overwrite images with the same label

See plotimage for documentation
"""
function plotimage!(plotobj::DisplazWindow, texturefile::String, vertices::AbstractArray; kwargs...)
    plotimage(plotobj, texturefile, vertices; _overwrite_label=true, kwargs...)
end

# Plot to current window
plotimage!(texturefile, vertices; kwargs...) = plotimage!(current(), texturefile, vertices; kwargs...)
plotimage(texturefile, vertices; kwargs...)  = plotimage(current(), texturefile, vertices; kwargs...)

#-------------------------------------------------------------------------------
# 3D camera control
"""
    viewplot([plotobj::DisplazWindow=current()], label)

Set the point of view of the 3D camera.  The camera model is designed to view
an object at a given `center` of rotation, with rotations pivoting around that
position.

If not specified, `plotobj` is the current plot window.  Keyword arguments
define the camera view as follows:

* The `center` argument may be a 3D point or the label of a data set.
  If a label is supplied the centroid of the associated dataset will be used.
* The `radius` argument should be a number giving the distance that the
  camera will be away from the center of rotation.
* The `rotation` argument specifies the angles at which the camera will view
  the scene.  This should be a matrix transforming points into the standard
  OpenGL camera coordinates (+x right, +y up, -z into the scene).

"""
function viewplot(plotobj::DisplazWindow;
                  center=nothing, radius=nothing, rotation=nothing)
    center_args   = viewplot_center_args(center)
    radius_args   = viewplot_radius_args(radius)
    rotation_args = viewplot_rotation_args(rotation)
    run(`$_displaz_cmd -script -server $(plotobj.name) $center_args $radius_args $rotation_args`)
    nothing
end
viewplot(; kwargs...) = viewplot(current(); kwargs...)

# viewplot() helper stuff
# center
viewplot_center_args(::Nothing) = []
viewplot_center_args(s::AbstractString) = ["-viewlabel", string(s)]
viewplot_center_args(pos) = ["-viewposition", string(pos[1]), string(pos[2]), string(pos[3])]
# rotation
viewplot_rotation_args(::Nothing) = []
viewplot_rotation_args(M) = vcat("-viewrotation", map(string, vec(Matrix(M)'))) # Generate row-major order
# radius
viewplot_radius_args(::Nothing) = []
viewplot_radius_args(r) = ["-viewradius", string(r)]

#-------------------------------------------------------------------------------
# Event loop

include("events.jl")

end


text:
function rename!(g::AbstractExGraph, name::Symbol, new_name::Symbol)
    return rename!(g, Dict(name => new_name))
end

function rename!(g::AbstractExGraph, st::Dict{Symbol,Symbol})
    for nd in g.tape
        var = getvar(nd)
        new_var = subs(var, st)
        setvar!(nd, subs(var, st))
        setexpr!(nd, subs(getexpr(nd), st))
        delete!(g.idx, var)
        g.idx[new_var] = nd
    end
    return g
end

function rename(g::AbstractExGraph, st::Dict{Symbol,Symbol})
    new_g = reset_tape(g)
    for nd in g.tape
        push!(new_g, copy(nd; var=subs(getvar(nd), st), ex=subs(getexpr(nd), st)))
    end
    return new_g
end

function mergeex(g1::AbstractExGraph, g2::AbstractExGraph)
    g2 = deepcopy(g2)
    @assert typeof(g1) == typeof(g2)
    gm = typeof(g1)()
    # find out what vars in g2 need to be renamed
    g1_vars = Set(varname(nd) for nd in g1.tape)
    g2_vars = Set(varname(nd) for nd in g2.tape)
    need_renaming = Symbol[]
    for nd in g2.tape
        vname = varname(nd)
        if haskey(g1, vname) && getexpr(g1[vname]) != getexpr(nd)
            push!(need_renaming, vname)
        end
    end
    # rename variables in g2
    new_names = gennames(length(need_renaming))
    for (name, new_name) in zip(need_renaming, new_names)
        rename!(g2, name, new_name)
    end
    # concat graphs
    for nd in g1.tape
        push!(gm, nd)
    end
    for nd in g2.tape
        if !haskey(gm, varname(nd))
            push!(gm, nd)
        end
    end
    return gm
end

function mergeex(gs::Vararg{AbstractExGraph})
    return reduce(mergeex, gs)
end

function mergeex(exs::Vararg{Expr})
    indexed = any(isindexed, exs)
    gs = indexed ? [EinGraph(ex) for ex in exs] : [ExGraph(ex) for ex in exs]
    gm = mergeex(gs...)
    block = to_expr(gm)
    res_vars = [g[end].var for g in gs]
    push!(block.args, Expr(:tuple, res_vars...))
    return block
end