text:
# ------------------------------------------------------------------------------------------
# # Delay Coordinates & Neighborhoods
#
#
# Topics:
# * Delay coordinates & why they ROCK!
# * A `Reconstruction` is a subtype of `AbstractDataset`
# * Estimating Parameters for `Reconstruction`s
# * Multiple-time, multiple timeseries `Reconstruction`s
# * Finding neighborhoods of points in a `Dataset`
# * Excluding temporal neighbors
#
# ---
#
# # Delay Coordinates Reconstruction
# Let's say you have a "real-world system" which you measure in an experimental set-up. You
# are assuming that the system is composed of several dynamic variables, but you can only
# measure one of them (or some function of the variable).
#
# You have a severe lack of recorded information for the system. What do you do?
# 1. Give up on science, it is a complete waste of time.
# 2. Use [Taken's theorem](https://en.wikipedia.org/wiki/Takens%27s_theorem), which is
# indistinguishable from magic.
#
# **DynamicalSystems.jl** suggests the third approach.
#
# From a timeseries $s$ one can *reconstruct* a state-space $\mathbf{z}$ simply by shifting
# $s$ in time, like
#
#   $$\mathbf{z}(n) = (s(n), s(n+\tau), s(n+2\tau), \dots, s(n+(D-1)\tau))$$
#
# This is done with the `Reconstruction(s, D, τ)` function
# ------------------------------------------------------------------------------------------

using DynamicalSystems

s = rand(100000)
D = 3 # reconstruction dimension
τ = 4 # reconstruction delay
R = Reconstruction(s, D, τ)

# ------------------------------------------------------------------------------------------
# ---
#
# Here are some nice examples of `Reconstruction`s of a 3D continuous chaotic system, using
# each of the variables of the system, different delay times and dimension of `2`:
# ------------------------------------------------------------------------------------------

using DynamicalSystems, PyPlot

ds = Systems.gissinger(ones(3)) # 3D continuous chaotic system, also shown in orbit diagrams tutorial
dt = 0.05
data = trajectory(ds, 1000.0, dt = dt)

xyz = columns(data)

figure(figsize = (12,10))
k = 1
for i in 1:3
    for τ in [5, 30, 100]
        R = Reconstruction(xyz[i], 2, τ)
        ax = subplot(3,3,k)
        plot(R[:, 1], R[:, 2], color = "C$(k-1)", lw = 0.8)
        title("var = $i, τ = $τ")
        k+=1
    end
end

tight_layout()
suptitle("2D Reconstructions")
subplots_adjust(top=0.9);

# ------------------------------------------------------------------------------------------
# How does this compare to the "real" two-dimensional representation of the system?
# ------------------------------------------------------------------------------------------

figure(figsize=(6,4))
plot(data[:, 2], data[:, 3], lw=1.0);

# ------------------------------------------------------------------------------------------
# ---
#
# # `Reconstruction <: AbstractDataset`
# A `Reconstruction` instance can be passed around and used exactly like a `Dataset`!
#
# Let's look at a `Reconstruction` of data from a gissinger system's trajectory above
# ------------------------------------------------------------------------------------------

R = Reconstruction(data[:, 1], 2, 30)
R[31:end, 1] == R[1:end-30, 2]

a = 0.0
for point ∈ R
    a += mean(point)
end
a/length(R)

# ------------------------------------------------------------------------------------------
# **Taken's theorem says that some quantities remain invariant under a reconstruction**
#
# We'll show this using a `Reconstruction` with the 3rd dimension of `data`.
#
# (D = 2 is not best for this system. D = 3 is better!)
# ------------------------------------------------------------------------------------------

R = Reconstruction(data[:, 1], 3, 30)

# ------------------------------------------------------------------------------------------
# Let's compare the information dimension of the `Reconstruction`
# ------------------------------------------------------------------------------------------

I1 = information_dim(R)

# ------------------------------------------------------------------------------------------
# and the information dimension of the attractor directly
# ------------------------------------------------------------------------------------------

I2 = information_dim(data)

println("|Reconstructed - original| dimension: $(abs(I1 - I2))")

# ------------------------------------------------------------------------------------------
# # Estimating Reconstruction Parameters
#
# It is important to understand that even though Taken's theorem is 99% magic, it is **not**
# 100%. One still has to choose "appropriately good" values for both the delay time as well
# as the reconstruction dimension! Thankfully, **DynamicalSystems.jl** has some support for
# that as well!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# * `estimate_delay` estimates delay time `τ` using the autocorrelation of the signal
# * `estimate_dimension` returns an estimator of the embedding dimension `D` using Cao's
# method
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Multiple-time, multiple-timeseries Reconstructions
#
# The `Reconstruction` we have seen so far is just a "Vanilla version"...
#
# One can also perform:
#
# 1. Reconstructions with multiple delay times, which tries to capture the effect of
# multiple timescales existing in a system.
# 2. Reconstructions with multiple timeseries.
# 3. Reconstructions with multiple timeseries *and* multiple delay times!
#
# See the documentation string of `Reconstruction` for more!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Neighborhoods
#
# A "neighborhood" is a collection of points that is near a given point. `Dataset`s
# interface the module
# [`NearestNeighbors`](https://github.com/KristofferC/NearestNeighbors.jl) in order to find
# this neighborhood.
#
# We use the function `neighborhood`. The call signature is:
# ```julia
# neighborhood(point, tree, ntype)
# ```
#
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# `point` is simply the query point. `tree` is the structure required by
# [`NearestNeighbors`](https://github.com/KristofferC/NearestNeighbors.jl), and is obtained
# simply by:
# ------------------------------------------------------------------------------------------

dataset = Dataset(rand(1000,3))
tree = KDTree(dataset)

# ------------------------------------------------------------------------------------------
# The third argument to `neighborhood` is the *type* of the neighborhood.
#
# * There are two types of neighborhoods!
#
# The first one is defined as the `k` nearest points to a given point. It is represented in
# code by:
# ------------------------------------------------------------------------------------------

mybuddies = FixedMassNeighborhood(3)

# ------------------------------------------------------------------------------------------
# *For experts: for a `FixedMassNeighborhood` a KNN search is done*
# ------------------------------------------------------------------------------------------

point = ones(3)
n = neighborhood(point, tree, mybuddies)

# ------------------------------------------------------------------------------------------
# Notice that the `neighborhood` function does not return the points themselves, but rather
# the indices of the points in the original data:
# ------------------------------------------------------------------------------------------

println("Fixed mass neighborhood of $(point) is:")

for i in n
    println(dataset[i])
end

# ------------------------------------------------------------------------------------------
# ---
#
# The second type of neighborhood contains all the points that are within some given
# distance `ε` from the query.
#
# In code, we represent this as:
# ------------------------------------------------------------------------------------------

where_u_at = FixedSizeNeighborhood(0.001)

# ------------------------------------------------------------------------------------------
#  *For experts: for `FixedSizeNeighborhood` an inrange search is done*
# ------------------------------------------------------------------------------------------

n2 = neighborhood(point, tree, where_u_at)

plz_come_closer = FixedSizeNeighborhood(0.2)
n2 = neighborhood(point, tree, plz_come_closer)

println("Fixed size neighborhood of $(point) is:")

for i in n2
    println(dataset[i])
end

# ------------------------------------------------------------------------------------------
# Okay, so points that have distance < ε are accepted as a neighborhood.
#
# How do we define the "distance" though? When defining a `tree`, you can optionally give a
# distance function. By default Euclidean distance is used, but others also work. For
# example, we can use the `Distances` package to get distance functions,
# ------------------------------------------------------------------------------------------

using Distances

# ------------------------------------------------------------------------------------------
# and define the distance as the `Distances`'s `Chebyshev` distance:
# ------------------------------------------------------------------------------------------

funky_tree = KDTree(dataset, Chebyshev())

n3 = neighborhood(point, funky_tree, plz_come_closer)

# ------------------------------------------------------------------------------------------
# # Excluding temporal neighbors
#
# Before moving on, let's see one last thing.
#
# In this example, the point I want the neighborhood of is now part of my dataset:
# ------------------------------------------------------------------------------------------

point = dataset[end]

# ------------------------------------------------------------------------------------------
# Let's calculate again the two neighborhoods
# ------------------------------------------------------------------------------------------

tree = KDTree(dataset)

# ------------------------------------------------------------------------------------------
# We'll find suuuuuuuper close neighbors with a **very** small $\epsilon$:
# ------------------------------------------------------------------------------------------

ε = 0.000001
where_u_at = FixedSizeNeighborhood(ε)
n2 = neighborhood(point, tree, where_u_at)

# ------------------------------------------------------------------------------------------
# and now we can find the nearest neighbor:
# ------------------------------------------------------------------------------------------

my_best_friend = FixedMassNeighborhood(1)
n3 = neighborhood(point, tree, my_best_friend)

println(n2)
println(n3)

length(dataset) == n2[1] == n3[1]

# ------------------------------------------------------------------------------------------
# **What is happening here is that the `neighborhood` also counted the `point` itself, since
# it is also part of the dataset.**
#
# * Almost always this behavior needs to be avoided. For this reason, there is a second
# method for `neighborhood`:
#
# ```julia
# neighborhood(point, tree, ntype, idx::Int, w::Int = 1)
# ```
#
# In this case, `idx` is the index of the point in the original data. `w` stands for the
# Theiler window (positive integer).
#
# Only points that have index
# `abs(i - idx) ≥ w` are returned as a neighborhood, to exclude close temporal neighbors.
#
# * The default `w=1` is the case of excluding the `point` itself.
#
# ---
#
# Let's revisit the last example (using the default value of `w = 1`):
# ------------------------------------------------------------------------------------------

point = dataset[end]
idx = length(dataset)

n2 = neighborhood(point, tree, where_u_at, idx)
n3 = neighborhood(point, tree, my_best_friend, idx)

println(n2)
println(n3)

# ------------------------------------------------------------------------------------------
# As you can see, there isn't *any* neighbor of `point` with distance `< 0.000001` in this
# dataset, but there is always a nearest neighbor:
# ------------------------------------------------------------------------------------------

println(dataset[n3[1]], " is the nearest neighbor of ", point)

# ------------------------------------------------------------------------------------------
# # Docstrings
# ------------------------------------------------------------------------------------------

?Reconstruction

?neighborhood

?AbstractNeighborhood

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/introduction-to-dynamicalsystems.jl/5. Delay Coordinates & Neighborhoods.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
function _mul!(threading::Threading, C, A, B, sz)
    n, k, m = size(C, 1), size(A, 2), size(C, 2)
    if n >= sz+8 && m >= sz+8 && k >= sz+8
        block_mat_mat_mul!(threading, C, A, B, sz)
    elseif n >= sz+8 && k >= sz+8 && m <  sz+8
        block_mat_vec_mul!(threading, C, A, B, sz)
    elseif n <  sz+8 && k >= sz+8 && m >= sz+8
        block_covec_mat_mul!(threading, C, A, B, sz)
    elseif n >= sz+8 && k <  sz+8 && m >= sz+8
        block_vec_covec_mul!(threading, C, A, B, sz)
    elseif n <  sz+8 && k >= sz+8 && m <  sz+8
        block_covec_vec_mul!(threading, C, A, B, sz)
    else
        gemm_kernel!(C, A, B)
    end
end

function _mul_add!(threading::Threading, C, A, B, sz, ::Val{factor} = Val(1)) where {factor}
    n, k, m = size(C, 1), size(A, 2), size(C, 2)
    if n >= sz+8 && m >= sz+8 && k >= sz+8
        block_mat_mat_mul_add!(threading, C, A, B, sz, Val(factor))
    elseif n >= sz+8 && k >= sz+8 && m <  sz+8
        block_mat_vec_mul_add!(threading, C, A, B, sz, Val(factor))
    elseif n <  sz+8 && k >= sz+8 && m >= sz+8
        block_covec_mat_mul_add!(threading, C, A, B, sz, Val(factor))
    elseif n >= sz+8 && k <  sz+8 && m >= sz+8
        block_vec_covec_mul_add!(threading, C, A, B, sz, Val(factor))
    elseif n <  sz+8 && k >= sz+8 && m <  sz+8
        block_covec_vec_mul_add!(threading, C, A, B, sz, Val(factor))
    else
        add_gemm_kernel!(C, A, B, Val(factor))
    end
end

function _mul!(threading::Threading, C::VecTypes{T}, A::MatTypes{T}, B::VecTypes{T}, sz) where {T<:Eltypes}
    n, k = size(A)
    if     n >= sz+8 && k >= sz+8
        block_mat_vec_mul!(threading, C, A, B, sz)
    elseif n <  sz+8 && k >= sz+8
        block_covec_vec_mul!(threading, C, A, B, sz)
    else
        gemm_kernel!(C, A, B)
    end
end

function _mul_add!(threading::Threading, C::VecTypes{T}, A::MatTypes{T}, B::VecTypes{T}, sz, ::Val{factor} = Val(1)) where {factor, T<:Eltypes}
    n, k = size(A)
    if     n >= sz+8 && k >= sz+8
        block_mat_vec_mul_add!(threading, C, A, B, sz, Val(factor))
    elseif n <  sz+8 && k >= sz+8
        block_covec_vec_mul_add!(threading, C, A, B, sz, Val(factor))
    else
        add_gemm_kernel!(C, A, B, Val(factor))
    end
end

function _mul!(threading::Threading, C::CoVecTypes{T}, A::CoVecTypes{T}, B::MatTypes{T}, sz) where {T<:Eltypes}
    n, k, m = 1, size(A, 1), length(C)
    if     k >= sz+8 && m >= sz+8
        block_covec_mat_mul!(threading, C, A, B, sz)
    elseif k >= sz+8 && m <  sz+8
        block_covec_vec_mul!(threading, C, A, B, sz)
    else
        gemm_kernel!(C, A, B)
    end
end

function _mul_add!(threading::Threading, C::CoVecTypes{T}, A::CoVecTypes{T}, B::MatTypes{T}, sz, ::Val{factor} = Val(1)) where {factor, T<:Eltypes}
    n, k, m = 1, size(A, 1), length(C)
    if     k >= sz+8 && m >= sz+8
        block_covec_mat_mul!(threading, C, A, B, sz, Val(factor))
    elseif k >= sz+8 && m <  sz+8
        block_covec_vec_mul!(threading, C, A, B, sz, Val(factor))
    else
        add_gemm_kernel!(C, A, B, Val(factor))
    end
end

ext: jl


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -5.483106975000851 <= q <= 6.193029789058127)

                     add_NL_constraint(m, :((0.7547513043737739 + -0.48113681936364117*(-0.10262045764995253 + -0.12743000480166566*(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))) + -0.1298782680512347*(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))))/(1 + exp(-(-0.10262045764995253 + -0.12743000480166566*(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))) + -0.1298782680512347*(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))))))) + -0.9820460002458025*(0.1263157800600383 + 0.1285425867505281*(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))) + -0.6128023058992138*(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))))/(1 + exp(-(0.1263157800600383 + 0.1285425867505281*(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))) + -0.6128023058992138*(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))))))))/(1 + exp(-(0.7547513043737739 + -0.48113681936364117*(-0.10262045764995253 + -0.12743000480166566*(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))) + -0.1298782680512347*(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))))/(1 + exp(-(-0.10262045764995253 + -0.12743000480166566*(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))) + -0.1298782680512347*(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))))))) + -0.9820460002458025*(0.1263157800600383 + 0.1285425867505281*(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))) + -0.6128023058992138*(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))))/(1 + exp(-(0.1263157800600383 + 0.1285425867505281*(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))) + -0.6128023058992138*(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))))))))) + (-0.6861294460251406 + 0.36167652828680286*(-0.10262045764995253 + -0.12743000480166566*(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))) + -0.1298782680512347*(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))))/(1 + exp(-(-0.10262045764995253 + -0.12743000480166566*(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))) + -0.1298782680512347*(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))))))) + -0.5440699588803759*(0.1263157800600383 + 0.1285425867505281*(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))) + -0.6128023058992138*(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))))/(1 + exp(-(0.1263157800600383 + 0.1285425867505281*(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))) + -0.6128023058992138*(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))))))))/(1 + exp(-(-0.6861294460251406 + 0.36167652828680286*(-0.10262045764995253 + -0.12743000480166566*(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))) + -0.1298782680512347*(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))))/(1 + exp(-(-0.10262045764995253 + -0.12743000480166566*(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))) + -0.1298782680512347*(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))))))) + -0.5440699588803759*(0.1263157800600383 + 0.1285425867505281*(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))) + -0.6128023058992138*(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))))/(1 + exp(-(0.1263157800600383 + 0.1285425867505281*(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.5778005478095731 + 0.8852363638482865*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + -0.8793054881197468*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))) + -0.6128023058992138*(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2])))))/(1 + exp(-(0.07154777431849713 + -0.5862593464659436*(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2]))/(1 + exp(-(0.4643971142728849 + -0.24024828850470392*$(x[1]) + -0.4967453535429427*$(x[2])))) + 0.9939436563912345*(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))/(1 + exp(-(-0.24055139829732441 + -0.12659057285603792*$(x[1]) + -0.42759160850418754*$(x[2]))))))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
"""
Module with utility functions for BLTandPantelides.

* Author: Hilding Elmqvist, Mogram AB
* Date: July-August 2016
* License: MIT

"""
module BLTandPantelidesUtilities

#using ..BLTandPantelides
using ..ModiaLogging

export buildExtendedSystem, addDependencies, buildFullIncidence
export invertDer, invertAssign
export createNames, printList, printAssignedEquations, printSortedEquations, printUnassigned, makeList

"""
  function buildExtendedSystem(A)

Extend a system according to Pantelides equation (15), i.e. return the incidence for function h(x, der(x)).

* `A`: A[j] = if V[k] = der(V[j]) then k else 0
* `return G`: bipartite graph

Example:
    julia> BLTandPantelidesUtilities.buildExtendedSystem([5,6,7,8,0,0,0,0,0])
    4-element Array{Any,1}:
     [1,5]
     [2,6]
     [3,7]
     [4,8]
"""
function buildExtendedSystem(A)
  G = []
	for i in 1:length(A)
	  a = A[i]
	  if a > 0
      push!(G, [i, a])		# h(x, der(x))
		end
	end
	return G
end

function addDependencies(G, Vindices)
  newG = []
  for g in G
    push!(newG, [g; Vindices])
  end
  return newG
end

"""
    buildFullIncidence(n,m)

Build a bipartite graph with full incidence, i.e. all of the n E-nodes refer to all of the m V-nodes.

* `n`: number of E-nodes
* `m`: number of V-nodes
* `return G`: bipartite graph

Example:
    julia> BLTandPantelidesUtilities.buildFullIncidence(2,3)
    2-element Array{Any,1}:
     [1,2,3]
     [1,2,3]
"""
function buildFullIncidence(n,m)
  G = []
	for i in 1:n
    push!(G, [j for j in 1:m])
	end
	return G
end

"""
    function invertDer(A)

Invert derivative relationships for variables and equations

* `A`: A[j] = if V[k] = der(V[j]) then k else 0 (or correspondingly for E-nodes)
* `return orgIndex`: index of original variable or equation
* `return derOrder`: derivative order

Note that invertDer can be used to invert from list of E-nodes to list of V-nodes as well.

Example:
 julia> BLTandPantelidesUtilities.invertDer([5,6,7,8,10,11,0,0,0,0,0])
([1,2,3,4,1,2,3,4,9,1,2],[0,0,0,0,1,1,1,1,0,2,2])
"""
function invertDer(A)
  orgIndex = [i for i in 1:length(A)]  # Index of original variable or equation
	derOrder = fill(0, length(A)) # Derivative order
  for i in 1:length(A)
	  a = A[i]
		if a > 0
			derOrder[a] = derOrder[i] + 1
			orgIndex[a] = orgIndex[i]
		end
	end
	return orgIndex, derOrder
end

"""
    invertAssign(assign, n=length(assign))

Invert assignment relationships for variables and equations.

* `assign`: assign[j] contains the E-node to which V-node j is assigned or 0 if V-node j not assigned
* `n`: number of E-nodes
* `return invAssign`: invAssign[i] contains the V-node to which E-node i is assigned or 0 if E-node i not assigned
* `return unAssigned`: unassigned V-nodes

Note that invertAssign can be used to invert from list of E-nodes to list of V-nodes as well.

Example:
julia> inv=BLTandPantelidesUtilities.invertAssign([0,0,0,0,1,2,7,4,3,9,8])
([5,6,9,8,0,0,7,11,10,0,0],[1,2,3,4])

julia> BLTandPantelides.invertAssign(inv[1])
([0,0,0,0,1,2,7,4,3,9,8],[5,6,10,11])
"""
function invertAssign(assign, n=length(assign))
  invAssign = fill(0, n)
	unAssigned::Vector{Int} = []
	for j in 1:length(assign)
 	  i = assign[j]
		if i > 0
	  	invAssign[i] = j
		else
		  push!(unAssigned, j)
		end
	end
	return invAssign, unAssigned
end

"""
    function createNames(infixes, A)

Creates names.

* `infixes`: infix strings for original variable
* `A`: A[j] = if V[k] = der(V[j]) then k else 0

Example:
julia> BLTandPantelidesUtilities.createNames(["x", "y", "w", "z", "", "", "", "", "T"], [5,6,7,8,10,11,0,0,0,0,0])
x, y, w, z, der(x), der(y), der(w), der(z), T, der2(x), der2(y)

"""
function createNames(infixes, A)
  names = []
	(orgIndex, derOrder) = invertDer(A)
	for j in 1:length(A)
    if derOrder[j] > 0
      d = "der"
      if derOrder[j] > 1
        d = d*string(derOrder[j])
      end
      d = "$d($(infixes[orgIndex[j]]))"
    else
      d = infixes[orgIndex[j]]
    end
    push!(names, d)
	end
  names
end

"""
    function printList(infixes, indices, A, vertical=false)

Print list of variables or equations.

* `infixes`: infix strings for original variable or equation
* `indices`: indices for the variables or equations to be printed
* `A`: A[j] = if V[k] = der(V[j]) then k else 0 (or correspondingly for E-nodes)
* `vertical`: if vertical then new line separation else comma separation

Example:
julia> BLTandPantelidesUtilities.printList(["x", "y", "w", "z", "", "", "", "", "T"], 1:11, [5,6,7,8,10,11,0,0,0,0,0])
x, y, w, z, der(x), der(y), der(w), der(z), T, der2(x), der2(y)

"""
function printList(infixes, indices, A, vertical=false)
	(orgIndex, derOrder) = invertDer(A)
	for ind in 1:length(indices)
	  j = indices[ind]
		if j > 0
			if ind > 1
			  if vertical
				  loglnModia()
				else
  				logModia(", ")
				end
			end
      if derOrder[j] > 0
				if vertical
					logModia("DER")
				else
					logModia("der")
				end
				if derOrder[j] > 1
					logModia(derOrder[j])
				end
				logModia("(")
			end
			logModia(infixes[orgIndex[j]])
			if derOrder[j] > 0
				logModia(")")
			end
		end
	end
	loglnModia()
end

function makeList(infixes, indices, A, vertical=false)
  l = String[]
	(orgIndex, derOrder) = invertDer(A)
	for ind in 1:length(indices)
    s = ""
    j = indices[ind]
		if j > 0
      if derOrder[j] > 0
				if vertical
					s = "DER"
				else
					s = "der"
				end
				if derOrder[j] > 1
					s *= string(derOrder[j])
				end
				s *= "_" # "("
#        s *= "this."
			end
			s *= string(infixes[orgIndex[j]])
			if derOrder[j] > 0
#				s *= ")"
			end
      push!(l, s)
		end
	end
  return l
end

"""
    printAssignedEquations(equations, variables, indices, assign, A, B)

Print assigned equations.

* `equations`: infix string for original equations
* `variables`: infix string for original variables
* `indices`: indices for the equations to be printed
* `assign`: assign[j] contains the E-node to which V-node j is assigned or 0 if V-node j not assigned
* `A`: A[j] = if V[k] = der(V[j]) then k else 0
* `B`: B[i] = if E[l] = der(E[l]) then l else 0

Example:
See testBLTandPantelides.testPendulum

"""
function printAssignedEquations(equations, variables, indices, assign, A, B)
	(orgIndexVar, derOrderVar) = invertDer(A)
	(orgIndexEqu, derOrderEqu) = invertDer(B)
	(assignedVar, unAssigned) = invertAssign(assign)
  if unAssigned != []
#    @show assignedVar unAssigned
  end
	for i in indices
    logModia(lpad(string(i)*":", 5, " "))
		j = assignedVar[i]
		if j > 0
			if derOrderVar[j] == 1
				prefix = "der("
        suffix = ")"
			elseif derOrderVar[j] > 1
				prefix = "der"*string(derOrderVar[j])*"("
        suffix = ")"
      else
        prefix = ""
        suffix = ""
			end
			logModia(lpad(prefix*string(variables[orgIndexVar[j]])*suffix, 25, " "))
    else
      logModia(" "^25)
		end
		logModia(":  ")

		if derOrderEqu[i] == 1
			prefix = "DER( "
      suffix = " )"
		elseif derOrderEqu[i] > 1
			prefix = "DER"*string(derOrderEqu[i])*"( "
      suffix = " )"
    else
      prefix = ""
      suffix = ""
		end
		loglnModia(prefix*string(equations[orgIndexEqu[i]])*suffix)
	end
end

"""
    printSortedEquations(equations, variables, components, assign, A, B)

Print sorted equations.

* `equations`: infix string for original equations
* `variables`: infix string for original variables
* `components`: cell array of components. Each component is a list of indices to E-nodes
* `assign`: assign[j] contains the E-node to which V-node j is assigned or 0 if V-node j not assigned
* `A`: A[j] = if V[k] = der(V[j]) then k else 0
* `B`: B[i] = if E[l] = der(E[l]) then l else 0

Example:
See testBLTandPantelides.testPendulum

"""
function printSortedEquations(equations, variables, components, assign, A, B)
  loglnModia("[assigned variable]: [differentiation] equation")
  loglnModia("Strongly connected components are enclosed in []")
	for c in components
	  if length(c) > 1
		  loglnModia("[")
		end
		printAssignedEquations(equations, variables, c, assign, A, B)
	  if length(c) > 1
		  loglnModia("]")
		end
	end
end

"""
    printUnassigned(equations, variables, components, assign, A, B)

Print unassigned variables and equations.

* `equations`: infix string for original equations
* `variables`: infix string for original variables
* `assign`: assign[j] contains the E-node to which V-node j is assigned or 0 if V-node j not assigned
* `A`: A[j] = if V[k] = der(V[j]) then k else 0
* `B`: B[i] = if E[l] = der(E[l]) then l else 0

Example:
See testBLTandPantelides.testPendulum

"""
function printUnassigned(equations, variables, assign, A, B, vActive=[])
	(invAssign, unAssignedVariables) = invertAssign(assign, length(B))
	(ass, unAssignedEquations) = invertAssign(invAssign, length(assign))
	if false # BLTandPantelides.log
    @show vActive
  	@show invAssign unAssignedVariables
	  @show ass unAssignedEquations
  end
  if vActive != []
    # Don't print not active variables
    unass = []
    for v in unAssignedVariables
      if vActive[v]
        push!(unass, v)
      end
    end
    unAssignedVariables = unass
  end
  loglnModia("\nUnassigned variables:")
	printList(variables, unAssignedVariables, A)

  loglnModia("\nUnassigned equations:")
	printList(equations, unAssignedEquations, B, true)
end

end


text:
# https://github.com/dcjones/Gadfly.jl

supported_args(::GadflyBackend) = merge_with_base_supported([
    :annotations,
    :background_color, :foreground_color, :color_palette,
    :group, :label, :seriestype,
    :seriescolor, :seriesalpha,
    :linecolor, :linestyle, :linewidth, :linealpha,
    :markershape, :markercolor, :markersize, :markeralpha,
    :markerstrokewidth, :markerstrokecolor, :markerstrokealpha,
    :fillrange, :fillcolor, :fillalpha,
    :bins, :n, :nc, :nr, :layout, :smooth,
    :title, :window_title, :show, :size,
    :x, :xguide, :xlims, :xticks, :xscale, :xflip,
    :y, :yguide, :ylims, :yticks, :yscale, :yflip,
    :z,
    :tickfont, :guidefont, :legendfont,
    :grid, :legend, :colorbar,
    :marker_z, :levels,
    :xerror, :yerror,
    :ribbon, :quiver,
    :orientation,
  ])
supported_types(::GadflyBackend) = [
        :path,
        :scatter, :hexbin,
        :bar,
        :contour, :shape
    ]
supported_styles(::GadflyBackend) = [:auto, :solid, :dash, :dot, :dashdot, :dashdotdot]
supported_markers(::GadflyBackend) = vcat(_allMarkers, Shape)
supported_scales(::GadflyBackend) = [:identity, :ln, :log2, :log10, :asinh, :sqrt]
is_subplot_supported(::GadflyBackend) = true

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

function _initialize_backend(::GadflyBackend; kw...)
    @eval begin
        import Gadfly, Compose
        export Gadfly, Compose
        include(joinpath(Pkg.dir("Plots"), "src", "backends", "gadfly_shapes.jl"))
    end
end

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

# immutable MissingVec <: AbstractVector{Float64} end
# Base.size(v::MissingVec) = (1,)
# Base.getindex(v::MissingVec, i::Integer) = 0.0

function createGadflyPlotObject(d::KW)
    gplt = Gadfly.Plot()
    gplt.mapping = Dict()
    gplt.data_source = Gadfly.DataFrames.DataFrame()
    # gplt.layers = gplt.layers[1:0]
    gplt.layers = [Gadfly.layer(Gadfly.Geom.point(tag=:remove), x=zeros(1), y=zeros(1));] # x=MissingVec(), y=MissingVec());]
    gplt.guides = Gadfly.GuideElement[Gadfly.Guide.xlabel(d[:xguide]),
                                   Gadfly.Guide.ylabel(d[:yguide]),
                                   Gadfly.Guide.title(d[:title])]
    gplt
end

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

function getLineGeom(d::KW)
    st = d[:seriestype]
    xbins, ybins = maketuple(d[:bins])
    if st == :hexb
        Gadfly.Geom.hexbin(xbincount = xbins, ybincount = ybins)
    elseif st == :histogram2d
        Gadfly.Geom.histogram2d(xbincount = xbins, ybincount = ybins)
    elseif st == :histogram
        Gadfly.Geom.histogram(bincount = xbins,
                              orientation = isvertical(d) ? :vertical : :horizontal,
                              position = d[:bar_position] == :stack ? :stack : :dodge)
    elseif st == :path
        Gadfly.Geom.path
    elseif st in (:bar, :sticks)
        Gadfly.Geom.bar
    elseif st == :steppost
        Gadfly.Geom.step
    elseif st == :steppre
        Gadfly.Geom.step(direction = :vh)
    elseif st == :hline
        Gadfly.Geom.hline
    elseif st == :vline
        Gadfly.Geom.vline
    elseif st == :contour
        Gadfly.Geom.contour(levels = d[:levels])
    # elseif st == :shape
    #     Gadfly.Geom.polygon(fill = true, preserve_order = true)
    else
        nothing
    end
end

function get_extra_theme_args(d::KW, k::Symbol)
    # gracefully handles old Gadfly versions
    extra_theme_args = KW()
    try
        extra_theme_args[:line_style] = Gadfly.get_stroke_vector(d[k])
    catch err
        if string(err) == "UndefVarError(:get_stroke_vector)"
            Base.warn_once("Gadfly.get_stroke_vector failed... do you have an old version of Gadfly?")
        else
            rethrow()
        end
    end
    extra_theme_args
end

function getGadflyLineTheme(d::KW)
    st = d[:seriestype]
    lc = convertColor(getColor(d[:linecolor]), d[:linealpha])
    fc = convertColor(getColor(d[:fillcolor]), d[:fillalpha])

    Gadfly.Theme(;
        default_color = (st in (:histogram,:histogram2d,:hexbin,:bar,:sticks) ? fc : lc),
        line_width = (st == :sticks ? 1 : d[:linewidth]) * Gadfly.px,
        # line_style = Gadfly.get_stroke_vector(d[:linestyle]),
        lowlight_color = x->RGB(fc),  # fill/ribbon
        lowlight_opacity = alpha(fc), # fill/ribbon
        bar_highlight = RGB(lc),      # bars
        get_extra_theme_args(d, :linestyle)...
    )
end

# add a line as a new layer
function addGadflyLine!(plt::Plot, numlayers::Int, d::KW, geoms...)
    gplt = getGadflyContext(plt)
    gfargs = vcat(geoms..., getGadflyLineTheme(d))
    kwargs = KW()
    st = d[:seriestype]

    # add a fill?
    if d[:fillrange] != nothing && st != :contour
        fillmin, fillmax = map(makevec, maketuple(d[:fillrange]))
        nmin, nmax = length(fillmin), length(fillmax)
        kwargs[:ymin] = Float64[min(y, fillmin[mod1(i, nmin)], fillmax[mod1(i, nmax)]) for (i,y) in enumerate(d[:y])]
        kwargs[:ymax] = Float64[max(y, fillmin[mod1(i, nmin)], fillmax[mod1(i, nmax)]) for (i,y) in enumerate(d[:y])]
        push!(gfargs, Gadfly.Geom.ribbon)
    end

    if st in (:hline, :vline)
        kwargs[st == :hline ? :yintercept : :xintercept] = d[:y]

    else
        if st == :sticks
            w = 0.01 * mean(diff(d[:x]))
            kwargs[:xmin] = d[:x] - w
            kwargs[:xmax] = d[:x] + w
        elseif st == :contour
            kwargs[:z] = d[:z].surf
            addGadflyContColorScale(plt, d[:linecolor])
        end

        kwargs[:x] = d[st == :histogram ? :y : :x]
        kwargs[:y] = d[:y]

    end

    # # add the layer
    Gadfly.layer(gfargs...; order=numlayers, kwargs...)
end

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

get_shape(sym::Symbol) = _shapes[sym]
get_shape(shape::Shape) = shape

# extract the underlying ShapeGeometry object(s)
getMarkerGeom(shapes::AVec) = gadflyshape(map(get_shape, shapes))
getMarkerGeom(other) = gadflyshape(get_shape(other))

# getMarkerGeom(shape::Shape) = gadflyshape(shape)
# getMarkerGeom(shape::Symbol) = gadflyshape(_shapes[shape])
# getMarkerGeom(shapes::AVec) = gadflyshape(map(gadflyshape, shapes)) # map(getMarkerGeom, shapes)
function getMarkerGeom(d::KW)
    if d[:seriestype] == :shape
        Gadfly.Geom.polygon(fill = true, preserve_order = true)
    else
        getMarkerGeom(d[:markershape])
    end
end

function getGadflyMarkerTheme(d::KW, attr::KW)
    c = getColor(d[:markercolor])
    α = d[:markeralpha]
    if α != nothing
        c = RGBA(RGB(c), α)
    end

    ms = d[:markersize]
    ms = if typeof(ms) <: AVec
        warn("Gadfly doesn't support variable marker sizes... using the average: $(mean(ms))")
        mean(ms) * Gadfly.px
    else
        ms * Gadfly.px
    end

    Gadfly.Theme(;
        default_color = c,
        default_point_size = ms,
        discrete_highlight_color = c -> RGB(getColor(d[:markerstrokecolor])),
        highlight_width = d[:markerstrokewidth] * Gadfly.px,
        line_width = d[:markerstrokewidth] * Gadfly.px,
        # get_extra_theme_args(d, :markerstrokestyle)...
    )
end

function addGadflyContColorScale(plt::Plot{GadflyBackend}, c)
    plt.attr[:colorbar] == :none && return
    if !isa(c, ColorGradient)
        c = default_gradient()
    end
    push!(getGadflyContext(plt).scales, Gadfly.Scale.ContinuousColorScale(p -> RGB(getColorZ(c, p))))
end

function addGadflyMarker!(plt::Plot, numlayers::Int, d::KW, attr::KW, geoms...)
    gfargs = vcat(geoms..., getGadflyMarkerTheme(d, attr), getMarkerGeom(d))
    kwargs = KW()

    # handle continuous color scales for the markers
    zcolor = d[:marker_z]
    if zcolor != nothing && typeof(zcolor) <: AVec
        kwargs[:color] = zcolor
        addGadflyContColorScale(plt, d[:markercolor])
    end

    Gadfly.layer(gfargs...; x = d[:x], y = d[:y], order=numlayers, kwargs...)
end

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

function addToGadflyLegend(plt::Plot, d::KW)
    if plt.attr[:legend] != :none && d[:label] != ""
        gplt = getGadflyContext(plt)

        # add the legend if needed
        if all(g -> !isa(g, Gadfly.Guide.ManualColorKey), gplt.guides)
            unshift!(gplt.guides, Gadfly.Guide.manual_color_key("", @compat(AbstractString)[], Color[]))
        end

        # now add the series to the legend
        for guide in gplt.guides
            if isa(guide, Gadfly.Guide.ManualColorKey)
                # TODO: there's a BUG in gadfly if you pass in the same color more than once,
                # since gadfly will call unique(colors), but doesn't also merge the rows that match
                # Should ensure from this side that colors which are the same are merged together

                c = getColor(d[d[:markershape] == :none ? :linecolor : :markercolor])
                foundit = false

                # extend the label if we found this color
                for i in 1:length(guide.colors)
                    if RGB(c) == guide.colors[i]
                        guide.labels[i] *= ", " * d[:label]
                        foundit = true
                    end
                end

                # didn't find the color, so add a new entry into the legend
                if !foundit
                    push!(guide.labels, d[:label])
                    push!(guide.colors, c)
                end
            end
        end
    end
end

getGadflySmoothing(smooth::Bool) = smooth ? [Gadfly.Geom.smooth(method=:lm)] : Any[]
getGadflySmoothing(smooth::Real) = [Gadfly.Geom.smooth(method=:loess, smoothing=float(smooth))]

function addGadflySeries!(plt::Plot, d::KW)
    layers = Gadfly.Layer[]
    gplt = getGadflyContext(plt)

    # add a regression line?
    # TODO: make more flexible
    smooth = getGadflySmoothing(d[:smooth])

    # lines
    geom = getLineGeom(d)
    if geom != nothing
        prepend!(layers, addGadflyLine!(plt, length(gplt.layers), d, geom, smooth...))
        smooth = Any[] # don't add a regression for markers too
    end

    # special handling for ohlc and scatter
    st = d[:seriestype]
    # if st == :ohlc
    #     error("Haven't re-implemented after refactoring")
    if st in (:histogram2d, :hexbin) && (isa(d[:fillcolor], ColorGradient) || isa(d[:fillcolor], ColorFunction))
        push!(gplt.scales, Gadfly.Scale.ContinuousColorScale(p -> RGB(getColorZ(d[:fillcolor], p))))
    elseif st == :scatter && d[:markershape] == :none
        d[:markershape] = :circle
    end

    # markers
    if d[:markershape] != :none || st == :shape
        prepend!(layers, addGadflyMarker!(plt, length(gplt.layers), d, plt.attr, smooth...))
    end

    st in (:histogram2d, :hexbin, :contour) || addToGadflyLegend(plt, d)

    # now save the layers that apply to this series
    d[:gadflylayers] = layers
    prepend!(gplt.layers, layers)
end

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

# NOTE: I'm leaving this here and commented out just in case I want to implement again... it was hacky code to create multi-colored line segments

#   # colorgroup
#   z = d[:z]

#   # handle line segments of different colors
#   cscheme = d[:linecolor]
#   if isa(cscheme, ColorVector)
#     # create a color scale, and set the color group to the index of the color
#     push!(gplt.scales, Gadfly.Scale.color_discrete_manual(cscheme.v...))

#     # this is super weird, but... oh well... for some reason this creates n separate line segments...
#     # create a list of vertices that go: [x1,x2,x2,x3,x3, ... ,xi,xi, ... xn,xn] (same for y)
#     # then the vector passed to the "color" keyword should be a vector: [1,1,2,2,3,3,4,4, ..., i,i, ... , n,n]
#     csindices = Int[mod1(i,length(cscheme.v)) for i in 1:length(d[:y])]
#     cs = collect(repmat(csindices', 2, 1))[1:end-1]
#     grp = collect(repmat((1:length(d[:y]))', 2, 1))[1:end-1]
#     d[:x], d[:y] = map(createSegments, (d[:x], d[:y]))
#     colorgroup = [(:linecolor, cs), (:group, grp)]

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

function addGadflyTicksGuide(gplt, ticks, isx::Bool)
    ticks == :auto && return

    # remove the ticks?
    if ticks in (:none, false, nothing)
        return addOrReplace(gplt.guides, isx ? Gadfly.Guide.xticks : Gadfly.Guide.yticks; label=false)
    end

    ttype = ticksType(ticks)

    # just the values... put ticks here, but use standard labels
    if ttype == :ticks
        gtype = isx ? Gadfly.Guide.xticks : Gadfly.Guide.yticks
        replaceType(gplt.guides, gtype(ticks = collect(ticks)))

    # set the ticks and the labels
    # Note: this is pretty convoluted, but I think it works.  We set the ticks using Gadfly.Guide,
    #   and then set the label function (wraps a dict lookup) through a continuous Gadfly.Scale.
    elseif ttype == :ticks_and_labels
        gtype = isx ? Gadfly.Guide.xticks : Gadfly.Guide.yticks
        replaceType(gplt.guides, gtype(ticks = collect(ticks[1])))

        # # TODO add xtick_label function (given tick, return label??)
        # # Scale.x_discrete(; labels=nothing, levels=nothing, order=nothing)
        # filterGadflyScale(gplt, isx)
        # gfunc = isx ? Gadfly.Scale.x_discrete : Gadfly.Scale.y_discrete
        # labelmap = Dict(zip(ticks...))
        # labelfunc = val -> labelmap[val]
        # push!(gplt.scales, gfunc(levels = collect(ticks[1]), labels = labelfunc))

        filterGadflyScale(gplt, isx)
        gfunc = isx ? Gadfly.Scale.x_continuous : Gadfly.Scale.y_continuous
        labelmap = Dict(zip(ticks...))
        labelfunc = val -> labelmap[val]
        push!(gplt.scales, gfunc(labels = labelfunc))

    else
        error("Invalid input for $(isx ? "xticks" : "yticks"): ", ticks)
    end
end

continuousAndSameAxis(scale, isx::Bool) = isa(scale, Gadfly.Scale.ContinuousScale) && scale.vars[1] == (isx ? :x : :y)
filterGadflyScale(gplt, isx::Bool) = filter!(scale -> !continuousAndSameAxis(scale, isx), gplt.scales)

function getGadflyScaleFunction(d::KW, isx::Bool)
    scalekey = isx ? :xscale : :yscale
    hasScaleKey = haskey(d, scalekey)
    if hasScaleKey
        scale = d[scalekey]
        scale == :ln && return isx ? Gadfly.Scale.x_log : Gadfly.Scale.y_log, hasScaleKey, log
        scale == :log2 && return isx ? Gadfly.Scale.x_log2 : Gadfly.Scale.y_log2, hasScaleKey, log2
        scale == :log10 && return isx ? Gadfly.Scale.x_log10 : Gadfly.Scale.y_log10, hasScaleKey, log10
        scale == :asinh && return isx ? Gadfly.Scale.x_asinh : Gadfly.Scale.y_asinh, hasScaleKey, asinh
        scale == :sqrt && return isx ? Gadfly.Scale.x_sqrt : Gadfly.Scale.y_sqrt, hasScaleKey, sqrt
    end
    isx ? Gadfly.Scale.x_continuous : Gadfly.Scale.y_continuous, hasScaleKey, identity
end

function addGadflyLimitsScale(gplt, d::KW, isx::Bool)
    gfunc, hasScaleKey, func = getGadflyScaleFunction(d, isx)

    # do we want to add min/max limits for the axis?
    limsym = isx ? :xlims : :ylims
    limargs = Any[]

    # map :auto to nothing, otherwise add to limargs
    lims = get(d, limsym, :auto)
    if lims == :auto
        lims = nothing
    else
        if limsType(lims) == :limits
            push!(limargs, (:minvalue, min(lims...)))
            push!(limargs, (:maxvalue, max(lims...)))
        else
            error("Invalid input for $(isx ? "xlims" : "ylims"): ", lims)
        end
    end

    # replace any current scales with this one
    if hasScaleKey || !isempty(limargs)
        filterGadflyScale(gplt, isx)
        push!(gplt.scales, gfunc(; limargs...))
    end

    lims, func
end

function updateGadflyAxisFlips(gplt, d::KW, xlims, ylims, xfunc, yfunc)
    if isa(gplt.coord, Gadfly.Coord.Cartesian)
        gplt.coord = Gadfly.Coord.cartesian(
            gplt.coord.xvars,
            gplt.coord.yvars;
            xmin = xlims == nothing ? gplt.coord.xmin : xfunc(minimum(xlims)),
            xmax = xlims == nothing ? gplt.coord.xmax : xfunc(maximum(xlims)),
            ymin = ylims == nothing ? gplt.coord.ymin : yfunc(minimum(ylims)),
            ymax = ylims == nothing ? gplt.coord.ymax : yfunc(maximum(ylims)),
            xflip = get(d, :xflip, gplt.coord.xflip),
            yflip = get(d, :yflip, gplt.coord.yflip),
            fixed = gplt.coord.fixed,
            aspect_ratio = gplt.coord.aspect_ratio,
            raster = gplt.coord.raster
        )
    else
        gplt.coord = Gadfly.Coord.Cartesian(
            xflip = get(d, :xflip, false),
            yflip = get(d, :yflip, false)
        )
    end
end

function findGuideAndSet(gplt, t::DataType, args...; kw...) #s::@compat(AbstractString))
    for (i,guide) in enumerate(gplt.guides)
        if isa(guide, t)
            gplt.guides[i] = t(args...; kw...)
        end
    end
end

function updateGadflyGuides(plt::Plot, d::KW)
    gplt = getGadflyContext(plt)
    haskey(d, :title) && findGuideAndSet(gplt, Gadfly.Guide.title, string(d[:title]))
    haskey(d, :xguide) && findGuideAndSet(gplt, Gadfly.Guide.xlabel, string(d[:xguide]))
    haskey(d, :yguide) && findGuideAndSet(gplt, Gadfly.Guide.ylabel, string(d[:yguide]))

    xlims, xfunc = addGadflyLimitsScale(gplt, d, true)
    ylims, yfunc = addGadflyLimitsScale(gplt, d, false)

    ticks = get(d, :xticks, :auto)
    if ticks == :none
        _remove_axis(plt, true)
    else
        addGadflyTicksGuide(gplt, ticks, true)
    end
    ticks = get(d, :yticks, :auto)
    if ticks == :none
        _remove_axis(plt, false)
    else
        addGadflyTicksGuide(gplt, ticks, false)
    end

    updateGadflyAxisFlips(gplt, d, xlims, ylims, xfunc, yfunc)
end

function updateGadflyPlotTheme(plt::Plot, d::KW)
    kwargs = KW()

    # colors
    insidecolor, gridcolor, textcolor, guidecolor, legendcolor =
        map(s -> getColor(d[s]), (
            :background_color_inside,
            :foreground_color_grid,
            :foreground_color_text,
            :foreground_color_guide,
            :foreground_color_legend
        ))

    # # hide the legend?
    leg = d[d[:legend] == :none ? :colorbar : :legend]
    if leg != :best
        kwargs[:key_position] = leg == :inside ? :right : leg
    end

    if !get(d, :grid, true)
        kwargs[:grid_color] = gridcolor
    end

    # fonts
    tfont, gfont, lfont = d[:tickfont], d[:guidefont], d[:legendfont]

    getGadflyContext(plt).theme = Gadfly.Theme(;
        background_color = insidecolor,
        minor_label_color = textcolor,
        minor_label_font = tfont.family,
        minor_label_font_size = tfont.pointsize * Gadfly.pt,
        major_label_color = guidecolor,
        major_label_font = gfont.family,
        major_label_font_size = gfont.pointsize * Gadfly.pt,
        key_title_color = guidecolor,
        key_title_font = gfont.family,
        key_title_font_size = gfont.pointsize * Gadfly.pt,
        key_label_color = legendcolor,
        key_label_font = lfont.family,
        key_label_font_size = lfont.pointsize * Gadfly.pt,
        plot_padding = 1 * Gadfly.mm,
        kwargs...
    )
end

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

function createGadflyAnnotationObject(x, y, val::@compat(AbstractString))
    Gadfly.Guide.annotation(Compose.compose(
                              Compose.context(),
                              Compose.text(x, y, val)
                            ))
end

function createGadflyAnnotationObject(x, y, txt::PlotText)
    halign = (txt.font.halign == :hcenter ? Compose.hcenter : (txt.font.halign == :left ? Compose.hleft : Compose.hright))
    valign = (txt.font.valign == :vcenter ? Compose.vcenter : (txt.font.valign == :top ? Compose.vtop : Compose.vbottom))
    rotations = (txt.font.rotation == 0.0 ? [] : [Compose.Rotation(txt.font.rotation, Compose.Point(Compose.x_measure(x), Compose.y_measure(y)))])
    Gadfly.Guide.annotation(Compose.compose(
                              Compose.context(),
                              Compose.text(x, y, txt.str, halign, valign, rotations...),
                              Compose.font(string(txt.font.family)),
                              Compose.fontsize(txt.font.pointsize * Gadfly.pt),
                              Compose.stroke(txt.font.color),
                              Compose.fill(txt.font.color)
                            ))
end

function _add_annotations{X,Y,V}(plt::Plot{GadflyBackend}, anns::AVec{@compat(Tuple{X,Y,V})})
    for ann in anns
        push!(plt.o.guides, createGadflyAnnotationObject(ann...))
    end
end

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

# create a blank Gadfly.Plot object
# function _create_plot(pkg::GadflyBackend, d::KW)
#     gplt = createGadflyPlotObject(d)
#     Plot(gplt, pkg, 0, d, KW[])
# end
function _create_backend_figure(plt::Plot{GadflyBackend})
    createGadflyPlotObject(plt.attr)
end

# plot one data series
# function _series_added(::GadflyBackend, plt::Plot, d::KW)
function _series_added(plt::Plot{GadflyBackend}, series::Series)
    # first clear out the temporary layer
    gplt = getGadflyContext(plt)
    if gplt.layers[1].geom.tag == :remove
        gplt.layers = gplt.layers[2:end]
    end

    addGadflySeries!(plt, series.d)
    # push!(plt.seriesargs, d)
    # plt
end

function _update_plot_object(plt::Plot{GadflyBackend}, d::KW)
    updateGadflyGuides(plt, d)
    updateGadflyPlotTheme(plt, d)
end

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

# accessors for x/y data

# TODO: need to save all the layer indices which apply to this series
function getGadflyMappings(plt::Plot, i::Integer)
    @assert i > 0 && i <= plt.n
    mappings = [l.mapping for l in plt.seriesargs[i][:gadflylayers]]
end

function getxy(plt::Plot{GadflyBackend}, i::Integer)
    mapping = getGadflyMappings(plt, i)[1]
    mapping[:x], mapping[:y]
end

function setxy!{X,Y}(plt::Plot{GadflyBackend}, xy::Tuple{X,Y}, i::Integer)
    for mapping in getGadflyMappings(plt, i)
        mapping[:x], mapping[:y] = xy
    end
    plt
end

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

# # create the underlying object (each backend will do this differently)
# function _create_subplot(subplt::Subplot{GadflyBackend}, isbefore::Bool)
#     isbefore && return false # wait until after plotting to create the subplots
#     subplt.o = nothing
#     true
# end

function _remove_axis(plt::Plot{GadflyBackend}, isx::Bool)
    gplt = getGadflyContext(plt)
    addOrReplace(gplt.guides, isx ? Gadfly.Guide.xticks : Gadfly.Guide.yticks; label=false)
    addOrReplace(gplt.guides, isx ? Gadfly.Guide.xlabel : Gadfly.Guide.ylabel, "")
end

function _expand_limits(lims, plt::Plot{GadflyBackend}, isx::Bool)
    for l in getGadflyContext(plt).layers
        _expand_limits(lims, l.mapping[isx ? :x : :y])
    end
end

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

getGadflyContext(plt::Plot{GadflyBackend}) = plt.o
# getGadflyContext(subplt::Subplot{GadflyBackend}) = buildGadflySubplotContext(subplt)

# # create my Compose.Context grid by hstacking and vstacking the Gadfly.Plot objects
# function buildGadflySubplotContext(subplt::Subplot)
#     rows = Any[]
#     row = Any[]
#     for (i,(r,c)) in enumerate(subplt.layout)
#
#         # add the Plot object to the row
#         push!(row, getGadflyContext(subplt.plts[i]))
#
#         # add the row
#         if c == ncols(subplt.layout, r)
#             push!(rows, Gadfly.hstack(row...))
#             row = Any[]
#         end
#     end
#
#     # stack the rows
#     Gadfly.vstack(rows...)
# end

setGadflyDisplaySize(w,h) = Compose.set_default_graphic_size(w * Compose.px, h * Compose.px)
setGadflyDisplaySize(plt::Plot) = setGadflyDisplaySize(plt.attr[:size]...)
# setGadflyDisplaySize(subplt::Subplot) = setGadflyDisplaySize(getattr(subplt, 1)[:size]...)
# -------------------------------------------------------------------------

function dowritemime{P<:Union{GadflyBackend,ImmerseBackend}}(io::IO, func, plt::AbstractPlot{P})
    gplt = getGadflyContext(plt)
    setGadflyDisplaySize(plt)
    Gadfly.draw(func(io, Compose.default_graphic_width, Compose.default_graphic_height), gplt)
end

getGadflyWriteFunc(::MIME"image/png") = Gadfly.PNG
getGadflyWriteFunc(::MIME"image/svg+xml") = Gadfly.SVG
# getGadflyWriteFunc(::MIME"text/html") = Gadfly.SVGJS
getGadflyWriteFunc(::MIME"application/pdf") = Gadfly.PDF
getGadflyWriteFunc(::MIME"application/postscript") = Gadfly.PS
getGadflyWriteFunc(::MIME"application/x-tex") = Gadfly.PGF
getGadflyWriteFunc(m::MIME) = error("Unsupported in Gadfly/Immerse: ", m)

for mime in (MIME"image/png", MIME"image/svg+xml", MIME"application/pdf", MIME"application/postscript", MIME"application/x-tex")
    @eval function Base.writemime{P<:Union{GadflyBackend,ImmerseBackend}}(io::IO, ::$mime, plt::AbstractPlot{P})
        func = getGadflyWriteFunc($mime())
        dowritemime(io, func, plt)
    end
end

function Base.display(::PlotsDisplay, plt::Plot{GadflyBackend})
    setGadflyDisplaySize(plt.attr[:size]...)
    display(plt.o)
end

# function Base.display(::PlotsDisplay, subplt::Subplot{GadflyBackend})
#     setGadflyDisplaySize(getattr(subplt,1)[:size]...)
#     ctx = buildGadflySubplotContext(subplt)
#
#     # taken from Gadfly since I couldn't figure out how to do it directly
#
#     filename = string(Gadfly.tempname(), ".html")
#     output = open(filename, "w")
#
#     plot_output = IOBuffer()
#     Gadfly.draw(Gadfly.SVGJS(plot_output, Compose.default_graphic_width,
#              Compose.default_graphic_height, false), ctx)
#     plotsvg = takebuf_string(plot_output)
#
#     write(output,
#       """
#       <!DOCTYPE html>
#       <html>
#         <head>
#           <title>Gadfly Plot</title>
#           <meta charset="utf-8">
#         </head>
#           <body>
#           <script charset="utf-8">
#               $(readall(Compose.snapsvgjs))
#           </script>
#           <script charset="utf-8">
#               $(readall(Gadfly.gadflyjs))
#           </script>
#           $(plotsvg)
#         </body>
#       </html>
#       """)
#     close(output)
#     Gadfly.open_file(filename)
# end


text:
using Metida
using SnoopCompile
path    = dirname(@__FILE__)
cd(path)

inf_timing = @snoopi tmin=0.00001 include("test.jl")
pc = SnoopCompile.parcel(inf_timing)
SnoopCompile.write("precompile", pc)

SnoopCompile.@snoopc "$path/precompile/metida_compiles.log" begin
    using Metida, Pkg
    include(joinpath(dirname(dirname(pathof(Metida))), "test", "test.jl"))
end
data = SnoopCompile.read("$path/precompile/metida_compiles.log")
pc = SnoopCompile.parcel(reverse!(data[2]))
SnoopCompile.write("$path/precompile", pc)

using SnoopCompileCore
invalidations = @snoopr begin
    using Metida
end
using SnoopCompile
ui = uinvalidated(invalidations)
trees = invalidation_trees(invalidations)
ftrees = filtermod(Metida, trees)


text:
using SwiftObjectStores
using Base.Test

# write your own tests here
include("integration_tests.jl")


text:
# %% Fractal dimension of experimental datasets
using DrWatson
@quickactivate :FractalDimension # uses DynamicalSystems, PyPlot
include(srcdir("style.jl"))

names = String[
    "electrochemical1", "electrochemical2", "shinriki1",
    "roessler_embed",  "doublependulum",
    # "roessler_all", # doesn't make sense to use 28 dimensional data...
]
labels = String[
    "electroch. 1", "electroch. 2", "Shinriki",
    "Rössler Net",  "Double pend.",
    # "Rössler all",
]

all_ros_range = MathConstants.e .^ (-4:0.25:1)
Cmethod = "standard" # bueno or standard. Decides εmax for correlation sum.

qH = 1
qC = 2
# Calculate values for H, C
eHs, eCs, Hs, Cs = [Vector{Float64}[] for i in 1:4]
for i in 1:length(labels)
    # These are the parameters that change between the different elements of the plot
    data = :experimental_data
    name = names[i]

    # Here we simply pack all parameters into a dictionary
    # (other parameters are (probably) globals)
    params = @strdict qH qC data name
    if Cmethod ≠ "standard"
        params["Cmethod"] = Cmethod
    end
    params["z"] = -1

    # This is the main call that calculates everything
    output, s = produce_or_load(
        datadir("main"), params, make_C_H;
        prefix = string(data), suffix = "jld2", force = false,
        ignores = ["data"], storepatch = false,
    )
    @unpack eH, eC, H, C = output
    push!(eHs, eH); push!(Hs, H); push!(eCs, eC); push!(Cs, C)
end

# Do the actual plot
legendtitle = "experimental data"

fig, axs = mainplot(
    Hs, Cs, eHs, eCs, labels, legendtitle;
    qH, qC, tol = 0.25,
    offsets = range(0; length = 6, step = 2),
    # dimension_fit_C = FractalDimension.logarithmic_corrected_fit_lsqfit,
    dimension_fit_C = FractalDimension.linear_regression_fit_glm,
)
wsave(plotsdir("paper", "experimental"), fig)

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
using Test
using DataStructure101
const DS = DataStructure101

@testset "Priority Queue" begin
    pq = DS.PriorityQueue{Char}()
    for i = 1:26
        push!(pq, '`'+i)
    end
    @test repr(pq) == "PriorityQueue{Char}([" *
        "'a', 'b', 'c', 'd', 'e', 'f', 'g', " *
        "'h', 'i', 'j', 'k', 'l', 'm', 'n', " *
        "'o', 'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z']"
    @test minimum(pq) == (1 => 'a')
    @test maximum(pq) == (26 => 'z')
    @test DS.pop_min!(pq) == 'a'
    @test DS.pop_max!(pq) == 'z'
    DS.update!(pq, 24, -5)
    @test pq.data[end-5] == 'y'
    DS.update!(pq, 1, 5)
    @test pq.data[6] == 'b'
end


text:
using RecipesBase, Measures
export InvMcmcSampler, McmcOutput

struct InvMcmcSampler{T<:InvMcmcMove}
    move::T
    desired_samples::Int
    burn_in::Int
    lag::Int
    K::Int
    init::McmcInitialiser
    pointers::InteractionSequence{Int}
    function InvMcmcSampler(
        move::T;
        desired_samples::Int=1000, burn_in::Int=0, lag::Int=1,
        K::Int=100, init=InitMode()
        ) where {T<:InvMcmcMove}
        pointers = [Int[] for i in 1:(2K)]
        new{T}(move, desired_samples, burn_in, lag, K, init, pointers)
    end
end

Base.show(io::IO, x::InvMcmcSampler{T}) where {T} = print(io, typeof(x))

struct McmcOutput{T<:Union{SIS,SIM}}
    sample::Vector{Vector{Path{Int}}}  # The sample
    model::T
end

Base.show(io::IO, x::McmcOutput) = print(io, typeof(x))

@recipe function f(output::McmcOutput)
    model = output.model
    sample = output.sample
    x = map(x -> model.dist(model.mode, x), sample)
    xguide --> "Sample"
    yguide --> "Distance from Mode"
    legend --> false
    size --> (800, 300)
    margin --> 5mm
    x
end

acceptance_prob(mcmc::InvMcmcSampler) = acceptance_prob(mcmc.move)

function eval_accept_prob(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    model::T,
    log_ratio::Float64
    ) where {T<:Union{SIS,SIM}}

    mode, γ, dist = (model.mode, model.γ, model.dist)

    # @show S_curr, S_prop
    log_lik_ratio = -γ * (
        dist(mode, S_prop)-dist(mode, S_curr)
    )

    if typeof(model)==SIM
        log_multinom_term = log_multinomial_ratio(S_curr, S_prop)
        return log_lik_ratio + log_ratio + log_multinom_term
    else
        return log_lik_ratio + log_ratio
    end

end

function accept_reject!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    pointers::InteractionSequence{Int},
    move::InvMcmcMove,
    model::T
    ) where {T<:Union{SIS,SIM}}

    move.counts[2] += 1

    log_ratio = prop_sample!(S_curr, S_prop, move, pointers, model.V)

    # Catch out of bounds proposals (reject them, i.e. 0 acc prob)
    if any(!(1 ≤ length(x) ≤ model.K_inner.u) for x in S_prop)
        log_α = -Inf
    elseif !(1 ≤ length(S_prop) ≤ model.K_outer.u)
        log_α = -Inf
    else
        log_α = eval_accept_prob(S_curr, S_prop, model, log_ratio)
    end

    # @show log_α
    if log(rand()) < log_α
        # We accept!
        move.counts[1] += 1
        enact_accept!(S_curr, S_prop, pointers, move)
    else
        # We reject!
        enact_reject!(S_curr, S_prop, pointers, move)
    end

end

function intialise_states!(
    pointers::InteractionSequence{Int},
    init::InteractionSequence{Int}
    )

    S_curr = InteractionSequence{Int}()
    S_prop = InteractionSequence{Int}()
    for init_path in init
        tmp = pop!(pointers)
        copy!(tmp, init_path)
        push!(S_curr, tmp)
        tmp = pop!(pointers)
        copy!(tmp, init_path)
        push!(S_prop, tmp)
    end

    return S_curr, S_prop

end

function return_states!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    pointers::InteractionSequence{Int}
    )
    for i in eachindex(S_curr)
        tmp = pop!(S_curr)
        pushfirst!(pointers, tmp)
        tmp = pop!(S_prop)
        pushfirst!(pointers, tmp)
    end
end

function draw_sample!(
    sample_out::Union{InteractionSequenceSample{Int}, SubArray},
    mcmc::InvMcmcSampler,
    model::T;
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    init::InteractionSequence{Int}=get_init(mcmc.init, model)
    ) where {T<:Union{SIS,SIM}}

    pointers = mcmc.pointers

    S_curr, S_prop = intialise_states!(pointers, init)

    sample_count = 1    # Keeps which sample to be stored we are working to get
    i = 0               # Keeps track all samples (included lags and burn_ins)
    reset_counts!(mcmc.move) # Reset counts for mcmc move (tracks acceptances)

    while sample_count ≤ length(sample_out)
        i += 1
        # Store value
        if (i > burn_in) & (((i-1) % lag)==0)
            @inbounds sample_out[sample_count] = deepcopy(S_curr)
            sample_count += 1
        end
        accept_reject!(
            S_curr, S_prop,
            pointers,
            mcmc.move,
            model
        )
        # println(S_prop)
    end
    return_states!(S_curr, S_prop, pointers)
end

function draw_sample(
    mcmc::InvMcmcSampler,
    model::T;
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    init::InteractionSequence{Int}=get_init(mcmc.init, model)
    ) where {T<:Union{SIS,SIM}}

    sample_out = InteractionSequenceSample{Int}(undef, desired_samples)
    draw_sample!(sample_out, mcmc, model, burn_in=burn_in, lag=lag, init=init)
    return sample_out
end

function (mcmc::InvMcmcSampler)(
    model::T;
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    init::InteractionSequence{Int}=get_init(mcmc.init, model)
    ) where {T<:Union{SIS,SIM}}

    sample_out = draw_sample(mcmc, model, desired_samples=desired_samples, burn_in=burn_in, lag=lag, init=init)

    return McmcOutput(sample_out, model)

end


text:
const pyjlvectortype = pynew()

function pyjlvector_resize(x::AbstractVector, size_::Py)
    size = pyconvertarg(Int, size_, "size")
    resize!(x, size)
    Py(nothing)
end

function pyjlvector_sort(x::AbstractVector, reverse_::Py, key_::Py)
    reverse = pyconvertarg(Bool, reverse_, "reverse")
    key = pyconvertarg(Any, key_, "size")
    if key === nothing
        sort!(x, rev=reverse)
    else
        sort!(x, rev=reverse, by=key)
    end
    Py(nothing)
end

function pyjlvector_reverse(x::AbstractVector)
    reverse!(x)
    Py(nothing)
end

function pyjlvector_clear(x::AbstractVector)
    empty!(x)
    Py(nothing)
end

function pyjlvector_reversed(x::AbstractVector)
    Py(reverse(x))
end

function pyjlvector_insert(x::AbstractVector, k_::Py, v_::Py)
    k = pyconvertarg(Int, k_, "index")
    a = axes(x, 1)
    k′ = k < 0 ? (last(a) + 1 + k) : (first(a) + k)
    if checkbounds(Bool, x, k′) || k′ == last(a)+1
        v = pyconvertarg(eltype(x), v_, "value")
        insert!(x, k′, v)
        return Py(nothing)
    else
        errset(pybuiltins.IndexError, "array index out of bounds");
        return pynew()
    end
end

function pyjlvector_append(x::AbstractVector, v_::Py)
    v = pyconvertarg(Int, v_, "value")
    push!(x, v)
    Py(nothing)
end

function pyjlvector_extend(x::AbstractVector, vs_::Py)
    for v_ in vs_
        v = pyconvert_and_del(eltype(x), v_)
        push!(x, v)
    end
    Py(nothing)
end

function pyjlvector_pop(x::AbstractVector, k_::Py)
    k = pyconvertarg(Int, k_, "index")
    a = axes(x, 1)
    k′ = k < 0 ? (last(a) + 1 + k) : (first(a) + k)
    if checkbounds(Bool, x, k′)
        if k′ == last(a)
            v = pop!(x)
        elseif k′ == first(a)
            v = popfirst!(x)
        else
            v = x[k′]
            deleteat!(x, k′)
        end
        return Py(v)
    else
        errset(pybuiltins.IndexError, "pop from empty array")
        return pynew()
    end
end

function pyjlvector_remove(x::AbstractVector, v_::Py)
    v = @pyconvert eltype(x) v_ begin
        errset(pybuiltins.ValueError, "value not in array")
        return pynew()
    end
    v = pyconvert_result(r)
    k = findfirst(==(v), x)
    if k === nothing
        errset(pybuiltins.ValueError, "value not in array")
        return pynew()
    end
    deleteat!(x, k)
    Py(nothing)
end

function pyjlvector_index(x::AbstractVector, v_::Py)
    v = @pyconvert eltype(x) v_ begin
        errset(pybuiltins.ValueError, "value not in array")
        return pynew()
    end
    k = findfirst(==(v), x)
    if k === nothing
        errset(pybuiltins.ValueError, "value not in array")
        return pynew()
    end
    Py(k - first(axes(x, 1)))
end

function pyjlvector_count(x::AbstractVector, v_::Py)
    v = @pyconvert(eltype(x), v_, (return Py(0)))
    Py(count(==(v), x))
end

function init_jlwrap_vector()
    jl = pyjuliacallmodule
    filename = "$(@__FILE__):$(1+@__LINE__)"
    pybuiltins.exec(pybuiltins.compile("""
    class VectorValue(ArrayValue):
        __slots__ = ()
        __module__ = "juliacall"
        def resize(self, size):
            return self._jl_callmethod($(pyjl_methodnum(pyjlvector_resize)), size)
        def sort(self, reverse=False, key=None):
            return self._jl_callmethod($(pyjl_methodnum(pyjlvector_sort)), reverse, key)
        def reverse(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlvector_reverse)))
        def clear(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlvector_clear)))
        def __reversed__(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlvector_reversed)))
        def insert(self, index, value):
            return self._jl_callmethod($(pyjl_methodnum(pyjlvector_insert)), index, value)
        def append(self, value):
            return self._jl_callmethod($(pyjl_methodnum(pyjlvector_append)), value)
        def extend(self, values):
            return self._jl_callmethod($(pyjl_methodnum(pyjlvector_extend)), values)
        def pop(self, index=-1):
            return self._jl_callmethod($(pyjl_methodnum(pyjlvector_pop)), index)
        def remove(self, value):
            return self._jl_callmethod($(pyjl_methodnum(pyjlvector_remove)), value)
        def index(self, value):
            return self._jl_callmethod($(pyjl_methodnum(pyjlvector_index)), value)
        def count(self, value):
            return self._jl_callmethod($(pyjl_methodnum(pyjlvector_count)), value)
    import collections.abc
    collections.abc.MutableSequence.register(VectorValue)
    del collections
    """, filename, "exec"), jl.__dict__)
    pycopy!(pyjlvectortype, jl.VectorValue)
end

pyjl(v::AbstractVector) = pyjl(pyjlvectortype, v)


text:
module SparseSensors
using LinearAlgebra: Matrix, conj
greet() = print("Hello World!")
using LinearAlgebra
include("qr.jl")
include("ccqr.jl")
include("basis.jl")
include("reconstruct.jl")

export QRPivot,CostQRPivot,fit,VandermondeBasis

export get_sensors

function get_sensors(sampler)
    return sampler.pivots
end

function get_sensors(pivots::AbstractArray,n_sensors::Int)
    return pivots[1:n_sensors]
end

function get_sensors(sampler,n_sensors)
    return sampler.pivots[1:n_sensors]
end
end


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

include("../utilities.jl")

mock = MOIU.MockOptimizer(MOIU.Model{Float64}())
config = MOIT.TestConfig()

@testset "SOCR" begin
    bridged_mock = MOIB.Constraint.SOCR{Float64}(mock)

    MOIT.basic_constraint_tests(
        bridged_mock, config,
        include = [(F, S)
                   for F in [MOI.VectorOfVariables, MOI.VectorAffineFunction{Float64},
                             MOI.VectorQuadraticFunction{Float64}]
                   for S in [MOI.SecondOrderCone]])

    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 1/√2, 1/√2],
                          (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone) => [[1 - 1/√2, 1 + 1/√2, -1]],
                          (MOI.VectorAffineFunction{Float64}, MOI.Zeros)           => [[-√2]])
    MOIT.soc1vtest(bridged_mock, config)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 1/√2, 1/√2],
                          (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone) => [[1 - 1/√2, 1 + 1/√2, -1]],
                          (MOI.VectorAffineFunction{Float64}, MOI.Zeros)           => [[-√2]])
    MOIT.soc1ftest(bridged_mock, config)
    ci = first(MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64},
                                    MOI.SecondOrderCone}()))
    test_delete_bridge(bridged_mock, ci, 3,
                       ((MOI.VectorAffineFunction{Float64},
                         MOI.RotatedSecondOrderCone, 0),))
end


text:
# include("/media/mat/HDD/EnKF/src/EnKF.jl")
Pkg.activate("/media/mat/HDD/EnKF/")
using EnKF
using PotentialFlow

# import EnKF

"""
Create a specific structure for EnKF-based Vortex model,
Darakananda, et al. Data-assimilated low-order vortex modeling of separated flows, APS (2018)

"""

"""
    VortexENKFprop

    Define propagation operator VortexENKFprop for VortexENKF
"""

function (::PropagationFunction)(t::Float64, ENS::EnsembleState{N, TS}) where {N, TS}

    empty_inds = Vector{Int}[]
    len = Int[]
    ensnew = []
    # predicted_states (with merging)
    states = deepcopy(ENS.S)
    for (i, state) in enumerate(states)
        new_state = state_equation(state, t, config)
        new_blobs, lesp, tesp = state_to_blobs(new_state, config.δ)
        push!(empty_inds, findall(b -> circulation(b) == 0, new_blobs))
        push!(len, length(new_blobs))
        push!(ensnew, new_state)
    end

    # trim zeros
    @assert all(len .== len[1])
    toremove = intersect(empty_inds...)
    tokeep = filter(i -> i ∉ toremove, 1:len[1])
    for (i, state) in enumerate(ensnew)
        new_blobs, lesp, tesp = state_to_blobs(state, config.δ)
        ensnew[i] = blobs_to_state(new_blobs[tokeep], lesp, tesp)
    end

    # Wrap into an EnsembleState variable
    ENS.S .= deepcopy(ensnew)

    return ENS
end

VortexENKFprop = PropagationFunction()

"""
    VortexENKFinflation

    Define multiplico-additive inflation for VortexENKF
"""

function (RI::RecipeInflation)(ENS::EnsembleState{N, TS}) where {N, TS}
    # unwrap parameters
    αp, αΓ, αlesp, αtesp, β, n = RI.p

    n = ceil(Int, n)

    # get size of a member
    NS = size(ENS.S[1])[1]

    Ŝ = mean(deepcopy(ENS))

    # Covariance inflation
    for (_,s) in enumerate(ENS.S)
        s .= Ŝ .+ β*(s .- Ŝ)
        # don't perturb the state at the 4 first time steps
        if n > 5
        for i in 1:3:size(Ŝ, 1)-3
            s[i:i+1] += αp*randn(2)
            s[i+2]   += αΓ*randn()
        end
            s[end-1] = s[end-1] + αlesp*randn()
            s[end] = s[end] + αtesp*randn()
        end
    end

    return ENS
end

"""
    VortexENKFfiltering

    Define state filtering for VortexENKF
    Keep LESPc of each ensemble member between 0.0 and 3.0

"""

function (::FilteringFunction)(ENS::EnsembleState{N, TS}) where {N, TS}

    for s in ENS.S
        clamp(s[end-1],0.0,3.0)
    end
    return ENS
end

VortexENKFfilter = FilteringFunction()

"""
    Define specific constructor for EnKF

"""

function VortexENKF(N, NZ, ri, m, z, ϵ)

#     return ENKF{N, NZ}(VortexENKFprop, ri, VortexENKFfilter, m, z, ϵ, true, true, true)

    return ENKF{N, NZ}(VortexENKFprop, ri, VortexENKFfilter, m, z, ϵ, true, true, true, true);
end

max_stars_repo_name: mleprovost/EnKF

max_issues_repo_name: mleprovost/EnKF

max_forks_repo_name: mleprovost/EnKF