text:
# # Vehicle Platoon Model
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/Platoon.ipynb)
#
#md # !!! note "Overview"
#md #     System type: affine with hybrid dynamics\
#md #     State dimension: 9 + 1\
#md #     Application domain: Autonomous Driving
#
# ## Model description
#
# This benchmark considers a platoon of three vehicles following each
# other. This benchmark considers loss of communication between vehicles. The model
# has two modes, connected and disconnected.
# Three scenarios are considered for the loss of communication:
#
# **PLAA01 (arbitrary loss)** The loss of communication can occur at any time. This
#       includes the possibility of no communication at all.
#
# **PLADxy (loss at deterministic times)** The loss of communication occurs at fixed
#       points in time, which are determined by clock constraints ``c_1`` and ``c_2``.
#       The clock t is reset when communication is lost and when it is re-established.
#       Note that the transitions have must-semantics, i.e., they take place as soon
#       as possible. We will consider `PLAD01: c1 = c_2 = 5`.
#
# **PLANxy (loss at nondeterministic times)** The loss of communication occurs at
#       any time ``t ∈ [t_b, t_c]``. The clock t is reset when communication is lost
#       and when it is reestablished. Communication is reestablished at any time
#       ``t ∈ [0, t_r]``. This scenario covers loss of communication after an
#       arbitrarily long time ``t ≥ t_c`` by reestablishing communication in zero time.
#       We will consider `PLAN01: tb = 10, tc = 20, tr = 20`.

# The Julia model is developed next. It's convenient to create two independent
# functions, `platoon_connected` and `platoon_disconnected`, which describe the
# dynamics of the connected (resp. disconnected) modes. These functions are
# connected in the `platoon` function using a hybrid automaton with two modes
# and two discrete transitions, see the following diagram.

# ![](https://github.com/JuliaReach/ReachabilityAnalysis.jl/blob/master/examples/Platoon/platoon_hybrid.png?raw=true)

# On the left, the system can switch arbitrarily between the modes.
# In the middle, mode switches are only possible at given points in time, i.e. deterministic switching.
# On the right, mode switches are only possible during given time intervals.
# In this notebook we will only consider the case of deterministic switching.

# ### Dynamics of the "connected" platoon

using ReachabilityAnalysis, SparseArrays, Symbolics

const var = @variables x[1:9] t

function platoon_connected(; deterministic_switching::Bool=true,
                             c1=5.0)  # clock constraints
    n = 9 + 1

    ## x' = Ax + Bu + c
    A = Matrix{Float64}(undef, n, n)
    A[1, :] = [0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0]
    A[2, :] = [0, 0, -1.0, 0, 0, 0, 0, 0, 0, 0]
    A[3, :] = [1.6050, 4.8680, -3.5754, -0.8198, 0.4270, -0.0450, -0.1942,  0.3626, -0.0946, 0.]
    A[4, :] = [0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0,]
    A[5, :] = [0, 0, 1.0, 0, 0, -1.0, 0, 0, 0, 0]
    A[6, :] = [0.8718, 3.8140, -0.0754,  1.1936, 3.6258, -3.2396, -0.5950,  0.1294, -0.0796, 0.]
    A[7, :] = [0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0]
    A[8, :] = [0, 0, 0, 0, 0, 1.0, 0, 0, -1.0, 0]
    A[9, :] = [0.7132, 3.5730, -0.0964,  0.8472, 3.2568, -0.0876,  1.2726,  3.0720, -3.1356, 0.]
    A[10, :] = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0]; # t' = 1

    if deterministic_switching
        invariant = HalfSpace(t <= c1, var)
    else
        invariant = Universe(n)
    end

    ## acceleration of the lead vehicle + time
    B = sparse([2], [1], [1.0], n, 1)
    U = Hyperrectangle(low=[-9.], high=[1.])
    c = [0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0]
    @system(x' = Ax + Bu + c, x ∈ invariant, u ∈ U)
end

# ### Dynamics of the "disconnected" platoon

function platoon_disconnected(; deterministic_switching::Bool=true,
                                c2=5.0)  # clock constraints
    n = 10 # 9 dimensions + time

    ## x' = Ax + Bu + c
    A = Matrix{Float64}(undef, n, n)
    A[1, :] = [0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0]
    A[2, :] = [0, 0, -1.0, 0, 0, 0, 0, 0, 0, 0]
    A[3, :] = [1.6050, 4.8680, -3.5754, 0, 0, 0, 0, 0, 0, 0]
    A[4, :] = [0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0,]
    A[5, :] = [0, 0, 1.0, 0, 0, -1.0, 0, 0, 0, 0]
    A[6, :] = [0, 0, 0,  1.1936, 3.6258, -3.2396, 0, 0, 0, 0.]
    A[7, :] = [0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0]
    A[8, :] = [0, 0, 0, 0, 0, 1.0, 0, 0, -1.0, 0]
    A[9, :] = [0.7132, 3.5730, -0.0964,  0.8472, 3.2568, -0.0876,  1.2726,  3.0720, -3.1356, 0.]
    A[10, :] = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0]; # t' = 1

    if deterministic_switching
        invariant = HalfSpace(t <= c2, var)
    else
        invariant = Universe(n)
    end

    ## acceleration of the lead vehicle + time
    B = sparse([2], [1], [1.0], n, 1)
    U = Hyperrectangle(low=[-9.], high=[1.])
    c = [0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0]
    @system(x' = Ax + Bu + c, x ∈ invariant, u ∈ U)
end

# ### Hybrid system

function platoon(; deterministic_switching::Bool=true,
                   c1=5.0,  # clock constraints
                   c2=5.0,  # clock constraints
                   tb=10.0,  # lower bound for loss of communication
                   tc=20.0, tr=20.0) # upper bound for loss of communication (tc) and reset time (tr)

    ## three variables for each vehicle, (ei, d(et)/dt, ai) for
    ## (spacing error, relative velocity, speed), and the last dimension is time
    n = 9 + 1

    ## transition graph
    automaton = LightAutomaton(2)
    add_transition!(automaton, 1, 2, 1)
    add_transition!(automaton, 2, 1, 2)

    ## modes
    mode1 = platoon_connected(deterministic_switching=deterministic_switching, c1=c1)
    mode2 = platoon_disconnected(deterministic_switching=deterministic_switching, c2=c2)
    modes = [mode1, mode2]

    ## common reset
    reset = Dict(n => 0.)

    ## transition l1 -> l2
    if deterministic_switching
        guard = Hyperplane(t == c1, var)
    else
        guard = HPolyhedron([tb <= t, t <= tc], var)
    end
    t1 = ConstrainedResetMap(n, guard, reset)

    ## transition l2 -> l1
    if deterministic_switching
        guard = Hyperplane(t == c2, var)
    else
        guard = HalfSpace(t <= tr, var)
    end
    t2 = ConstrainedResetMap(n, guard, reset)
    resetmaps = [t1, t2]

    H = HybridSystem(automaton, modes, resetmaps, [AutonomousSwitching()])

    ## initial condition is at the orgin in mode 1
    X0 = BallInf(zeros(n), 0.0)
    initial_condition = [(1, X0)]

    return IVP(H, initial_condition)
end

# ## Safety specifications
#
# The verification goal is to check whether the minimum distance between
# vehicles is preserved. The choice of the coordinate system is such that the
# minimum distance is a negative value.
#
# **BNDxy**: Bounded time (no explicit bound on the number of transitions): For all ``t ∈ [0, 20] [s]``,
#
#    - ``x_1(t) ≥ −d_{min} [m]``
#    - ``x_4(t) ≥ −d_{min} [m]``
#    - ``x_7(t) ≥ −d_{min} [m]``
#
# **BND42**: ``d_{min} = 42``.
#
# **BND30**: ``d_{min} = 30``.

function dmin_specification(sol, dmin)
    return (-ρ(sparsevec([1], [-1.0], 10), sol) > -dmin) &&
           (-ρ(sparsevec([4], [-1.0], 10), sol) > -dmin) &&
           (-ρ(sparsevec([7], [-1.0], 10), sol) > -dmin)
end

# ## Results

# We will only consider the case of deterministic switching in this notebook. We consider two
# scenarios of increasing difficulty, BND42 and BND30.

prob_PLAD01 = platoon();

const boxdirs = BoxDirections(10)
const octdirs = OctDirections(10);

#  ### PLAD01 - BND42

# This scenario can be solved using a hyperrectangular set representation with step
# size ``δ = 0.01``.  We use a template that contains all box (i.e., canonical)
# directions in the ambient space of the state-space, ``\mathbb{R}^{10}``.
# There are ``20`` such directions, two for each coordinate:

length(boxdirs)

#-

alg = BOX(δ=0.01)
sol_PLAD01_BND42 = solve(prob_PLAD01,
                         alg=alg,
                         clustering_method=BoxClustering(1),
                         intersection_method=TemplateHullIntersection(boxdirs),
                         intersect_source_invariant=false,
                         tspan = (0.0 .. 20.0));

# Let's verify that the specification holds:

dmin_specification(sol_PLAD01_BND42, 42)

# In more detail we can check how is the flowpipe from violating the property.
# The specification requires that each of the follwing quantities is greater
# than `-dmin = -42`.

# Minimum of ``x_1(t)``:
-ρ(sparsevec([1], [-1.0], 10), sol_PLAD01_BND42)
#-
# Minimum of ``x_4(t)``:
-ρ(sparsevec([4], [-1.0], 10), sol_PLAD01_BND42)
#-
# Minimum of ``x_7(t)``:
-ρ(sparsevec([7], [-1.0], 10), sol_PLAD01_BND42)

# We plot variable ``x_1`` vs time.

using Plots, LaTeXStrings

plot(sol_PLAD01_BND42, vars=(0, 1), xlab=L"t", ylab=L"x_1", title="PLAD01 - BND42", lw=0.1)
plot!(x->x, x->-42., 0., 20., linewidth=2, color="red", linestyle=:dash, leg=nothing)

#  ### PLAD01 - BND30

# Note that the previous solution obtained for `PLAD01 - BND42` does not verify the
# `BND30` specifications since, for example,
# the minimum of variable ``x_4(t)`` is ``≈ -35.52``, which
# is below the given bound ``-d_{min} = -30``. As a consequence, to prove the safety
# properties for this scenario we have to use a solver with more precision.
# Instead of the `BOX` algorithm, we propose to use `LGG09` with step-size
# ``δ=0.03`` and octagonal template directions. There are 200 such directions:

length(octdirs)

#md # !!! tip "Performance tip"
#md # The increase in the number of directions implies an increase in runtime.
#md # Since evaluating the 200 directions of the template is quite expensive,
#md # we try using a concrete set after the discretization, instead of using a
#md # lazy discretization. This is achieved by passing the option
#md # `approx_model=Forward(setops=octdirs)` to the `LGG09` algorithm, specifying that
#md # we want to oveapproximate the initial set of the set-based recurrence, that we call `Ω0`,
#md # with an octagonal template. It turns out in this example that the option gives a
#md # gain in runtime of ``~30\%``, without a noticeable loss in precision.

alg = LGG09(δ=0.03, template=octdirs, approx_model=Forward(setops=octdirs));
sol_PLAD01_BND30 = solve(prob_PLAD01,
                         alg=alg,
                         clustering_method=LazyClustering(1),
                         intersection_method=TemplateHullIntersection(octdirs),
                         intersect_source_invariant=false,
                         tspan = (0.0 .. 20.0));

# Verifying that the specification holds:

dmin_specification(sol_PLAD01_BND30, 30)

# Let's check in more detail how close is the flowpipe to the safety conditions:

# Minimum of ``x_1(t)``:
-ρ(sparsevec([1], [-1.0], 10), sol_PLAD01_BND30)
#-
# Minimum of ``x_4(t)``:
-ρ(sparsevec([4], [-1.0], 10), sol_PLAD01_BND30)
#-
# Minimum of ``x_7(t)``:
-ρ(sparsevec([7], [-1.0], 10), sol_PLAD01_BND30)

#-

plot(sol_PLAD01_BND30, vars=(0, 1), xlab=L"t", ylab=L"x_1", title="PLAD01 - BND30", lw=0.1)
plot!(x->x, x->-30., 0., 20., linewidth=2, color="red", linestyle=:dash, leg=nothing)

# ## References

# [^]:


text:
struct SymEGRSSCholesky{T,UT<:AbstractMatrix,WT<:AbstractMatrix} <: AbstractMatrix{T}
    Ut::UT
    Wt::WT
    n::Int
    p::Int
    function SymEGRSSCholesky{T,UT,WT}(Ut,Wt,n,p) where
			 {T,UT<:AbstractMatrix,WT<:AbstractMatrix}
        Un, Um = size(Ut)
		Wn, Wm = size(Wt)
		(Un == Wn && Um == Wm) || throw(DimensionMismatch())
        new(Ut,Wt,n,p)
    end
end

# Creating generator such that L = tril(U*W')
function ss_create_w(Ut,Vt)
    p,n = size(Ut);
    wTw = zeros(p,p);
    Wt = similar(Vt);
    @inbounds for j = 1:n
        tmpu = Ut[:,j];
        tmp  = Vt[:,j] - wTw*tmpu;
        w = tmp/sqrt(abs(tmpu'*tmp))
        Wt[:,j] = w;
        wTw = wTw + w*w';
    end
    return Wt
end

cholesky(K::SymEGRSSMatrix{T,UT,VT}) where {T,UT<:AbstractMatrix,VT<:AbstractMatrix} =
        SymEGRSSCholesky{T,typeof(K.Ut),typeof(K.Vt)}(K.Ut,ss_create_w(K.Ut,K.Vt),K.n,K.p)

########################################################################
#### Helpful properties. Not nessecarily computionally efficient    ####
########################################################################
Matrix(K::SymEGRSSCholesky) = getproperty(K,:L)

size(K::SymEGRSSCholesky) = (K.n, K.n)
size(K::SymEGRSSCholesky, d::Int) = (1 <= d && d <=2) ? size(K)[d] : throw(ArgumentError("Invalid dimension $d"))

function getindex(K::SymEGRSSCholesky, i::Int, j::Int)
	i >= j && return dot(K.Ut[:,i], K.Wt[:,j])
	return 0
end

function getproperty(K::SymEGRSSCholesky, d::Symbol)
    if d === :U
        return UpperTriangular(K.Wt'*K.Ut)
    elseif d === :L
        return LowerTriangular(K.Ut'*K.Wt)
    else
        return getfield(K, d)
    end
end

Base.propertynames(F::SymEGRSSCholesky, private::Bool=false) =
    (:U, :L, (private ? fieldnames(typeof(F)) : ())...)

function Base.show(io::IO, mime::MIME{Symbol("text/plain")},
		 K::SymEGRSSCholesky{<:Any,<:AbstractMatrix,<:AbstractMatrix})
    summary(io, K); println(io)
    show(io, mime, K.L)
end

########################################################################
#### Solving systems using the Cholesky factorization    			####
########################################################################

#### Backward substitution (solve Lx = b) ####
function ss_forward!(X::AbstractArray, Ut::AbstractArray,
                     Wt::AbstractArray, B::AbstractArray)
    p, n = size(Ut)
    mx = size(B,2)
    Wbar = zeros(p, mx);
    @inbounds for i = 1:n
        tmpU = Ut[:,i]
        tmpW = Wt[:,i]
        X[i,:] = (B[i:i,:] - tmpU'*Wbar)./(tmpU'*tmpW);
        Wbar += tmpW .* X[i:i,:];
    end
end

#### Backward substitution (solve L'x = b) ####
function ssa_backward!(X::AbstractArray, Ut::AbstractArray,
                       Wt::AbstractArray, B::AbstractArray)
    p, n = size(Ut)
    mx = size(B,2)
    Ubar = zeros(p,mx);
    @inbounds for i = n:-1:1
        tmpU = Ut[:,i];
        tmpW = Wt[:,i];
        X[i,:] = (B[i:i,:] - tmpW'*Ubar)/(tmpU'*tmpW);
        Ubar += tmpU .* X[i:i,:];
    end
end

#### Triangular product (Lx) ####
function ss_tri_mul!(Y::AbstractArray,K::SymEGRSSCholesky,X::AbstractArray)
     p, n = size(K.Ut);
     mx = size(X,2);
     Wbar = zeros(p,mx);
     for i = 1:n
         tmpW   = K.Wt[:,i]
         tmpU   = K.Ut[:,i]
         tmpX   = X[i:i,:]
         Wbar  += tmpW .* tmpX;
         Y[i,:] = Wbar'*tmpU;
     end
	 return Y
end
#### Adjoint triangular product (L'x) ####
function ssa_tri_mul!(Y::AbstractArray,K::SymEGRSSCholesky,X::AbstractArray)
     p, n = size(K.Ut);
	 mx = size(X,2);
     Ubar = zeros(p,mx);
     Ubar = Ubar + K.Ut*X;
     for i = 1:n
         tmpW = K.Wt[:,i]
         tmpU = K.Ut[:,i]
         tmpX = X[i:i,:]
         Y[i,:] = Ubar'*tmpW;
         Ubar  -= tmpU .* tmpX
     end
	 return Y
end

########################################################################
#### More linear Algebra routines 							   		####
########################################################################
function det(L::SymEGRSSCholesky)
    dd = one(eltype(L))
    @inbounds for i in 1:L.n
        dd *= dot(L.Ut[:,i],L.Wt[:,i])
    end
    return dd
end

function logdet(L::SymEGRSSCholesky)
    dd = zero(eltype(L))
    @inbounds for i in 1:L.n
        dd += log(dot(L.Ut[:,i],L.Wt[:,i]))
    end
    dd
end

#
mul!(y::AbstractVecOrMat, K::SymEGRSSCholesky, x::AbstractVecOrMat) =
	ss_tri_mul!(y,K,x)
mul!(y::AbstractVecOrMat, K::Adjoint{<:Any,<:SymEGRSSCholesky}, x::AbstractVecOrMat) =
	ssa_tri_mul!(y,K.parent,x)
function (\)(F::SymEGRSSCholesky, B::AbstractVecOrMat)
	X = similar(B)
	ss_forward!(X,F.Ut,F.Wt,B)
	return X
end
function (\)(F::Adjoint{<:Any,<:SymEGRSSCholesky}, B::AbstractVecOrMat)
	Y = similar(B)
	ssa_backward!(Y,F.parent.Ut,F.parent.Wt,B)
	return Y
end

ext: jl

max_stars_repo_name: mipals/SymEGRSSMatrices


text:
module Benchmarkplotting

using BenchmarkTools
using Gadfly
using DataFrames
using Printf
using StatsBase
using Statistics

export bcompare, report

get_names(::NamedTuple{Names, T}) where {Names, T} = Names

function bcompare(criterion :: Function,
                  data :: Vector{Pair{Symbol, T}},
                  implementations :: Vector{Pair{Symbol, Function}};
                  repeat :: Int = 2,
                  quiet :: Bool = false) where T
    rows = []
    for (impl_name, impl) in implementations
        for (case_name, case) in data
            # avoid the evaluation of case(if it's of AST types)
            @inline function action()
                impl(case)
            end
            if !quiet
                @info impl_name case_name action()
            end
            fst_test = criterion(@benchmark $action())
            fields = get_names(fst_test)
            repeats = [criterion(@benchmark $action) for _ in 1:repeat]
            res = Dict([field => fst_test[field] for field in fields])
            for field in fields
                res[field] = (res[field] + sum([it[field] for it in repeats])) ./ (1 + repeat)
            end
            row = (implementation = impl_name,
                   case = case_name,
                   res...)
            push!(rows, row)
        end
    end
    DataFrame(rows)
end

function report(field :: Symbol,
                result :: DataFrame,
                layouts...)
    result = copy(result)
    impls = []
    means = []

    case_names  = unique(result.case)
    case_names_with_unit = Dict()
    for case_name in case_names
        idx = result[:case] .== case_name
        tmp = result[idx, field]
        minval = minimum(tmp)
        tmp = tmp ./ minval
        result[idx, field] = tmp
        gmean = geomean(tmp)
        push!(means, gmean)
        minval_repr = @sprintf "%.3f" minval
        case_names_with_unit[case_name] = Symbol(case_name, "(min: $minval_repr)")
    end

    casemean = DataFrame(case=case_names, geomean=means)
    benchmarks = join(result, casemean, on=:case)

    benchmarks[:, :case] = map(x -> case_names_with_unit[x], benchmarks[:, :case])
    casemean[:, :case] = map(x -> case_names_with_unit[x], casemean[:, :case])

    plot(benchmarks,
        x = :case,
        y = field,
        color=:implementation,
        Guide.ylabel(nothing),
        Guide.xlabel(nothing),
        Theme(
            guide_title_position = :left,
            colorkey_swatch_shape = :circle,
            minor_label_font = "Consolas",
            major_label_font = "Consolas"
         ),
         layouts...,
    ), casemean[:, 1:2]

end

end # module


text:
using DynamicalSystemsBase
using DelayEmbeddings
using Test, DelimitedFiles

println("\nTesting garcia_almeida.jl...")
@testset "Garcia & Almeida method" begin

# Test Lorenz example
# For comparison reasons using Travis CI we carry out the integration on a UNIX
# OS and save the resulting time series
# See https://github.com/JuliaDynamics/JuliaDynamics for the storage of
# the time series used for testing
#
# u0 = [0, 10.0, 0.0]
# lo = Systems.lorenz(u0; σ=10, ρ=28, β=8/3)
# tr = trajectory(lo, 100; Δt = 0.01, Ttr = 100)
tr = readdlm(joinpath(tsfolder, "test_time_series_lorenz_standard_N_10000_multivariate.csv"))
tr = Dataset(tr)

x = tr[:, 1]
Y = Dataset(x)
Y = standardize(Y)

@testset "N-statistic" begin
    τs = 0:100
    N , NN_distances = n_statistic(Y, x, w=0, T=1, τs = τs)

    T = 17
    N2 , NN_distances2 = n_statistic(Y, x, w=0, T=T, τs = τs)

    # check whether the `d_E1`-statistic is the same
    @test NN_distances[1][1:50] == NN_distances2[1][1:50]
    @test NN_distances[5][1:100] == NN_distances2[5][1:100]

    (max_1_idx,_) = DelayEmbeddings.findlocalextrema(N)
    (max_2_idx,min_2_idx) = DelayEmbeddings.findlocalextrema(N2)

    @test 22 ≤ max_1_idx[5] ≤ 24
    @test 35 ≤ max_1_idx[6] ≤ 37
    @test 79 ≤ max_1_idx[12] ≤ 81

    @test 11 ≤ min_2_idx[1] ≤ 13
    @test 36 ≤ max_2_idx[3] ≤ 38
    @test 79 ≤ max_2_idx[7] ≤ 81

    # # plot N-Statistic for the Lorenz system as in Fig. 2(a) in [^Garcia2005a]
    # using Plots
    # plot(N,linewidth = 2, label = "T=1 (as in the Paper)", xaxis=:log)
    # plot!(N2,linewidth = 2, label = "T=$T", xaxis=:log)
    # plot!(N,seriestype = :scatter, label = "",xaxis=:log)
    # plot!(N2,seriestype = :scatter, label = "",xaxis=:log)
    # plot!(title = "Lorenz System as in Fig. 2(a) in the Paper")
    # xlabel!("τ")
    # ylabel!("N-Statistic")
end

@testset "garcia_embed univariate" begin
    Y_act, τ_vals, ts_vals, FNNs, NS = garcia_almeida_embedding(x; τs=0:100,  w = 17, T = 17)
    Y_act2, τ_vals2, ts_vals2, FNNs2, NS2 = garcia_almeida_embedding(x; τs=0:100,  w = 1, T = 1)
    Y_act3, τ_vals3, ts_vals3, FNNs3, NS3 = garcia_almeida_embedding(x; τs=0:100,  w = 17, T = 17, fnn_thres=0.1)

    @test size(Y_act,2) == 3
    @test size(Y_act2,2) == 3
    @test 15 ≤ τ_vals[2] ≤ 17
    @test 1 ≤ τ_vals2[2] ≤ 3
    @test FNNs[1] ≥ FNNs[2] ≤ FNNs[3]
    @test FNNs2[end]>FNNs2[end-1]
    @test size(Y_act3,2) == 2
    @test τ_vals3 == τ_vals[1:2]

    # using Plots
    # plot3d(Y_act[:,1],Y_act[:,2],Y_act[:,3], marker=2, camera = (6, 4))
    # plot!(title = "Reconstruction of Lorenz attractor using Garcia&Almeida univariate embedding")
    #
    # plot3d(Y_act2[:,1],Y_act2[:,2],Y_act2[:,3], marker=2, camera = (6, 4))
    # plot!(title = "Reconstruction of Lorenz attractor using Garcia&Almeida univariate embedding as in the paper")
end

@testset "garcia_embed multivariate" begin

    Y_act, τ_vals, ts_vals, FNNs, NS = garcia_almeida_embedding(tr; τs=0:100,  w = 17, T = 17)
    Y_act2, τ_vals2, ts_vals2, FNNs2, NS2 = garcia_almeida_embedding(tr; τs=0:100,  w = 1, T = 1)

    @test size(Y_act,2) == 5
    @test size(Y_act2,2) == 3
    @test ts_vals[1] == ts_vals[4] == 1
    @test ts_vals[2] == ts_vals[5] == 3
    @test ts_vals[3] == 2
    @test τ_vals == [0, 5, 5, 38, 55]

    @test ts_vals2 == [1, 1, 1]
    @test τ_vals2 == [0, 1, 2]

    # try to reproduce Fig.2a in [^Garcia2005b]
    tra = Dataset(hcat(tr[:,1], tr[:,3]))
    tra = standardize(tra)
    taus = 0:100

    N , _ = n_statistic(tra, tra[:,1], w=17, T=17, τs = taus)
    N2 , _ = n_statistic(tra, tra[:,2], w=17, T=17, τs = taus)

    _, minis2 = findmin(N2)

    @test 20 ≤ taus[minis2] ≤ 30

    # using Plots
    # NN = hcat(N,N2)
    # plot(NN, xaxis=:log)
    # plot(NN)

end

end


text:
function _precompile_()
    #    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(searchsortedfirst, (Array{Int64,1}, Base.StepRange{Int64,Int64}, Int64, Int64))
    precompile(searchsortedfirst, (Array{Int64,1}, Int64, Int64, Int64))
    precompile(searchsortedfirst, (Array{Int64,1}, Array{Int64,1}, Int64, Int64))
    precompile(length, (RLEVector{Int64,Int64},))
    precompile(numruns, (Array{Int64,1}, Array{Int64,1}))
    precompile(getindex, (RLEVector{Int64,Int64}, Base.UnitRange{Int64}))
    precompile(ind2run, (RLEVector{Int64,Int64}, Base.UnitRange{Int64}))
    precompile(numruns, (Array{Int64,1},))
    precompile(ree, (Array{Int64,1}, Array{Int64,1}, Int64))
    precompile(ind2runcontext, (RLEVector{Int64,Int64}, Int64))
    precompile(getindex, (RLEVector{Int64,Int64}, Base.StepRange{Int64,Int64}))
    precompile(inverse_ree, (Array{Int64,1}, Array{Int64,1}))
    precompile(ree, (Array{Int64,1},))
    precompile(ind2runcontext, (RLEVector{Int64,Int64}, Base.UnitRange{Int64}))
    precompile(setindex!, (RLEVector{Int64,Int64}, Int64, Int64))
    precompile(setindex!, (RLEVector{Int64,Int64}, Int64, Base.UnitRange{Int64}))
    precompile(Type, (Type{RLEVector}, Array{Int64,1}, Array{Int64,1}))
    precompile(setindex!, (RLEVector{Int64,Int64}, Array{Int64,1}, Array{Bool,1}))
    precompile(ind2run, (RLEVector{Int64,Int64}, Array{Int64,1}))
    precompile(ind2run, (RLEVector{Int64,Int64}, Int64))
    precompile(getindex, (RLEVector{Int64,Int64}, Array{Int64,1}))
    precompile(ind2run, (RLEVector{Int64,Int64}, Base.StepRange{Int64,Int64}))
    precompile(Type, (Type{RLEVector}, Array{Int64,1}))
end


text:
using CUDA

# GPU implementation
for S in (32,64)
    T = Symbol("Float$S")
    forw = Symbol("dropout_$S")
    back = Symbol("dropback_$S")
    @eval begin
        function dropout!(p::Number, x::CuArray{$T}, y::CuArray{$T})
            rand!(y)
            @knet8($forw,(Cint,$T,CuPtr{$T},CuPtr{$T}),length(y),$T(p),x,y)
            return y
        end
        function dropback!(p::Number, x::CuArray{$T}, y::CuArray{$T}, dy::CuArray{$T}, dx::CuArray{$T})
            @knet8($back,(Cint,$T,CuPtr{$T},CuPtr{$T},CuPtr{$T},CuPtr{$T}),length(dx),$T(p),x,y,dy,dx)
            return dx
        end
    end
end


text:
export LandDomainBC,
    LandComponentBC,
    Dirichlet,
    Neumann,
    NoBC,
    SurfaceDrivenWaterBoundaryConditions,
    SurfaceDrivenHeatBoundaryConditions

"""
   AbstractBoundaryConditions
"""
abstract type AbstractBoundaryConditions end

"""
    NoBC <: AbstractBoundaryConditions

This type is used for dispatch when no boundary condition needs
to be enforced - for example, if no prognostic variables are included
for a subcomponent, or for lateral faces when a 1D vertical setup
is used.
"""
struct NoBC <: AbstractBoundaryConditions end

"""
    Dirichlet{Fs} <: AbstractBoundaryConditions

A concrete type to hold the state variable
function, if Dirichlet boundary conditions are desired.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct Dirichlet{Fs} <: AbstractBoundaryConditions
    "state boundary condition"
    state_bc::Fs
end

"""
    Neumann{Ff} <: AbstractBoundaryConditions

A concrete type for specifying the magnitude of an inward diffusive flux,
normal to the domain boundary, from which a Neumann boundary condition
is determined and applied.

For the soil water model, the Darcy flux is defined as F⃗(Darcy) = -K∇h.
With the `Neumann` boundary
condition type, the user supplies a scalar `f` such that the applied flux
boundary conditon is F⃗(Darcy, boundary) = -f n̂,
where n̂ is the normal vector pointing out of the domain. For the soil heat model,
the same rule applies, except now we have a conductive heat flux
F⃗(Conductive, boundary) = -κ∇T = -f n̂.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct Neumann{Ff} <: AbstractBoundaryConditions
    "Scalar flux boundary condition"
    scalar_flux_bc::Ff
end

"""
    SurfaceDrivenWaterBoundaryConditions{FT, PD, RD} <: AbstractBoundaryConditions

Boundary condition type to be used when the user wishes to
apply physical fluxes of water at the top of the domain, normal to the surface.
The flux can be inwards or outwards, depending on the magnitude of infiltration,
evaporation, and precipitation.

Precipitation is assumed to be
in the ẑ direction, P⃗ = Pẑ, while evaporation is along the normal
direction, E⃗ = En̂. The applied boundary flux is F⃗ = f n̂, and the magnitude f is
determined internally from P⃗, E⃗, and soil conditions.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct SurfaceDrivenWaterBoundaryConditions{FT, PD, RD} <:
       AbstractBoundaryConditions where {FT, PD, RD}
    "Precipitation model"
    precip_model::PD
    "Runoff model"
    runoff_model::RD
end

"""
    SurfaceDrivenWaterBoundaryConditions(
        ::Type{FT};
        precip_model::AbstractPrecipModel{FT} = DrivenConstantPrecip{FT}(),
        runoff_model::AbstractSurfaceRunoffModel{FT} = NoRunoff(),
    ) where {FT}

Constructor for the SurfaceDrivenWaterBoundaryConditions object. The default
is a constant precipitation rate on the subgrid scale, and no runoff.
"""
function SurfaceDrivenWaterBoundaryConditions(
    ::Type{FT};
    precip_model::AbstractPrecipModel{FT} = DrivenConstantPrecip{FT}(
        (t) -> (0.0),
    ),
    runoff_model::AbstractSurfaceRunoffModel = NoRunoff(),
) where {FT}
    args = (precip_model, runoff_model)
    return SurfaceDrivenWaterBoundaryConditions{FT, typeof.(args)...}(args...)
end

"""
    SurfaceDrivenHeatBoundaryConditions{FT, SWD} <: AbstractBoundaryConditions

Boundary condition type to be used when the user wishes to
apply normal fluxes of heat at the top of the domain
(according to radiative energy fluxes).

# Fields
$(DocStringExtensions.FIELDS)
"""
struct SurfaceDrivenHeatBoundaryConditions{FT, SWD} <:
       AbstractBoundaryConditions where {FT, SWD}
    "Net short wave flux model"
    nswf_model::SWD
end

"""
    SurfaceDrivenHeatBoundaryConditions(
        ::Type{FT};
        nswf_model::AbstractNetSwFluxModel{FT} = PrescribedNetSwFlux{FT}()
    ) where {FT}

Constructor for the SurfaceDrivenHeatBoundaryConditions object. The default
is no shortwave flux.
"""
function SurfaceDrivenHeatBoundaryConditions(
    ::Type{FT};
    nswf_model::AbstractNetSwFluxModel{FT} = PrescribedNetSwFlux(
        FT;
        nswf = t -> eltype(t)(0),
    ),
) where {FT}
    return SurfaceDrivenHeatBoundaryConditions{FT, typeof(nswf_model)}(
        nswf_model,
    )
end

"""
    LandDomainBC{TBC, BBC, LBC}

A container for the land boundary conditions, with options for surface
boundary conditions, bottom boundary conditions, or lateral face
boundary conditions.

The user should supply an instance of `LandComponentBC` for each piece, as needed.
If none is supplied, the default for is `NoBC` for each subcomponent. At a minimum, both top and bottom boundary conditions should be supplied. Whether or not to include the lateral faces depends on the configuration of the domain.
"""
Base.@kwdef struct LandDomainBC{TBC, BBC, MinXBC, MaxXBC, MinYBC, MaxYBC}
    "surface boundary conditions"
    surface_bc::TBC = LandComponentBC()
    "bottom boundary conditions"
    bottom_bc::BBC = LandComponentBC()
    "lateral boundary conditions"
    minx_bc::MinXBC = LandComponentBC()
    maxx_bc::MaxXBC = LandComponentBC()
    miny_bc::MinYBC = LandComponentBC()
    maxy_bc::MaxYBC = LandComponentBC()
end

"""
    LandComponentBC{SW, SH, SF}

An object that holds the boundary conditions for each of the subcomponents
of the land model.

The boundary conditions supplied should be of type `AbstractBoundaryConditions`.
The default is `NoBC` for each component, so that the user only
needs to define the BC for the components they wish to model.
"""
Base.@kwdef struct LandComponentBC{
    SW <: AbstractBoundaryConditions,
    SH <: AbstractBoundaryConditions,
    SF <: AbstractBoundaryConditions,
}
    soil_water::SW = NoBC()
    soil_heat::SH = NoBC()
    surface::SF = NoBC()
end

"""
    function boundary_conditions(land::LandModel)

Unpacks the `boundary_conditions` field of the land model, and
puts into the correct order based on the integers used to identify
faces, as defined in the Driver configuration.
"""
function boundary_conditions(land::LandModel)
    bc = land.boundary_conditions
    mytuple = (
        bc.bottom_bc,
        bc.surface_bc,
        bc.minx_bc,
        bc.maxx_bc,
        bc.miny_bc,
        bc.maxy_bc,
    )
    # faces labeled integer 1,2 are bottom, top, lateral sides are 3, 4, 5, 6
    return mytuple
end

function boundary_state!(
    nf,
    bc::LandComponentBC,
    land::LandModel,
    state⁺::Vars,
    aux⁺::Vars,
    n,
    state⁻,
    aux⁻,
    t,
    args...,
)
    land_boundary_state!(
        nf,
        bc,
        land,
        state⁺,
        aux⁺,
        n,
        state⁻,
        aux⁻,
        t,
        args...,
    )
end

function land_boundary_state!(nf, bc::LandComponentBC, land, args...)
    soil_boundary_state!(nf, bc.soil_water, land.soil.water, land, args...)
    soil_boundary_state!(nf, bc.soil_heat, land.soil.heat, land, args...)
    surface_boundary_state!(nf, bc.surface, land.surface, land, args...)
end

function boundary_state!(
    nf,
    bc::LandComponentBC,
    land::LandModel,
    state⁺::Vars,
    diff⁺::Vars,
    hyperdiff⁺::Vars,
    aux⁺::Vars,
    n,
    state⁻,
    diff⁻,
    hyperdiff⁻,
    aux⁻,
    t,
    args...,
)
    land_boundary_flux!(
        nf,
        bc,
        land,
        state⁺,
        diff⁺,
        aux⁺,
        n,
        state⁻,
        diff⁻,
        aux⁻,
        t,
        args...,
    )
end

function land_boundary_flux!(nf, bc::LandComponentBC, land, args...)
    soil_boundary_flux!(nf, bc.soil_water, land.soil.water, land, args...)
    soil_boundary_flux!(nf, bc.soil_heat, land.soil.heat, land, args...)
    surface_boundary_flux!(nf, bc.surface, land.surface, land, args...)
end


text:
const _GW_TASK_README_FILE_NAME = "README.md"
_readme_file(taskdir::String) = joinpath(taskdir, _GW_TASK_README_FILE_NAME)
_readme_file(gwt::GWTask) = _readme_file(task_dir(gwt))

function _readme_txt(;
        tid, src, desc, lang,
    )
    return string(
        "# ", tid, "\n\n",
        desc, "\n\n",
        "## Usage", "\n\n",
        "- Run `runme.jl` for a full reproduction", "\n",
        "- Check `out.log` for last run output", "\n",
        "- Check `run.log` for run history", "\n",
        "\n",
        "## Source code", "\n\n",
        "```", lang, "\n",
        src, "\n",
        "```",
        "\n"
    )
end

# gen a readme from task data
function _do_readme!(gwt::GWTask)
    readme = _readme_txt(;
        tid  = task_id(gwt),
        src = _task_src(gwt),
        desc = _task_desc(gwt),
        lang = _task_lang(gwt)
    )
    mdfile = _readme_file(gwt)
    write(mdfile, readme)
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
using YAML
using LightGraphs, MetaGraphs
using DataFrames, CSV

import ModelParameters: update, update!, Model
using Waterflow

using Infiltrator

network = YAML.load_file("../tests/data/AWRA_R_Network/campaspe_network.yml")
g, mg = create_network("Example Network", network)

inlets, outlets = find_inlets_and_outlets(g)

@info "Network has the following inlets and outlets:" inlets outlets

# climate_data = DataFrame!(CSV.File("../tests/data/climate/climate_historic.csv",
#                           comment="#",
#                           dateformat="YYYY-mm-dd"))
# climate = Climate(climate_data, "_rain", "_evap")

# @info "Running example stream..."
# timesteps = sim_length(climate)
# for ts in (1:timesteps)
#     for outlet in outlets
#         run_node!(mg, g, outlet, climate, ts)
#     end
# end

max_stars_repo_name: ConnectedSystems/camgraph

max_issues_repo_name: ConnectedSystems/camgraph

max_forks_repo_name: ConnectedSystems/camgraph


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

"""
    EigenAnalysis(proj; maxdim=nothing, pratio=1.0)

The eigenanalysis of the covariance with a given projection `proj`.
Optionally specify the maximum number of dimensions in the output `maxdim`
and the percentage of variance to retain `pratio`. Default to all dimensions of
the input.

## Projections

* `:V` - Uncorrelated variables (PCA transform)
* `:VD` - Uncorrelated variables and variance one (DRS transform)
* `:VDV` - Uncorrelated variables and variance one (SDS transformation)

The `:V` projection used in the PCA transform projects the data on the eigenvectors
V of the covariance matrix.

The `:VD` projection used in the DRS transform. Similar to the `:V` projection,
but the eigenvectors are multiplied by the squared inverse of the eigenvalues D.

The `:VDV` projection used in the SDS transform. Similar to the `:VD` transform,
but the data is projected back to the basis of the original variables using the Vᵀ matrix.

See [https://geostatisticslessons.com/lessons/sphereingmaf](https://geostatisticslessons.com/lessons/sphereingmaf)
for more details about these three variants of eigenanalysis.

# Examples

```julia
EigenAnalysis(:V)
EigenAnalysis(:VD)
EigenAnalysis(:VDV)
EigenAnalysis(:V, maxdim=3)
EigenAnalysis(:VD, pratio=0.99)
EigenAnalysis(:VDV, maxdim=3, pratio=0.99)
```
"""
struct EigenAnalysis <: Transform
  proj::Symbol
  maxdim::Union{Int,Nothing}
  pratio::Float64

  function EigenAnalysis(proj, maxdim, pratio)
    @assert proj ∈ (:V, :VD, :VDV) "Invalid projection."
    @assert 0 ≤ pratio ≤ 1 "Invalid pratio."
    new(proj, maxdim, pratio)
  end
end

EigenAnalysis(proj; maxdim=nothing, pratio=1.0) =
  EigenAnalysis(proj, maxdim, pratio)

assertions(::Type{EigenAnalysis}) = [assert_continuous]

isrevertible(::Type{EigenAnalysis}) = true

function apply(transform::EigenAnalysis, table)
  # basic checks
  for assertion in assertions(transform)
    assertion(table)
  end

  # original columns names
  cols = Tables.columns(table)
  onames = Tables.columnnames(cols)

  # table as matrix
  X = Tables.matrix(table)

  # center the data
  μ = mean(X, dims=1)
  Y = X .- μ

  # eigenanalysis of covariance
  S, S⁻¹ = eigenmatrices(transform, Y)

  # project the data
  Z = Y * S

  # column names
  names = Symbol.(:PC, 1:size(Z, 2))

  # table with transformed columns
  𝒯 = (; zip(names, eachcol(Z))...)
  newtable = 𝒯 |> Tables.materializer(table)

  newtable, (μ, S, S⁻¹, onames)
end

function revert(::EigenAnalysis, newtable, cache)
  # table as matrix
  Z = Tables.matrix(newtable)

  # retrieve cache
  μ, S, S⁻¹, onames = cache

  # undo projection
  Y = Z * S⁻¹

  # undo centering
  X = Y .+ μ

  # table with original columns
  𝒯 = (; zip(onames, eachcol(X))...)
  𝒯 |> Tables.materializer(newtable)
end

function reapply(transform::EigenAnalysis, table, cache)
  # basic checks
  for assertion in assertions(transform)
    assertion(table)
  end

  # table as matrix
  X = Tables.matrix(table)

  # retrieve cache
  μ, S, S⁻¹, onames = cache

  # center the data
  Y = X .- μ

  # project the data
  Z = Y * S

  # column names
  names = Symbol.(:PC, 1:size(Z, 2))

  # table with transformed columns
  𝒯 = (; zip(names, eachcol(Z))...)
  𝒯 |> Tables.materializer(table)
end

_maxdim(maxdim::Int, Y) = maxdim
_maxdim(::Nothing, Y) = size(Y, 2)

function outdim(transform, Y, λ)
  pratio = transform.pratio
  csums  = cumsum(λ)
  ratios = csums ./ last(csums)
  mdim   = _maxdim(transform.maxdim, Y)
  pdim   = findfirst(≥(pratio), ratios)
  min(mdim, pdim)
end

function eigenmatrices(transform, Y)
  proj = transform.proj

  Σ = cov(Y)
  λ, V = eigen(Σ, sortby=λ -> -real(λ))

  if proj == :V
    S   = V
    S⁻¹ = transpose(V)
  elseif proj == :VD
    Λ   = Diagonal(sqrt.(λ))
    S   = V * inv(Λ)
    S⁻¹ = Λ * transpose(V)
  elseif proj == :VDV
    Λ   = Diagonal(sqrt.(λ))
    S   = V * inv(Λ) * transpose(V)
    S⁻¹ = V * Λ * transpose(V)
  end

  d = outdim(transform, Y, λ)

  S[:, 1:d], S⁻¹[1:d, :]
end

"""
    PCA(; maxdim=nothing, pratio=1.0)

The PCA transform is a shortcut for
`ZScore() → EigenAnalysis(:V; maxdim, pratio)`.

See also: [`ZScore`](@ref), [`EigenAnalysis`](@ref).

# Examples

```julia
PCA(maxdim=2)
PCA(pratio=0.86)
PCA(maxdim=2, pratio=0.86)
```
"""
PCA(; maxdim=nothing, pratio=1.0) =
  ZScore() → EigenAnalysis(:V, maxdim, pratio)

"""
    DRS(; maxdim=nothing, pratio=1.0)

The DRS transform is a shortcut for
`ZScore() → EigenAnalysis(:VD; maxdim, pratio)`.

See also: [`ZScore`](@ref), [`EigenAnalysis`](@ref).

# Examples

```julia
DRS(maxdim=3)
DRS(pratio=0.87)
DRS(maxdim=3, pratio=0.87)
```
"""
DRS(; maxdim=nothing, pratio=1.0) =
  ZScore() → EigenAnalysis(:VD, maxdim, pratio)

"""
    SDS(; maxdim=nothing, pratio=1.0)

The SDS transform is a shortcut for
`ZScore() → EigenAnalysis(:VDV; maxdim, pratio)`.

See also: [`ZScore`](@ref), [`EigenAnalysis`](@ref).

# Examples

```julia
SDS()
SDS(maxdim=4)
SDS(pratio=0.88)
SDS(maxdim=4, pratio=0.88)
```
"""
SDS(; maxdim=nothing, pratio=1.0) =
  ZScore() → EigenAnalysis(:VDV, maxdim, pratio)


text:
module EasyEvents

include("AbstractEvent.jl")
export AbstractEvent
export old_state, new_state, update!, reset!, trigger, states, istraking
export has_event, has_event!, on_event, on_event!

include("CustomEvent.jl")
export CustomEvent
export extras

include("FileMTimeEvent.jl")
export FileMTimeEvent

include("FileSizeEvent.jl")
export FileSizeEvent

include("FileContentEvent.jl")
export FileContentEvent

include("PropertyChangedEvent.jl")
export PropertyChangedEvent

end


text:
module Ops20_gpu

include("cudnn_retry.jl")
include("activation.jl")
include("batchnorm.jl")
include("bmm.jl")
include("conv.jl")
include("dropout.jl")
include("softmax.jl")
include("loss.jl")
include("rnn.jl")

end


text:
using GenericLinearAlgebra, DelimitedFiles,Random
# State of the
mutable struct State
    f   # function
    discr
    eltype
    #
    graph
    cref
    #
    params # Dict with keywords
end
function State(f,eltype)
    return State(f,[],eltype,Nothing,Nothing,Dict{Symbol,Any}(:opt_kwargs => Dict()));
end

function scale_and_square!(state)
    degopt=Degopt(state.graph)
    scale!(degopt,1/2)
    square!(degopt)
    (g,c)=graph_degopt(degopt)
    state.graph=g
    state.cref=c
end

function get_disc_discr(state;n=state.params[:n])
    rho=state.params[:rho]
    nn=round(Int,n);
    discr=rho*exp.(1im*range(0,2*pi,length=nn)[1:end-1]);
    discr=convert.(state.eltype,discr);
    return discr;
end
function get_taylor(eltype,n)
    c=1 ./factorial.(convert.(big(real(eltype)),0:(n-1)));
    c=convert.(eltype,c);
end

# Corresponds to a call to opt_gaussnewton!
struct OptSimulation
    kwargs
    cref
    special_case
end

function OptSimulation(s::State;special_case=0)
    return OptSimulation(s.params[:opt_kwargs],s.cref,special_case);
end

# Levenberg-Marquard simulations using LsqFit.curve_fit
struct LMSimulation
    kwargs
end

struct NoOp
end

function LMSimulation(s::State)
    kwargs=Dict();
    if (haskey(s.params,:curve_fit_kwargs))
        kwargs=s.params[:curve_fit_kwargs]
    end
    return LMSimulation(kwargs)
end

# Corresponds to changing a parameter
struct ModifyParam
    param
    factor
end

# Corresponds to a change of all variables
struct KickIt
    factor::Number
    mode
end
function KickIt(state::State,factor)

    if haskey(state.params,:kickit_mode)
        mode=state.params[:kickit_mode]
    else
        mode=0;
    end
    @show mode
    return KickIt(factor,mode);

end
struct ZeroNormalize
end

struct ShowInfo
    theta_guess;
end

struct LineSearch
end

struct SetSmallToZero
    k
end
struct PurgeZeroOpt
    k
end

struct DegoptNormalize
end

# Initialize State with a graph based on number of multiplications.
function init_state_mult(f,rho,n,mult; eltype=Complex{BigFloat})
    init_state_mult(nothing,f,rho,n,mult; eltype=eltype);
end

function init_state_mult(graphname,f,rho,n,mult; eltype=Complex{BigFloat})
    if (f != exp)
        error("Only exp implemented");
    end
    state=State(f,eltype);
    state.params[:n]=n;
    state.params[:target_n]=n;
    state.params[:rho]=rho;
    state.params[:target_n]=n;

   state.discr=get_disc_discr(state);

    init_state_mult!(state,graphname,mult);
    return state;
end

function init_state_mult!(state,graphname,mult;showmeta=false)

    if (!isnothing(graphname))
        print("Generating $graphname ");
        if (graphname == :sid)
            (graph,cref)=graph_sid_exp(mult)
        elseif (graphname == :bbc)
            (graph,cref)=graph_bbc_exp(m)
        elseif (graphname == :sastre)
            method=:y1s;
            if (m==6)
                method=:h2m
            elseif (m==8)
                method=:z1ps
            end
            (graph,cref)=graph_sastre_exp(m,method)
        end

        if (graphname in [:ps,:mono])

            taylorcoeffs=get_taylor(state.eltype,200);

            # Automatically detect degree
            graph_mult = 0;
            deg0 = m-1;
            graph=nothing; oldgraph=nothing;
            cref=nothing; oldcref=nothing;
            while graph_mult<mult+1
                deg0 += 1;
                c=taylorcoeffs[1:(deg0+1)];
                oldgraph=graph;
                oldcref=cref;
                if (graphname == :ps)
                    (graph,cref)=graph_ps_degopt(c)
                elseif (graphname == :mono)
                    (graph,cref)=graph_monomial_degopt(c)
                end
                graph_mult=count(values(graph.operations) .==:mult);
            end
            # We go one step too far in the automatic determination so
            # use the old graph
            graph=oldgraph;
            cref=oldcref;

            println(" degree $(deg0-1)");
        end

        state.params[:graphname]="$(graphname)_m$(mult)"
        graph = Compgraph(state.eltype,graph);
        state.graph=graph;
        state.cref=cref;
        if (showmeta)
            showerr(state,showmeta=true)
        end

    end
    return state;
end
function init_state_file!(state,graphname,filename;showmeta=true,scale_and_square=false);

    if (filename isa String)
        state.graph=import_compgraph(filename);
    else
        state.graph=filename;
    end

    mult=count(values(state.graph.operations) .==:mult);
    if (scale_and_square)
        mult += 1;
    end
    state.params[:graphname]="$(graphname)_m$(mult)";
    if (scale_and_square)
        scale_and_square!(state)
    else
        (_,cref)=graph_degopt(Degopt(state.graph))
        state.cref=cref;
    end

    if (showmeta)
        println("Generating $graphname based on $filename");
        showerr(state,showmeta=true)
    end

    return state;
end

# initcommand(control::String,state)  --> command object
function initcommand(control,state)
    if (control == "s")
        return OptSimulation(state);
    elseif (control == "S")
        return OptSimulation(state;special_case=1);
    elseif (control == "2")
        return OptSimulation(state;special_case=2);
    elseif (control == "3")
        return OptSimulation(state;special_case=3);
    elseif (control == "4")
        return OptSimulation(state;special_case=4);
    elseif (control == "5")
        return OptSimulation(state;special_case=5);
    elseif (control == "l")
        return LMSimulation(state);
    elseif (control == "L")
        return LineSearch();
    elseif (control == "g")
        return ModifyParam(:γ0,0.8);
    elseif (control == "G")
        return ModifyParam(:γ0,1/0.8);
    elseif (control=="n")
        return ModifyParam(:n,1/2);
    elseif (control=="Z")
        return SetSmallToZero(1);
    elseif (control=="N")
        return ModifyParam(:n,2);
    elseif (control=="D")
        return ModifyParam(:droptol,2);
    elseif (control=="d")
        return ModifyParam(:droptol,1/2);
    elseif (control=="k")
        return KickIt(state,1e-9);
    elseif (control=="K")
        return KickIt(state,1e-6);
    elseif (control=="z")
        return ZeroNormalize();
    elseif (control=="p") # Purge crefs
        return PurgeZeroOpt(1);
    elseif (control=="0") #
        return DegoptNormalize();
    elseif (control=="t")
        return ShowInfo(Nothing);
    elseif (control=="T")
        return ShowInfo(0.3);
    else
        return NoOp();
    end
end

# runcommand(commandobj, state) -> State
function runcommand(s::OptSimulation,state)
    # Run Gauss-Newton with params in s.params[:opt_kwargs]
    f=state.f;
    graph=deepcopy(state.graph);
    discr=state.discr;
    cref=deepcopy(s.cref);
    if (s.special_case == 1)
        discr=deepcopy(discr);
        push!(discr,zero(eltype(discr)));
    end
    if (s.special_case == 2)
        discr0=discr;
        discr=deepcopy(discr);
        append!(discr,discr0/2);
        append!(discr,0.1*discr0/maximum(abs.(discr0)));
    end
    if (s.special_case == 3)
        discr=discr/2;
    end
    if (s.special_case == 4)
        discr=discr/10;
    end
    if (s.special_case == 5)
        c=state.params[:cref_count]
        cref=[cref[mod(c,size(cref,1))+1]]
        state.params[:cref_count] += 1;
        println("count c:$c $(cref[1])");
    end

    opt_gauss_newton!(graph,f,discr;logger=1,
                      stoptol=1e-40,cref=cref,
                      s.kwargs...);

    ff=f.(discr);
    ftilde=eval_graph(graph,discr);
    errv=abs.((ff-ftilde)./ftilde);

    state=deepcopy(state);

    state.graph=graph # Set the output
    currentname=state.params[:graphname]
    if (!contains(currentname,"_opt"))
        state.params[:graphname]=currentname*"_opt";
    end

    return state;
end

function runcommand(s::ModifyParam,state)
    # Update a parameter some opt_kwargs treated separately
    state=deepcopy(state);
    param = s.param;
    factor = s.factor;
    if (param in [:droptol,:γ0])
        o=state.params[:opt_kwargs]
    else
        o=state.params;
    end
    T=typeof(o[param])
    o[param]=o[param]*factor;
    println("$param=$(o[param])");
    if (param == :n)
        state.discr=get_disc_discr(state);
    end

    return state;
end

function runcommand(s::LineSearch,state)

    graph0=deepcopy(state.graph);
    smallest_err=Inf;
    smallest_k=Inf;
    smallest_graph=deepcopy(state.graph);
    for (k,cref) in enumerate(state.cref)
        cref=state.cref[1];
        err0=showerr(state,output=false);

        v0=get_coeffs(state.graph,[cref])[1];

        vv=[v0];
        append!(vv,v0*(1 .+ 10 .^(-20.0:2:-1)))
        append!(vv,v0*(1 .- 10 .^(-20.0:2:-1)))
        #vv=v0*(1 .+ [0;])
        errv=zeros(typeof(err0),size(vv));
        for (i,v) in enumerate(vv)
            state.graph=deepcopy(graph0);
            set_coeffs!(state.graph,v,[cref]);
            errv[i]=showerr(state,output=false);
        end
        j=argmin(errv);
        state.graph=deepcopy(graph0);
        set_coeffs!(state.graph,vv[j],[cref]);
        if (errv[j] <= smallest_err)
            #@show errv
            smallest_err=errv[j];
            smallest_k=k;
            smallest_graph=deepcopy(state.graph);
        end
    end
    println("Modifying $(state.cref[smallest_k])");
    state.graph=smallest_graph;
    return state;

end
function runcommand(s::SetSmallToZero,state)
    num_elements=s.k
    state=deepcopy(state);
    graph=deepcopy(state.graph);
    #crefs=get_all_cref(graph);
    crefs=deepcopy(state.cref);
    vals=abs.(get_coeffs(graph,crefs));
    J=sortperm(vals)
    removed=Float64(norm(vals[J[1:num_elements]]));
    println("zeroing: $(removed)");
    set_coeffs!(graph,zero(vals[J[1:num_elements]]),
                crefs[J[1:num_elements]]);
    state.graph=graph;

    vals=abs.(get_coeffs(graph,crefs));
    J=findall(vals .== 0);

    return state;
end

function runcommand(s::PurgeZeroOpt,state)

    num_elements=s.k
    state=deepcopy(state);
    graph=deepcopy(state.graph);
    #crefs=get_all_cref(graph);
    crefs=deepcopy(state.cref);
    vals=abs.(get_coeffs(graph,crefs));
    J=findall(vals .== 0);
    if (size(J,1)>0)
        println("Purging cref $(crefs[J[1:s.k]])");
        setdiff!(state.cref,crefs[J[1:s.k]]);
    else
        println("Found nothing to purge");
    end
    return state;
end
function runcommand(s::DegoptNormalize,state)
    state=deepcopy(state);
    degopt=Degopt(state.graph);
    normalize!(degopt);
    (g,_)=graph_degopt(degopt);
    state.graph=g;
    return state;
end

function curve_fit_model(x,coeffs,graph,cref)
    p=size(x,1)/2;
    org_coeffs=get_coeffs(graph,cref);
    set_coeffs!(graph,coeffs,cref)
    z=eval_graph(graph,x);
    #set_coeffs!(graph,org_coeffs,cref);
    return [real.(z);imag.(z)];
end

function curve_fit_jacobian_model(x,coeffs,graph,cref)
    org_coeffs=get_coeffs(graph,cref);
    set_coeffs!(graph,coeffs,cref)
    J=eval_jac(graph,x,cref);
    #set_coeffs!(graph,org_coeffs,cref);
    return [real.(J);imag.(J)]
end

function runcommand(s::LMSimulation,state)

    discr=state.discr;
    T=real(state.eltype);
    xdata=complex(T).(discr);
    ydata_c=(state.f.(discr));
    ydata=T.([real.(ydata_c);imag(ydata_c)])

    graph=Compgraph(T,state.graph);
    cref=state.cref;
    p0=get_coeffs(graph,cref);

    model = (x,coeffs) -> curve_fit_model(x,coeffs,graph,cref)
    jacobian_model =
        (x,coeffs) -> curve_fit_jacobian_model(x,coeffs,graph,cref)

    fit = curve_fit(model, jacobian_model, xdata, ydata, p0,#show_trace=true,
                    x_tol=big(1e-20), g_tol=1e-20, maxIter=10,
                    s.kwargs...)

    state.graph=graph # Set the output
    @show typeof(graph)
    currentname=state.params[:graphname]
    if (!contains(currentname,"_opt"))
        state.params[:graphname]=currentname*"_opt";
    end
    return state;

end
function runcommand(s::NoOp,state)
    println("No op");
    return state;
end

function runcommand(s::KickIt,state)
    factor=s.factor;
    (x,y)=get_degopt_crefs(state.graph);
    random_state=FakeRandomState(0);
    #Random.seed!(0);
    if (s.mode == 0)
        cref=[x[end][1][1]];
        vals=get_coeffs(state.graph,cref);
        vals=vals.*(1 .+ fake_randn(random_state,size(vals,1)));
        set_coeffs!(state.graph,vals,cref);

    elseif (s.mode == 1)
        cref=[x[end][1][1]];
        vals=get_coeffs(state.graph,cref);
        vals=vals .+ fake_randn(random_state,size(vals,1));
        set_coeffs!(state.graph,vals,cref);
    elseif (s.mode == 2)

        cref=[x[1][2][1]];
        vals=get_coeffs(state.graph,cref);
        vals .+= 0.0001
        set_coeffs!(state.graph,vals,cref);
    elseif (s.mode == 3)
        cref=[x[end][2][1]];
        vals=get_coeffs(state.graph,cref);
        vals=vals .+ 0.0001*fake_randn(random_state,size(vals,1));
        set_coeffs!(state.graph,vals,cref);
    elseif s.mode ==4
        y_cref=get_degopt_crefs(state.graph)[2];
        # Kick the y-degopt vector
        Ha=[ -0.000599134   1.00256       0.0          0.0           0.0         0.0         0.0
 -0.000118639  -0.000778662   0.999709     0.0           0.0         0.0         0.0
 -1.09012f-5   -0.000246103  -0.00126751   0.999822      0.0         0.0         0.0
 -1.08452f-6   -3.92744f-5   -0.00017262  -0.00189144    0.999981    0.0         0.0
  5.84573f-7   -4.80433f-6   -2.1741f-5   -0.000177471  -0.00239572  0.999999    0.0
  1.25109f-7   -1.05495f-7   -2.48853f-6  -4.14857f-5    4.32646f-5  0.00116684  1.0
             ];
        Hb=[-0.000599134  1.00256   0.0        0.0         0.0          0.0         0.0
 -0.000838943  1.0011    0.0166752  0.0         0.0          0.0         0.0
 -0.0019134    1.00034   0.0377923  0.00100261  0.0          0.0         0.0
 -0.00242811   0.999187  0.0604852  0.00308486  6.53079f-5   0.0         0.0
  0.00116408   0.997916  0.0755767  0.0051393   0.000163854  6.79643f-8  0.0
 -8.09922f-8   0.997063  0.0954523  0.00867858  0.000495314  1.32541f-5  3.94245f-8]
        crand=[0.999999803981589
               1.0006179487544273
               0.4979281278966493
               0.16592558923333828
               0.03871151836039848
               0.0061833894649526416
               0.0006066317390421412
               3.0487323320460782e-5]
        (gtmp,cref)=graph_degopt(Degopt(Ha,Hb,crand));
        cvals=get_coeffs(gtmp,cref);
        set_coeffs!(state.graph,cvals,cref);
        #p=min(size(y_cref,1),size(crand,1));
        #set_coeffs!(state.graph,crand[1:p],y_cref);

    elseif (s.mode == 5)
        crefs=state.cref;
        vals=get_coeffs(state.graph,crefs);
        vals=vals .+ 0.0001*fake_randn(random_state,size(vals,1));
        set_coeffs!(state.graph,vals,crefs);
    else

        println("Unknown Kickit Mode");
    end

    return state;
end

function runcommand(s::ZeroNormalize,state)
    graph=state.graph;
    print("Diff before:"*string(eval_graph(graph,0)-state.f(0)));
    degopt=Degopt(graph);
    x=degopt.x;
    y=degopt.y;
    y .*= state.f(0)/eval_graph(graph,0);
    (graph,_)=graph_degopt(x,y);
    set_coeffs!(state.graph,get_coeffs(graph));
    err_after=eval_graph(graph,0)-state.f(0);
    println(" after $err_after");
    return state;
end

function runcommand(s::ShowInfo,state)

    graph=state.graph;
    try

        theta_guess=state.params[:rho];
        if (s.theta_guess != Nothing)
            theta_guess = s.theta_guess
        end

        (ff,theta)=compute_bwd_theta_exponential(graph,coefftype=BigFloat,
                                                tolerance=eps()/2,theta_init=big(theta_guess));

        @show theta
        @show ff(0)
        @show ff(1e-10)
        @show ff(theta_guess)
    catch e
        println("Unable to compute theta: $(s.theta_guess)");
    end

    return state;

end

function run_sequence(state,predefsims="");
    max_j=1000;
    statelist=Vector{Any}(undef,max_j);
    commandlist=Vector{String}(undef,max_j);
    statelist[1]=state;
    pre_commandlist=split(predefsims,"");

    err=showerr(state,output=false);
    println("Initialerr $(state.params[:graphname]):$err");
    j=1;
    command="";
    while (command != "q")
        state=statelist[j] # Base it on the old state
        print("Commandlist: "*String(join(commandlist[1:(j-1)]))* " ");

        x="";
        if (j <= size(pre_commandlist,1)) # It's stored in
            x=pre_commandlist[j];
            pre_commandlist[j] = "";
        end
        if x == ""
            x=String(read_one_key());
        end

        commandlist[j]=x;
        println("$x");
        if (x=="b") # Back
            j=j-2;
        elseif (x == "q")
            break
        else
            command = initcommand(x,state)
            newstate = runcommand(command,state)
            statelist[j+1] = newstate;
        end

        showerr(statelist[j+1]);

        j += 1;

    end

    state=statelist[j];
    return (state,statelist[1:j],join(commandlist[1:j]))

end

function showmetainfo(state)
    println("rho: $(state.params[:graphname]): ρ=$(state.params[:rho]) ");
end

function showerr(state;output=true,showmeta=false,n=state.params[:target_n])
    f=state.f;
    discr=get_disc_discr(state,n=n);

    graph=state.graph;

    err=norm((f.(discr)-eval_graph(graph,discr))./f.(discr),Inf)
    imagnorm=Float64(norm(imag.(get_coeffs(graph))));
    if (output)
        if (showmeta)
            showmetainfo(state)
        end

        if (imagnorm>0)
            imagstr="norm(imag(coeffs))= $imagnorm"
        else
            imagstr="";
        end
        println("Target error: $(Float64(err)) $imagstr");
    end

    return err
end

function read_one_key(; io = stdin)
    setraw!(raw) = ccall(:jl_tty_set_mode, Int32, (Ptr{Cvoid},Int32), io.handle, raw)
    setraw!(true)
    s=read(io, 1)
    setraw!(false)
    return s
end

## Some generation requires random "kick-starting". The following
## is a tabulation of the random numbers generated by
## Julia Version 1.7.0-DEV.1241 such that the generation
## can work also if the random number generators change in later versions.
mutable struct FakeRandomState
    k
end
function fake_randn(f::FakeRandomState,s)

    #Random.seed!(f.k);
    data=Dict(1=> [0.20154471885274503],
              2=> [0.20154471885274503;1.0340551364982935],
              3=> [0.20154471885274503;1.0340551364982935;-0.5035805877040641],
              28=> [0.359629402423705
 -0.06242791627632372
 -1.399567921508011
  0.4475847763985716
 -0.8809915030490731
 -0.2017627452986302
 -0.5526516215876736
 -0.35572577945591155
 -0.1607797466665553
  1.1475981747910613
  1.5035230388445053
  0.2812185565541045
 -1.6580378604283104
 -0.5724292571431209
  0.11417065387535505
 -0.1722032235811361
  0.10954150466894817
 -0.20957790746640031
  0.9978114970316212
 -0.33835189120333353
  0.32711196563255096
  1.3252610879493163
 -0.6792960939508351
 -0.3575966370924379
  2.1577149397345545
 -0.9170456413152324
 -1.021037073006077
                    -1.300754560220488],
              70=>
              [  0.359629402423705
 -0.06242791627632372
 -1.399567921508011
  0.4475847763985716
 -0.8809915030490731
 -0.2017627452986302
 -0.5526516215876736
 -0.35572577945591155
 -0.1607797466665553
  1.1475981747910613
  1.5035230388445053
  0.2812185565541045
 -1.6580378604283104
 -0.5724292571431209
  0.11417065387535505
 -0.1722032235811361
  0.10954150466894817
 -0.20957790746640031
  0.9978114970316212
 -0.33835189120333353
  0.32711196563255096
  1.3252610879493163
 -0.6792960939508351
 -0.3575966370924379
  1.8019224806064675
 -0.1963681193189925
 -0.3211451010184779
  1.2571922200069134
  1.7786053509127615
  0.14664320175311135
 -0.7581751321548498
  0.04258439124049732
  0.3352535968623635
 -0.5458913912138869
 -0.22338451264854198
 -0.34276021437488724
  0.3055002759810978
 -1.7936388425860077
 -1.3858822457781843
  0.19432728408125308
  0.9846027231525627
  0.15797132251418744
 -1.133614480683449
 -1.4271913217697576
  0.681897286305289
 -0.5715336031075239
 -0.4494921999102889
  1.85227626772481
 -0.7005686939516828
  1.2121103292882787
 -0.48418549713360953
  1.2973921022483264
  1.2098214008295365
 -0.9687971894606597
 -1.0648969106849027
  0.935342888589825
  0.8011405444777072
  0.6074278447599545
 -1.3830376414176226
  0.7946559924063759
 -0.9880332785507824
 -1.4752050514434305
 -0.5193534241313768
  2.048090887179751
  2.1577149397345545
 -0.9170456413152324
 -1.021037073006077
 -1.300754560220488
  0.10149275111695119
                 -0.11990415249480861],
              31 => [  0.359629402423705
 -0.06242791627632372
 -1.399567921508011
  0.4475847763985716
 -0.8809915030490731
 -0.2017627452986302
 -0.5526516215876736
 -0.35572577945591155
 -0.1607797466665553
  1.1475981747910613
  1.5035230388445053
  0.2812185565541045
 -1.6580378604283104
 -0.5724292571431209
  0.11417065387535505
 -0.1722032235811361
  0.10954150466894817
 -0.20957790746640031
  0.9978114970316212
 -0.33835189120333353
  0.32711196563255096
  1.3252610879493163
 -0.6792960939508351
 -0.3575966370924379
  2.1577149397345545
 -0.9170456413152324
 -1.021037073006077
 -1.300754560220488
  0.10149275111695119
 -0.11990415249480861
  0.7959522954870297]);

    #x=randn(s);
    #writedlm("/tmp/random_$(f.k)_$s.csv",x);
    x=data[s];
    f.k += s;
    return x;
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
module WebBLAST

  using FastaIO
  using HttpCommon
  using Requests
  import Requests: get, post, put, delete, options
  using LightXML
  using BioSeq

  # include types, web interfaces

  include(joinpath(dirname(@__FILE__), "hit.jl"))
  include(joinpath(dirname(@__FILE__), "ncbi_blast.jl"))
  include(joinpath(dirname(@__FILE__), "ebi_blast.jl"))

  export webblast

  # searches webblast for fasta and returns results as Hit types
  function webblast(sequence::AbstractString, threshold::Float64=0.005, provider::Symbol=:ncbi, cached=false)
      info("Searching for sequence: $(sequence)")

      # try to use cached version
      if cached
        # create cache if not present
        if !isdir(joinpath(dirname(@__FILE__), "..", "..", ".cache"))
          mkdir(joinpath(dirname(@__FILE__), "..", "..", ".cache"))
        end
        # return the cached results
        cached_result_filename = joinpath(dirname(@__FILE__), "..", "..", ".cache", base(16,hash(sequence)))
        if isreadable(cached_result_filename)
          f = open(cached_result_filename, "r")
          return deserialize(f)
        end
      end

      if provider == :ncbi
        rid, rtoe = ncbi_blast_put(sequence)
        info("NCBI rid: $(rid)")
        if ncbi_blast_search_info(rid)
          results = ncbi_blast_get_results(rid, threshold)

          # cache the results
          if cached
            f = open(cached_result_filename, "w")
            serialize(f, results)
          end
          return results
        end
    end
  end

  # webblast with BioSeq AminoAcid Array as input
  function webblast(sequence::Array{AminoAcid,1}, threshold::Float64=0.005, provider::Symbol=:ncbi, cached=false)
    return webblast(convert(ASCIIString, sequence), threshold, provider, cached)
  end

end


text:
@block AnthonyWang [documentation] begin

    @cell "Axis theming" [stepper, axis, lines, stepper, "2d"] begin
        using GeometryBasics
        scene = Scene()
        points = decompose(Point2f0, Circle(Point2f0(10), 10f0), 9)
        lines!(scene, points, linewidth = 8, color = :black)
        axis = scene[Axis] # get axis
        axis[:frame][:linewidth] = 5
        axis[:grid][:linewidth] = (1, 5)
        axis[:grid][:linecolor] = ((:red, 0.3), (:blue, 0.5))
        axis[:names][:axisnames] = ("x", "y   ")
        axis[:ticks][:title_gap] = 1
        axis[:names][:rotation] = (0.0, -3/8*pi)
        axis[:names][:textcolor] = ((:red, 1.0), (:blue, 1.0))
        axis[:ticks][:font] = ("Dejavu Sans", "Helvetica")
        axis[:ticks][:rotation] = (0.0, -pi/2)
        axis[:ticks][:textsize] = (3, 7)
        axis[:ticks][:gap] = 5
        scene
    end

    @cell "Legend" [legend, linesegments, vbox, "2d"] begin
        scene = Scene(resolution = (500, 500))
        x = map([:dot, :dash, :dashdot], [2, 3, 4]) do ls, lw
            linesegments!(
                range(1, stop = 5, length = 100), RNG.rand(100), RNG.rand(100),
                linestyle = ls, linewidth = lw,
                color = RNG.rand(RGBAf0)
            )[end]
        end
        x = [x..., scatter!(range(1, stop=5, length=100), RNG.rand(100), RNG.rand(100))[end]]
        center!(scene)
        ls = AbstractPlotting.legend(x, ["attribute $i" for i in 1:4], camera = campixel!, raw = true)
        l = ls[end]
        l[:strokecolor] = RGBAf0(0.8, 0.8, 0.8)
        l[:gap] = 15
        l[:textsize] = 12
        l[:linepattern] = Point2f0[(0,-0.2), (0.5, 0.2), (0.5, 0.2), (1.0, -0.2)]
        l[:scatterpattern] = decompose(Point2f0, Circle(Point2f0(0.5, 0), 0.3f0), 9)
        l[:markersize] = 2f0
        vbox(scene, ls)
    end

end