text:
using Kinetic, Plots, Flux
using KitBase.JLD2
using KitBase.ProgressMeter: @showprogress
using Base.Threads: @threads
cd(@__DIR__)

function sbr!(
    fwL::T1,
    ffL::T2,
    w::T3,
    prim::T3,
    f::T4,
    fwR::T1,
    ffR::T2,
    uVelo::T5,
    vVelo::T5,
    wVelo::T5,
    γ,
    μᵣ,
    ω,
    Kn_bz,
    nm,
    phi,
    psi,
    phipsi,
    dx,
    dt,
    RES,
    AVG,
    collision = :fsm::Symbol,
) where {T1,T2,T3,T4,T5}

    w_old = deepcopy(w)
    @. w += (fwL - fwR) / dx
    prim .= conserve_prim(w, γ)

    M = maxwellian(uVelo, vVelo, wVelo, prim)
    τ = vhs_collision_time(prim, μᵣ, ω)

    @. RES += (w - w_old)^2
    @. AVG += abs(w)

    Q = zero(f[:, :, :])
    boltzmann_fft!(Q, f, Kn_bz, nm, phi, psi, phipsi)

    for k in axes(f, 3), j in axes(f, 2), i in axes(f, 1)
        f[i, j, k] =
            (
                f[i, j, k] +
                (ffL[i, j, k] - ffR[i, j, k]) / dx +
                dt * (M[i, j, k] / τ * (1 - exp(-dt / τ)) + Q[i, j, k] * exp(-dt / τ))
            ) / (1.0 + dt / τ * (1 - exp(-dt / τ)))
    end

end

function up!(ks, ctr, face, dt, regime, p)
    kn_bzm, nm, phi, psi, phipsi = p

    res = zeros(5)
    avg = zeros(5)

    @inbounds Threads.@threads for i = 1:ks.ps.nx
        KitBase.step!(
            face[i].fw,
            face[i].ff,
            ctr[i].w,
            ctr[i].prim,
            ctr[i].f,
            face[i+1].fw,
            face[i+1].ff,
            ks.gas.γ,
            kn_bzm,
            nm,
            phi,
            psi,
            phipsi,
            ks.ps.dx[i],
            dt,
            res,
            avg,
            :fsm,
        )
        #=sbr!(
            face[i].fw,
            face[i].ff,
            ctr[i].w,
            ctr[i].prim,
            ctr[i].f,
            face[i+1].fw,
            face[i+1].ff,
            ks.vs.u,
            ks.vs.v,
            ks.vs.w,
            ks.gas.γ,
            ks.gas.μᵣ,
            ks.gas.ω,
            kn_bzm,
            nm,
            phi,
            psi,
            phipsi,
            ks.ps.dx[i],
            dt,
            res,
            avg,
            :fsm,
        )=#
    end
end

begin
    set = Setup(
        case = "sod",
        space = "1d1f3v",
        maxTime = 0.15,
        collision = "fsm",
        boundary = ["fix", "fix"],
    )
    ps = PSpace1D(0, 1, 200, 1)
    vs = VSpace3D(-6.0, 6.0, 64, -6.0, 6.0, 28, -6.0, 6.0, 28)
    gas = begin
        Kn = 1e-2
        Gas(Kn = Kn, K = 0.0, fsm = fsm_kernel(vs, ref_vhs_vis(Kn, 1.0, 0.5)))
    end
    ib = IB1F(ib_sod(set, ps, vs, gas)...)
    ks = SolverSet(set, ps, vs, gas, ib)
    ctr, face = init_fvm(ks)
end

t = 0.0
dt = timestep(ks, ctr, t)
nt = Int(ks.set.maxTime ÷ dt)
res = zero(ctr[1].w)
regime = ones(Int, axes(ks.ps.x))
regime0 = deepcopy(regime)

@time @showprogress for iter = 1:nt
    @inbounds @threads for i = 1:ks.ps.nx+1
        flux_kfvs!(
            face[i].fw,
            face[i].ff,
            ctr[i-1].f .+ 0.5 .* ctr[i-1].sf .* ks.ps.dx[i-1],
            ctr[i].f .- 0.5 .* ctr[i].sf .* ks.ps.dx[i],
            ks.vs.u,
            ks.vs.v,
            ks.vs.w,
            ks.vs.weights,
            dt,
            1.0,
            ctr[i-1].sf,
            ctr[i].sf,
        )
    end

    up!(ks, ctr, face, dt, regime, ks.gas.fsm)
end

"""
Kn = 1e-4
422.120548 seconds (41.79 M allocations: 1.842 TiB, 2.94% gc time, 0.11% compilation time)

Kn = 1e-3
420.792971 seconds (41.64 M allocations: 1.842 TiB, 2.93% gc time, 0.09% compilation time)

Kn = 1e-2
421.556153 seconds (41.96 M allocations: 1.842 TiB, 2.45% gc time, 0.08% compilation time)
1649.024909 seconds (164.98 M allocations: 7.351 TiB, 2.56% gc time, 0.02% compilation time)

NS:
0.552949 seconds (6.36 M allocations: 591.885 MiB, 17.68% gc time, 71.81% compilation time)
"""

function up!(ks, ctr, face, dt, regime, p)
    kn_bzm, nm, phi, psi, phipsi = p

    res = zeros(5)
    avg = zeros(5)

    @inbounds Threads.@threads for i = 1:ks.ps.nx
        KitBase.step!(
            face[i].fw,
            ctr[i].w,
            ctr[i].prim,
            face[i+1].fw,
            ks.gas.γ,
            ks.ps.dx[i],
            res,
            avg,
        )
    end
end

begin
    set = Setup(
        case = "sod",
        space = "1d1f3v",
        maxTime = 0.15,
        collision = "fsm",
        boundary = ["fix", "fix"],
    )
    ps = PSpace1D(0, 1, 200, 1)
    vs = VSpace3D(-6.0, 6.0, 64, -6.0, 6.0, 28, -6.0, 6.0, 28)
    gas = begin
        Kn = 1e-4
        Gas(Kn = Kn, K = 0.0, fsm = fsm_kernel(vs, ref_vhs_vis(Kn, 1.0, 0.5)))
    end
    ib = IB1F(ib_sod(set, ps, vs, gas)...)
    ks = SolverSet(set, ps, vs, gas, ib)
    ctr, face = init_fvm(ks)
end

t = 0.0
dt = timestep(ks, ctr, t)
nt = Int(ks.set.maxTime ÷ dt)
res = zero(ctr[1].w)
regime = ones(Int, axes(ks.ps.x))
regime0 = deepcopy(regime)

@time @showprogress for iter = 1:nt
    @inbounds @threads for i = 1:ks.ps.nx+1
        flux_gks!(
            face[i].fw,
            ctr[i-1].w .+ 0.5 .* ctr[i-1].sw .* ks.ps.dx[i-1],
            ctr[i].w .- 0.5 .* ctr[i].sw .* ks.ps.dx[i],
            ks.gas.K,
            ks.gas.γ,
            ks.gas.μᵣ,
            ks.gas.ω,
            dt,
            ks.ps.dx[i-1] / 2,
            ks.ps.dx[i] / 2,
            1.0,
            ctr[i-1].sw,
            ctr[i].sw,
        )
    end

    up!(ks, ctr, face, dt, regime, ks.gas.fsm)
end

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
using TimeSeriesResampler
using Test

include("resample.jl")

ext: jl


text:
using Printf

nobjs = 4
a = collect(1:nobjs)
println("All permutations of ", nobjs, " objects:")
for i in 1:factorial(nobjs)
    p = nthperm(a, i)
    prank = nthperm(p)
    print(@sprintf("%5d => ", i))
    println(p, " (", prank, ")")
end

nobjs = 12
nsamp = 4
ptaken = Int[]
println()
println(nsamp, " random permutations of ", nobjs, " objects:")
for i in 1:nsamp
    p = randperm(nobjs)
    prank = nthperm(p)
    while prank in ptaken
        p = randperm(nobjs)
        prank = nthperm(p)
    end
    push!(ptaken, prank)
    println("         ", p, " (", prank, ")")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#
# product
#

"""
    sitedict(::Vector{<:Index})

Return a dictionary that maps a Vector of indices to
the integer position of the Index in the Vector.
"""
function sitedict(sites::Vector{IndexT}) where {IndexT <: Index}
  d = Dict{IndexT, Int}()
  for (n, s) in enumerate(sites)
    d[s] = n
  end
  return d
end

"""
    product(o::ITensor, ψ::MPS)

Apply the ITensor `o` to the MPS `ψ`.
"""
function ITensors.product(o::ITensor{N},
                          ψ::Union{MPS, MPO};
                          kwargs...) where {N}
  ns = findcommonsiteinds(o, ψ)
  return product(o, ψ, ns; kwargs...)
end

"""
    product(o::ITensor, ψ::MPS, ns::Tuple;
            move_sites_back::Bool = true)

Get the product of the operator `o` with the MPS `ψ`,
where the operator is applied to the sites `ns`.

If `ns` are non-contiguous, the sites of the MPS are
moved to be contiguous. By default, the sites are moved
back to their original locations. You can leave them where
they are by setting the keyword argument `move_sites_back`
to false.
"""
function ITensors.product(o::ITensor,
                          ψ::Union{MPS, MPO},
                          ns;
                          move_sites_back::Bool = true,
                          apply_dag::Bool = false,
                          kwargs...)
  N = length(ns)
  ns = sort(ns)

  # TODO: make this smarter by minimizing
  # distance to orthogonalization.
  # For example, if ITensors.orthocenter(ψ) > ns[end],
  # set to ns[end].
  # TODO: make a movesites(ψ, ns, ns′) function.
  ψ = orthogonalize(ψ, ns[1])
  diff_ns = diff(ns)
  ns′ = ns
  if any(!=(1), diff_ns)
    ns′ = [ns[1] + n - 1 for n in 1:N]
    ψ = movesites(ψ, ns, ns′; kwargs...)
  end
  ϕ = ψ[ns′[1]]
  for n in 2:N
    ϕ *= ψ[ns′[n]]
  end
  ϕ = product(o, ϕ; apply_dag = apply_dag)
  ψ = replacesites(ψ, ϕ;
                   firstsite = ns′[1],
                   lastsite = ns′[end],
                   kwargs...)

  if move_sites_back
    error("move_sites_back is not supported yet")
  end

  return ψ
end

function orthogonalize(ψ0::Union{MPS, MPO}, args...; kwargs...)
  ψ = copy(ψ0)
  orthogonalize!(ψ, args...; kwargs...)
  return ψ
end

"""
    product(ops::Vector{<:ITensor}, ψ::MPS;
            apply_dag::Bool = false)

Apply the ITensors `ops` to the MPS `ψ`.
"""
function ITensors.product(ops::Vector{<:ITensor},
                          ψ::Union{MPS, MPO};
                          apply_dag::Bool = false,
                          kwargs...)
  ψ0 = ψ
  s0 = siteinds(ψ0)
  for o in ops
    ψ = product(o, ψ; move_sites_back = false,
                      apply_dag = apply_dag,
                      kwargs...)
  end
  s = siteinds(ψ)
  ns = 1:length(ψ)
  ns′ = findsiteinds(ψ0, s)
  # Move the sites back to their original positions
  ψ = movesites(ψ, ns, ns′; kwargs...)
  return ψ
end

"""
    product(A::ITensor, B::ITensor)

Get the product of ITensor `A` and ITensor `B`.

`A` should have pairs of unprimed and primed indices,
and `B` should either have pairs of primed or unprimed
indices that are shared with `A` or just unprimed indices
shared with `A`.

If `B` has pairs of unprimed and primed indices shared with `A`,
it is treated as a matrix, and `A` is contracted as a matrix-matrix
product, with a result that has the same unprimed and primed
indices as `A` and `B`.

If `B` just has unprimed indices shared with `A`, then
`A` is contracted as a matrix-vector product, and the
result is unprimed so that it has the unprimed indices
that were shared between `A` and `B`.

`A` and `B` can both have unshared indices, which are not
contracted or effected by the operation. This would correspond
to batched matrix-matrix or matrix-vector products.
"""
function ITensors.product(A::ITensor,
                          B::ITensor;
                          apply_dag::Bool = false)
  commonABis = commoninds(A, B)
  if all(hasplev(0), commonABis)
    # B is a state
    apply_dag && error("Cannot apply to both sides of a state")
    return noprime(A * B; inds = commonABis')
  end
  # B is an operator
  commonABis = filterinds(commonABis; plev = 0)
  uniqueAis = uniqueinds(A, B)
  A′ = prime(A; inds = not(uniqueAis))
  AB = setprime(A′ * B, 1; inds = commonABis'')
  if apply_dag
    uniqueBis = uniqueinds(B, A)
    AB = prime(AB; inds = not(unioninds(uniqueAis, uniqueBis)))
    Adag = swapprime(dag(A), 0, 1; inds = not(uniqueAis))
    AB = setprime(AB * Adag, 1; inds = commonABis'')
  end
  return AB
end

"""
    product(ops::Vector{<:ITensor}, ψ::ITensor)

Apply the ITensors `ops` to the ITensor `ψ`.
"""
function ITensors.product(ops::Vector{<:ITensor},
                          ψ::ITensor;
                          apply_dag::Bool = false)
  for o in ops
    ψ = product(o, ψ; apply_dag = apply_dag)
  end
  return ψ
end

const apply = product


text:
struct DummySystem <: AbstractSystem; end

"""
Generic struct for a system
"""
mutable struct System <: AbstractSystem
  dim::Int64
  t::Float64
  x::Vector{Float64}
  f::Vector{Float64}
  function System(dim::Int64,t::Float64,x::Vector{Float64},f::Vector{Float64})
    new(dim,t,x,f)
  end
end

function System(model::AbstractPotential)
  dim::Int64 = getdimensionality(model)
  x::Vector{Float64} = zeros(dim)
  f::Vector{Float64} = zeros(dim)
  System(dim,0.0,x,f)
end

"""
Generic struct for a system with a thermostat
"""
mutable struct ThermostattedSystem <: AbstractThermostattedSystem
  dim::Int64
  t::Float64
  x::Vector{Float64}
  f::Vector{Float64}
  thermostat::AbstractThermostat
  function ThermostattedSystem(dim::Int64,t::Float64,x::Vector{Float64},f::Vector{Float64},thermostat::AbstractThermostat)
    new(dim,t,x,f,thermostat)
  end
end

function ThermostattedSystem(model::AbstractPotential,thermostat::AbstractThermostat)
  dim::Int64 = getdimensionality(model)
  x::Vector{Float64} = zeros(dim)
  f::Vector{Float64} = zeros(dim)
  ThermostattedSystem(dim,0.0,x,f,thermostat)
end

#"""
#Struct for an overdamped system with a langevin thermostat
#"""
#mutable struct OverdampedLangevinSystem <: AbstractThermostattedSystem
#  dim::Int64
#  t::Float64
#  x::Vector{Float64}
#  f::Vector{Float64}
#  thermostat::AbstractLangevin
#  function OverdampedLangevinSystem(dim::Int64,t::Float64,x::Vector{Float64},f::Vector{Float64},thermostat::AbstractLangevin)
#    new(dim,t,x,f,thermostat)
#  end
#end
#
#function OverdampedLangevinSystem(model::AbstractPotential,thermostat::AbstractLangevin)
#  dim::Int64 = getdimensionality(model)
#  x::Vector{Float64} = zeros(dim)
#  f::Vector{Float64} = zeros(dim)
#  ThermostattedSystem(dim,0.0,x,f,thermostat)
#end

"""
Functions for interacting with the force
"""
function zeroforce!(system::AbstractSystem)
  system.f .= 0.0
end

function addforce!(dt::Float64,system::AbstractSystem,model::AbstractSinglePotential)
  system.x .+= dt .* model.f
end

function addforce!(dt::Float64,system::AbstractSystem,model::AbstractTrainablePotential)
  if model.condition(system)
    system.x .+= dt .* model.f
  end
end

# TODO maybe take this function out since it's a repeat of the one above
function addforce!(dt::Float64,system::ThermostattedSystem,model::AbstractSinglePotential)
  system.x .+= dt .* model.f ./ system.thermostat.gamma
#  system.x .+= dt .* model.f
end

function addforce!(dt::Float64,system::ThermostattedSystem,model::AbstractTrainablePotential)
#  println(system.t)
  if model.condition(system)
#    println("hey")
    system.x .+= dt .* model.f ./ system.thermostat.gamma
  end
end

# TODO do unit tests on Overdamped system before allowing this
#function addforce!(dt::Float64,system::OverdampedLangevinSystem,model::AbstractPotential)
#  system.x .+= dt .* model.f ./ system.thermostat.gamma
#end


text:
using ModelingToolkit
using Graphs
using DiffEqBase
using Test
using UnPack

# Define some variables
@parameters t L g
@variables x(t) y(t) w(t) z(t) T(t) xˍt(t) yˍt(t) xˍˍt(t) yˍˍt(t)
D = Differential(t)

eqs2 = [D(D(x)) ~ T * x,
    D(D(y)) ~ T * y - g,
    0 ~ x^2 + y^2 - L^2]
pendulum2 = ODESystem(eqs2, t, [x, y, T], [L, g], name = :pendulum)
lowered_sys = ModelingToolkit.ode_order_lowering(pendulum2)

lowered_eqs = [D(xˍt) ~ T * x,
    D(yˍt) ~ T * y - g,
    D(x) ~ xˍt,
    D(y) ~ yˍt,
    0 ~ x^2 + y^2 - L^2]
@test ODESystem(lowered_eqs, t, [xˍt, yˍt, x, y, T], [L, g], name = :pendulum) ==
      lowered_sys
@test isequal(equations(lowered_sys), lowered_eqs)

# Simple pendulum in cartesian coordinates
eqs = [D(x) ~ w,
    D(y) ~ z,
    D(w) ~ T * x,
    D(z) ~ T * y - g,
    0 ~ x^2 + y^2 - L^2]
pendulum = ODESystem(eqs, t, [x, y, w, z, T], [L, g], name = :pendulum)

state = TearingState(pendulum)
@unpack graph, var_to_diff = state.structure
@test StructuralTransformations.maximal_matching(graph, eq -> true,
                                                 v -> var_to_diff[v] === nothing) ==
      map(x -> x == 0 ? StructuralTransformations.unassigned : x,
          [1, 2, 3, 4, 0, 0, 0, 0, 0])

using ModelingToolkit
@parameters t L g
@variables x(t) y(t) w(t) z(t) T(t) xˍt(t) yˍt(t)
D = Differential(t)
idx1_pendulum = [D(x) ~ w,
    D(y) ~ z,
    #0 ~ x^2 + y^2 - L^2,
    D(w) ~ T * x,
    D(z) ~ T * y - g,
    # intermediate 1: 0 ~ 2x*D(x) + 2y*D(y) - 0,
    # intermediate 2(a): 0 ~ 2x*w + 2y*z - 0, (substitute D(x) and D(y))
    #0 ~ 2x*w + 2y*z,
    # D(D(x)) ~ D(w) and substitute the rhs
    D(xˍt) ~ T * x,
    # D(D(y)) ~ D(z) and substitute the rhs
    D(yˍt) ~ T * y - g,
    # 2x*D(D(x)) + 2*D(x)*D(x) + 2y*D(D(y)) + 2*D(y)*D(y) and
    # substitute the rhs
    0 ~ 2x * (T * x) + 2 * xˍt * xˍt + 2y * (T * y - g) + 2 * yˍt * yˍt]
@named idx1_pendulum = ODESystem(idx1_pendulum, t, [x, y, w, z, xˍt, yˍt, T], [L, g])
first_order_idx1_pendulum = ode_order_lowering(idx1_pendulum)

using OrdinaryDiffEq
using LinearAlgebra
prob = ODEProblem(ODEFunction(first_order_idx1_pendulum),
                  #  [x, y, w, z, xˍt, yˍt, T]
                  [1, 0, 0, 0, 0, 0, 0.0],# 0, 0, 0, 0],
                  (0, 10.0),
                  [1, 9.8],
                  mass_matrix = calculate_massmatrix(first_order_idx1_pendulum))
sol = solve(prob, Rodas5());
#plot(sol, vars=(1, 2))

new_sys = dae_index_lowering(ModelingToolkit.ode_order_lowering(pendulum2))

prob_auto = ODEProblem(new_sys,
                       [D(x) => 0,
                           D(y) => 0,
                           x => 1,
                           y => 0,
                           T => 0.0],
                       (0, 100.0),
                       [1, 9.8])
sol = solve(prob_auto, Rodas5());
#plot(sol, vars=(x, y))

# Define some variables
@parameters t L g
@variables x(t) y(t) T(t)
D = Differential(t)

eqs2 = [D(D(x)) ~ T * x,
    D(D(y)) ~ T * y - g,
    0 ~ x^2 + y^2 - L^2]
pendulum2 = ODESystem(eqs2, t, [x, y, T], [L, g], name = :pendulum)

# Turn into a first order differential equation system
first_order_sys = ModelingToolkit.ode_order_lowering(pendulum2)

# Perform index reduction to get an Index 1 DAE
new_sys = dae_index_lowering(first_order_sys)

u0 = [
    D(x) => 0.0,
    D(y) => 0.0,
    x => 1.0,
    y => 0.0,
    T => 0.0,
]

p = [
    L => 1.0,
    g => 9.8,
]

prob_auto = ODEProblem(new_sys, u0, (0.0, 10.0), p)
sol = solve(prob_auto, Rodas5());
#plot(sol, vars=(D(x), y))

let pss_pendulum2 = partial_state_selection(pendulum2)
    # This currently selects `T` rather than `x` at top level. Needs tearing priorities to fix.
    @test_broken length(equations(pss_pendulum2)) == 3
    @test_broken length(equations(ModelingToolkit.ode_order_lowering(pss_pendulum2))) == 4
end

eqs = [D(x) ~ w,
    D(y) ~ z,
    D(w) ~ T * x,
    D(z) ~ T * y - g,
    0 ~ x^2 + y^2 - L^2]
pendulum = ODESystem(eqs, t, [x, y, w, z, T], [L, g], name = :pendulum)

let pss_pendulum = partial_state_selection(pendulum)
    # This currently selects `T` rather than `x` at top level. Needs tearing priorities to fix.
    @test_broken length(equations(pss_pendulum)) == 3
end

sys = structural_simplify(pendulum2)
@test length(equations(sys)) == 5
@test length(states(sys)) == 5

u0 = [
    D(x) => 0.0,
    D(y) => 0.0,
    x => sqrt(2) / 2,
    y => sqrt(2) / 2,
    T => 0.0,
]
p = [
    L => 1.0,
    g => 9.8,
]

prob_auto = DAEProblem(sys, zeros(length(u0)), u0, (0.0, 0.2), p)
sol = solve(prob_auto, DFBDF())
@test norm(sol[x] .^ 2 + sol[y] .^ 2 .- 1) < 1e-2


text:
push!( LOAD_PATH, "../" )
import InformationInequalities  # Import the module to  have access to functions

"""
List of Elemental Information measures (EIM) for a given `n` number of random variables. EIM comprise of conditional entropies `H(X1...,Xn|Y1...Ym)` and conditional mutual information `I(X₁...Xₙ;Y₁....Yₘ|Z₁...Zₖ)` measures.

```julia-repl
julia> ElementalMeasures(2)
["H(X1)","H(X1|X2)","H(X2)","H(X2|X1)","I(X1;X2)","I(X2;X1)"]

julia> ElementalMeasures(3,"🐘")
["H(🐘1)"
"H(🐘1|🐘2)"
"H(🐘1|🐘2,🐘3)"
"H(🐘1|🐘3)"
"H(🐘2)"
"H(🐘2|🐘1)"
"H(🐘2|🐘1,🐘3)"
"H(🐘2|🐘3)"
"H(🐘3)"
"H(🐘3|🐘1)"
"H(🐘3|🐘1,🐘2)"
"H(🐘3|🐘2)"
"I(🐘1;🐘2)"
"I(🐘1;🐘2|🐘3)"
"I(🐘1;🐘3)"
"I(🐘1;🐘3|🐘2)"
"I(🐘2;🐘1)"
"I(🐘2;🐘1|🐘3)"
"I(🐘2;🐘3)"
"I(🐘2;🐘3|🐘1)"
"I(🐘3;🐘1)"
"I(🐘3;🐘1|🐘2)"
"I(🐘3;🐘2)"
"I(🐘3;🐘2|🐘1)"]
```
"""
function ElementalMeasures(n::Int64=2,RV::AbstractString="X")
	yI=ConditionalMutualInformationList(n,RV)
	yH=ConditionalEntropyList(n,RV)
	y=[yH;yI][:]
	return y
end

"""
For a given number `n` of random variables, list down the maximum number of elemental information measures in minimal canonical form.
```julia-repl
julia> u,v,m,n=numEIM(2)
```
"""
function numEIM(n::Int64=2)
	~,~,k,l=minimal_EIM_list_canonical(n)
	return k,l
end

"""
    MutualInformation_ElemToCanon(s::String="I(Xi;X🐬|π)")
Convert a mutual information expression s to elemental form.
If `s` is unspecified, it performs a default expression.
##### Examples
```julia-repl
julia> Elemental2Canonical("I(Xi;X🐬|π)")
"H(Xi,π)+H(X🐬,π)-H(Xi,X🐬,π)-H(π)"
```
"""
function Elemental2Canonical_MI(s::String="I(Xi;X🐬|π)")
	s0=replace(s," " =>"")
	if (occursin(r"\|", s0))
		m1=match(r"(I.+;)",s0)
		m2=match(r"(;.+\|)",s0)
		m3=match(r"(\|.+\))",s0)
		a1=replace(m1.captures[1],"I("=>"",";"=>"")
		a2=replace(m2.captures[1],";"=>"","|"=>"")
		a3=replace(m3.captures[1],")"=>"","|"=>"")

		P1="H(" * a1* "," * a3 * ")"
		P2="H(" * a2* "," * a3 * ")"
		P3="H(" * a1 * "," * a2 * "," * a3 * ")"
		P4="H(" * a3 * ")"
		C= P1 * "+" * P2 * "-" * P3 * "-" * P4
	else
		m1=match(r"(I.+;)",s0)
		m2=match(r"(;.+\))",s0)
		a1=replace(m1.captures[1],"I("=>"",";"=>"")
		a2=replace(m2.captures[1],";"=>"",")"=>"")

		P1="H(" * a1* ")"
		P2="H(" * a2*  ")"
		P3="H(" * a1 * "," * a2 *")"
		C= P1 * "+" * P2  * "-" * P3
	end
	return C
end

"""
Elemental2Canonical_H(s::String="I(Xi;X🐬|π)")
Convert a Entropy or Conditional entropy expression s to elemental form.
If `s` is unspecified, it performs a default expression.
##### Examples
```julia-repl
julia> Elemental2Canonical_H("H(Xi,X🐬|π,γ)")
"H(Xi,X🐬,π,γ)-H(π,γ)"
```
"""
function Elemental2Canonical_H(s::String="H(ρ12|Xπ)")
s0=replace(s," " =>"")
if (occursin(r"\|", s0))
    m = match(r"([H-I]\(.*\|([^\)]*)\))", s0)
    y=replace(m.captures[1], '|' => ',') * "-H($(m.captures[2]))"
else
    y=s
end
return y
end

"""
    Elemental2Canonical(s::String="I(Xi;X🐬|π)")
Convert elemental Information measure to canonical form.
If `s` is unspecified, it performs a default expression.
##### Examples
```julia-repl
julia> Elemental2Canonical("I(Xi;X🐬|π)")
"H(Xi,π)+H(X🐬,π)-H(Xi,X🐬,π)-H(π)"

julia> Elemental2Canonical("H(Xi,X2|π,β)")
"H(Xi,X2,π,β) - H(π,β)"

julia>Elemental2Canonical("I(Xi;X2,ρ,🍎|Zπ,β,🍩)")
"H(Xi,Zπ,β,🍩)+H(X2,ρ,🍎,Zπ,β,🍩)-H(Xi,X2,ρ,🍎,Zπ,β,🍩)-H(Zπ,β,🍩)"
```
"""
function Elemental2Canonical(s::String="I(Xi;X🐬|π)")
	s0=replace(s," " =>"")
	if(startswith(s,"I"))
		C=Elemental2Canonical_MI(s)
	elseif(startswith(s,"H"))
		C=Elemental2Canonical_H(s)
	else
		error("$(s) is an invalid Information measure. Information expression must be of the form H(...) or I(...); e.g. H(Xi,X🐬|π)")
	end
	return C
end

#=
"""
	InformationExpressionToCanonical(F::AbstractString="-2.32I(X;Y|Z1,Z2)")
Convert a single scaled Information measure to its canonnical decomposition form. elemental Information measure to canonical form.
If `F` is unspecified, it performs on a default expression.

##### Examples
```julia-repl
julia> InformationExpressionToCanonical("2.32I(X;Y🍎|Z1,Z2)")
"2.32I(X;Y🍎|Z1,Z2)"

julia> InformationExpressionToCanonical("3H(β😄,Y🍎|Z1,Z2)")
"3H(β😄,Y🍎|Z1,Z2)"
```
"""
=#
function InformationExpressionToCanonical(F::AbstractString="-2.32I(X;Y|Z1,Z2)")
	# Do check and pick only one Iexpr H() or I() with scaling

	 kk=match(r"(^[+|-]?[H|I][(])",F)  #Check if starts with bare I() or H(
	if kk != nothing # starts with bare I() or H(
		sExpr=replace(string(F),"-I("=>"I(","-H("=>"H(")
		# sCoef="1"
		# kk.match[1]
		ρ = (kk.match[1]=='-') ? -1 : 1

		sCoef = string(ρ)
		 Elemental2Canonical(sExpr)
		# @show sCoef=(kk.match[1]=='-') ? "-1" : "1"
	else

	strIH=r"([H|I]+\([\w](.)*[\)])"
	strN=r"([-|+]?[\d]*[\.]?[\d ?]*)"
	Mexp=match(strIH,F)
	Mcoef=match(strN,F)

	sExpr=string(Mexp.captures[1])
	sCoef=string(Mcoef.captures[1])
	# scalar=parse(Float64,sCoef)
	end
	CanExpr=Elemental2Canonical(sExpr)
	Z0=replace(CanExpr,"H(" =>sCoef *"H(")
	Z1=replace(Z0,"--"=>"+")
	Z=replace(Z1,"+-"=>"-")

	return Z
end

function InformationMeasureToCanonical(S::AbstractString="-2.32I(X;Y|Z1,Z2)")
	# Do check and pick only one Iexpr H() or I() with scaling
    F=replace(S," "=>"")
	# Return an error if there are more than one IMeasures
	F1=replace(F,"-"=>"+-")
	F2= split(F1,"+")
	F3 = filter(!isempty, F2)
	if(length(F3)>1)
		error("A single Information Meausure expected; I(X...;Y...|Z...) or H(X,Y...|U,V...)")
	end

	 kk=match(r"(^[+|-]?[H|I][(])",F)  #Check if starts with bare ±I() or ±H(
	if kk != nothing # starts with bare I() or H(
		sExpr=replace(string(F),"-I("=>"I(","-H("=>"H(","+I("=>"I(","+H("=>"H(")
		# sCoef="1"
		# kk.match[1]
		ρ = (kk.match[1]=='-') ? -1 : 1
		sCoef = string(ρ)
		 # Elemental2Canonical(sExpr)
		# @show sCoef=(kk.match[1]=='-') ? "-1" : "1"
	else

		strIH=r"([H|I]+\([\w](.)*[\)])"
		strN=r"([-|+]?[\d]*[\. | \/]?[\d ?]*)"
		Mexp=match(strIH,F)
		Mcoef=match(strN,F)

		sExpr=string(Mexp.captures[1])
		sCoef=string(Mcoef.captures[1])
		# @show parse(Float64,sCoef)
		# @show scalar=parse(Float64,sCoef)
	end
	CanExpr=Elemental2Canonical(sExpr)
	Z0=replace(CanExpr,"H(" =>sCoef *"H(")
	Z1=replace(Z0,"--"=>"+")
	Z2=replace(Z1,"+-"=>"-")
	Z3=replace(Z2,"-+"=>"-")
	Z4= replace(Z3," "=>"")
	Z=replace(Z4,"++"=>"+")
	return Z
end

"""
LinearInformationExpressionToCanonical(A)
```julia-repl
julia>LinearInformationExpressionToCanonical("I(X;Y|Z)-2.3H(U,V)-2H(u)")
"1H(X,Z)+1H(Y,Z)-1H(X,Y,Z)-1H(Z)-2.3H(U,V)-2H(u)"
```
"""
function LinearInformationExpressionToCanonical(A)
    A1=replace(A,"-"=>"+-","+I("=>"+1I(","-I("=>"+-1I(","+H("=>"+1H(","-H("=>"+-1H(")
    A1a =replace(A1," "=>"")
    A2=split.(A1a,"+")
	A2a = filter(!isempty, A2)
    # A3=InformationExpressionToCanonical.(A2)
	A3=InformationMeasureToCanonical.(A2a)
	# A4=prod(A3)
	A4=""
	for ii=1:length(A3)
		A4= A4 * "+" * A3[ii]
	end
	A5=replace(A4,"+-"=>"-")
	return A5
end

"""
For a given number `n` of random variables, list down all the elemental information measures and their corresponding entropic decompositions. The entropic vectors are identified with the prefix `h` and follows lexicographic mapping. e.g., `H(X1,X3,X7)=h137`. Note that, for now this lexicographic mapping works only for `n < 10`.
"""
function entropy_vector(n::Int64,RV::AbstractString="X")
	X=Elemental2Canonical.(ElementalMeasures(n,RV))
	Y=sort(replace.(X,","=>"","X"=>"",")"=>"","("=>"","H"=>"h"))
	return Y
end

"""
For a given number `n<10` of random variables, list down all the unique elemental information measures and their corresponding entropic decompositions. The entropic vectors are identified with the prefix `h` and follows lexicographic mapping. e.g., `H(X1,X3,X7)=h137`. Note that, for now this lexicographic mapping works only for `n < 10`.
"""
function unique_entropy_vector(n::Int64=3)
    EMeasures=Elemental2Canonical.(ElementalMeasures(n))
    EMnumber=(replace.(EMeasures,","=>"","X"=>"",")"=>"","("=>"","H"=>"","+"=>",","-"=>","))
    y=[];
    for kk in 1:length(EMnumber)
        y=[y;"h" .* unique(sort(join.(sort.(digits.(parse.(Int64,split(EMnumber[kk],",")))))))]
    end
    Hfull=y;
    Hminimal = unique(Hfull)
    return EMeasures,Hfull,Hminimal
end

function find_entropic_vectorInt(n::Int64=3,entropicVector="1h12-1h2")
	# vv=toy_entropic_matrix(n)[5]
	x= entropicVector * "+"
		listE=unique_entropy_vector(n)[3]
	   sortedListE="h" .* string.(sort(parse.(Int64,replace.(listE,"h"=>"")),rev=true))
	G=zeros(1,length(listE))

	for idx=1:length(listE)
		mx=match(Regex("([-|+]?\\d$(sortedListE[idx])[-|+]+)"),x) # special quantifier Z to identify the end word

	if(mx == nothing)
		G[idx]=0
	else
		mz=match(r"([-|+]?\dh)",mx.captures[1])
		G[idx]= parse(Int64,replace(mz.captures[1],"h"=>""))
	end
		end
	return G[:]
end

"""
Find the Entropic matrix `G`` for a given `n`
```julia-repl
julia> entropic_matrix(3)
```
"""
function entropic_matrix(m::Int64=2)
	Q=minimal_EIM_list_canonical(m)[2]
	Qtmp=replace.(Q,r"^\+h" => s"1h") # Substitute the first +h
	Qtmp1=replace.(Qtmp,r"\+h" => s"+1h") # Substitute all +h
	Qout=replace.(Qtmp1,r"\-h" => s"-1h")
	return Qout
end

"""
For a given number `n` of random variables, list down all the elemental information measures in minimal canonical form.
```julia-repl
julia> u,v,m,n=minimal_EIM_list_canonical(2)
```
"""
function minimal_EIM_list_canonical(n::Int64=2)
	x=entropy_vector(n);
	x2=order_entropic_expression.(order_entropic.(x))
	y=unique(x2)
	return x,y,length(x),length(y)
end

"""
Each entropy word in an entropic vector is sorted.

#### Examples
```julia-repl
julia> order_entropic("h24-h32-h132-h2")
"h24 - h23 - h123 - h2"

```
"""
function order_entropic(s::AbstractString="Hello world")

	s0=replace(s,"h"=>"h ","+"=>" + ","-"=>" - ")
	x=split(s0," ") # Split into words delim=" "

	Y=""
	for wIdx = 1:length(x)
		# w=x[wIdx] # word
		y=string(sort([a for a in x[wIdx]])...)  # sort each word
		Y= Y * " " * y
	end
	Yout=replace(Y,"h "=>"h"," " =>"")

	return Yout
end

"""
Express an entropic expression in the lexicographic order of entropic vectors
```julia-repl
julia> order_entropic_expression("h12-h32-123")
```
"""
function order_entropic_expression(X::AbstractString="h12-h32-123")
	Y=join.(sort.(split.(replace.(X,"-h"=>" -h","+h"=>" +h","h"=>" +h")," ")))
	return Y
end

"""
Find the Entropic matrix `G` for a given `n`
```julia-repl
julia> find_matrixG(3)
```
"""
function find_matrixG(n::Int64=3)
	L=length(entropic_matrix(n))
	K=length(unique_entropy_vector(n)[3]) # Number of unique entropic dimensions
	Gv=[]
	# @show L
	for idx=1:L
		U = entropic_matrix(n)[idx]
	    gg = find_entropic_vector(n,U)
		# @show n,L,idx, gg
		Gv=[Gv;gg]
		# Gm[idx,:]=gg
	end
	G=reshape(Gv,K,L)
	H0=Int.(parse.(Float64,string.(G)))'
    H=Int.(H0)
	return H,K,L
end

function find_entropic_vector_int(n::Int64=3,entropicVector="1h12-1h2")
	x= entropicVector * "+"
		listE=unique_entropy_vector(n)[3]
	   sortedListE="h" .* string.(sort(parse.(Int64,replace.(listE,"h"=>"")),rev=true))
	G=zeros(1,length(listE))

	for idx=1:length(listE)
		mx=match(Regex("([-|+]?\\d$(sortedListE[idx])[-|+]+)"),x) # special quantifier Z to identify the end word

	if(mx == nothing)
		G[idx]=0
	else
		mz=match(r"([-|+]?\dh)",mx.captures[1])
		G[idx]= parse(Int64,replace(mz.captures[1],"h"=>""))
	end
		end
	return G[:]
end

"""
Given a linear expression in canonical form, it finds the co-ordinates in the entropic space `Γ`
"""
function find_entropic_vector(n::Int64=3,entropicVector="1h12-1h2")
	x= entropicVector * "+"
	listE=unique_entropy_vector(n)[3]
	sortedListE="h" .* string.(sort(parse.(Int64,replace.(listE,"h"=>"")),rev=true))
    G=zeros(1,length(listE))
	for idxx=1:length(listE)
	    mx=match(Regex("([-|+]?\\d*[.]?\\d*$(sortedListE[idxx])[-|+]+)"),x)

	    if(mx == nothing)
		    G[idxx]=0.0
	    else
		# mz=match(r"([-|+]?\dh)",mx.captures[1])
		    mz=match(r"([-|+]?[\d]*[.]?[\d]*h)",mx.captures[1])
		    G[idxx]= parse(Float64,replace(mz.captures[1],"h"=>""))
	    end
	end
    return G[:]
end

"""
Γ₂ geometry. This is the simplest case with two random variables (say `X,Y`) forming a geometry in three dimension. The geometric space is spanned by entropy vectors `H(X)`, `H(Y)` and `H(X,Y)`. Γ₂ is a 3D cone in the positive orthant.
This function is used for visualizing the entropic space in 3D.
"""
function GeometryConeGamma2(x,y,z)
	H=[ 0.0   0   1
 		-1   1   1
  		1  -1   0
  		0   1   0
  		1   0  -1]
	X=(0.0,0.0,0)
	Y=[]
	Z=[]
	p = 0.0
	idx=1
	if(-z+y+x ≥0 && z -y ≥ 0 && z -x ≥ 0 && z ≥0 && y≥0 && x ≥ 0)
		p=1.0

		idx,X=[X;x]
		 Y=[Y;y]
		# Z=[Z;z]
		# (x,y,z)
		# Z=push!(Z,z)
	end
	if(-z+y+x ≥0)
		if ( z -y ≥ 0)
			if(z -x ≥ 0)
				if(z ≥0)
					if (y≥0)
						if (x ≥ 0)
							Y=push!(Z,tuple(x,y,z))
						end
					end
				end
			end
		end
	end

	return Y
end

"""
Set of discrete points in `Γ₂`	 confined within a hypercube
"""
function elementsGamma2(;rmin::Float64=0.0,rmax::Float64=1.0,rstep::Float64 = 0.1)
	x = rmin:rstep:rmax
    y = rmin:rstep:rmax
    z = rmin:rstep:rmax
	vol = [GeometryConeGamma2(ix,iy,iz) for ix in x, iy in y, iz in z]
	U=vol[findall( x -> x != [], vol )]
	V=hcat(U...)
	W=hcat(collect.(V)...)
	return W
end

"""
    E,UE,V,λ =entropic_terms(S::AbstractString)

Find the additive entropic terms in the canonical expression. This function is internally used in `simplifyH(S::AbstractString)`.
#### Arguments
- S  = Any linear Information Expression. These are linear combination of ``I(X_1,...X_k;Y_1,..Y_l|Z_1,...Z_n)`` and ``H(X_1,X_2,...X_m|Z_1,..Z_n)``
#### Output
- E  = Constituent Information measures in S
- U = The distinct (unique) elements from E
- V = The individual elements of S
- λ = The scaling coefficients of U (i.e., ``V=λ^T E``)

#### Examples
```julia-repl
julia> S="I(X;Y)-H(X,Y|Z)-3H(X,Y)+2I(X;Y)"
julia> E,UE,V,λ=entropic_terms(S)
["I(X;Y)" "H(X,Y|Z)" "H(X,Y)" "I(X;Y)"],
["I(X;Y)" "H(X,Y|Z)" "H(X,Y)"],
["1I(X;Y)" "-1H(X,Y|Z)" "-3H(X,Y)" "2I(X;Y)"],
[1.0 -1.0 -3.0 2.0]
```
"""
function entropic_terms(S::AbstractString)
	A0= "+" * S
	A1= replace(A0,"+H"=>"+1H","-H"=>"-1H")
	A2= replace(A1,"+I"=>"+1I")
	A2a= replace(A2,"-I"=>"-1I")
	A2b= replace(A2a,"-"=>"+-")
	A3=replace(A2b,"++"=>"+")
	A4=split(A3,"+")
	V=A4[2:end]

	E=[match(r"([H|I]?\(.*)",i).captures[1] for i in V]
	λ=parse.(Float64,[match(r"([-|+]*\d*[.]?(\d*)?)",i).captures[1] for i in V])
	UE=unique(E)

	return E,UE,V,λ
end

"""
    C=minimalE(V,U,λ)

This function is internally used in `simplifyH(S::AbstractString)`.
### Arguments
- ``E``  = Constituent Information measures in ``S``
- ``U`` = The distinct (unique) elements from ``E``
- ``V`` = The individual elements of ``S``
- ``λ`` = The scaling coefficients of ``U`` (i.e., ``V=λ^T E``)

### Output
- C = Simplified Information Expression as a linear combination of ``U``

#### Examples
```julia-repl
julia> S="7H(X,Y|Z1,Z2)+2H(X,Y)-4H(X,Y)+H(Z)-3H(X,Y|Z1,Z2)"
julia> E,U,V,λ=entropic_terms(S)
julia> Z= minimalE(V,U,λ)
"+4.0H(X,Y|Z1,Z2)-2.0H(X,Y)+H(Z)"
```
"""
function minimalE(V,UE,λ)
	Z=""
	for u in UE
		Λ0=sum(λ[occursin.(u, V)])
		Λ = round(Λ0,digits=3)
		if Λ ≠ 0
			Y="+" * string(Λ) * u
		else
			Y=""
		end
		Z=Z*Y
	end
	P=replace(Z," "=>"","+-"=>"-","++"=>"+")
	C=replace(P,"-1.0H"=>"-H","+1.0H"=>"+H","-1.0I"=>"-I","+1.0I"=>"+I")
	return C
end

"""
Express an Information Expression in the simplest form (algebraically).
#### Examples
```julia-repl
julia> simplify("I(X;Y)-H(X,Y|Z)-3H(X,Y)+2I(X;Y)")
"+3.0I(X;Y)-H(X,Y|Z)-3.0H(X,Y)"

julia> simplify("3I(X;Y|Z)-H(X,Y|Z)-3I(X;Y|Z)-3H(X,Y|Z)+H(X1,X2)")
"-4.0H(X,Y|Z)+H(X1,X2)"
```
"""
function simplify(S)
	E,UE,V,λ =	entropic_terms(S)
	M=minimalE(V,UE,λ)
	return M
end

"""
Express an Information Expression in the simplest form (algebraically). This is used in conjunction with the canonical expression.
#### Examples
```julia-repl
julia>
```
"""
function simplifyH(S)
	E,UE,V,λ =	entropic_terms(S)
	M=minimalE(V,UE,λ)
	return M
end


text:
#=##############################################################################
# DESCRIPTION
    Types defining maneuvers of flight vehicles.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Oct 2019
  * License   : MIT
=###############################################################################

################################################################################
# ABSTRACT MANEUVER TYPE
################################################################################

using LinearAlgebra: norm, cross
using Plots

"""
    `AbstractManeuver{N, M}`

`N` indicates the number of tilting systems in this maneuver, while and `M`
indicates the number of rotor systems.

Every implementation of `AbstractManeuver` must have the properties:

 * `angle::NTuple{N, Function}` where `angle[i](t)` returns the angle of the
        i-th tilting system at time `t` (t is nondimensionalized by the total
        time of the maneuver, from 0 to 1, beginning to end).
 * `RPM::NTuple{M, Function}` where `RPM[i](t)` returns the normalized RPM of
        the i-th rotor system at time `t`. This RPM values are normalized by the
        an arbitrary RPM value (usually RPM in hover or cruise).
"""
abstract type AbstractManeuver{N, M} end

##### FUNCTIONS REQUIRED IN IMPLEMENTATIONS ####################################
"""
    `calc_dV(maneuver::AbstractManeuver, vehicle::Vehicle, t, dt, ttot, Vref)`

Returns the change in velocity `dV=[dVx, dVy, dVz]` (m/s) of `vehicle` performing
`maneuver` at time `t` (s) after a time step `dt` (s).  `Vref` and `tot` are the
reference velocity and the total time at which this maneuver is being performed,
respectively. `dV` is in the global reference system.
"""
function calc_dV(self::AbstractManeuver, vehicle::AbstractVehicle, t::Real,
                                            dt::Real, ttot::Real, Vref::Real)
    error("$(typeof(self)) has no implementation yet!")
end

"""
    `calc_dw(maneuver::AbstractManeuver, vehicle::Vehicle, t, dt, Vref, ttot)`

Returns the change in angular velocity `dW=[dWx, dWy, dWz]` (about global
axes, in radians) of `vehicle` performing `maneuver` at time `t` (s) after a
time step `dt` (s). `ttot` is the total time at which this maneuver is to be
performed.
"""
function calc_dW(self::AbstractManeuver, vehicle::AbstractVehicle, t::Real,
                                                        dt::Real, ttot::Real)
    error("$(typeof(self)) has no implementation yet!")
end

##### COMMON FUNCTIONS  ########################################################
"""
    `get_ntltsys(self::AbstractManeuver)`
Return number of tilting systems.
"""
get_ntltsys(self::AbstractManeuver) = typeof(self).parameters[1]

"""
    `get_nrtrsys(self::AbstractManeuver)`
Return number of rotor systems.
"""
get_nrtrsys(self::AbstractManeuver) = typeof(self).parameters[2]

"""
    `get_angle(maneuver::AbstractManeuver, i::Int, t::Real)`

Returns the angle (in degrees) of the i-th tilting system at the non-dimensional
time t.
"""
function get_angle(self::AbstractManeuver, i::Int, t::Real)
    if i<=0 || i>get_ntltsys(self)
        error("Invalid tilting system #$i (max is $(get_ntltsys(self))).")
    end
    if t<0 || t>1
        @warn("Got non-dimensionalized time $(t).")
    end
    return self.angle[i](t)
end

"""
    `get_angles(maneuver::AbstractManeuver, t::Real)`

Returns the angle (in degrees) of every tilting systems at the non-dimensional
time t.
"""
get_angles(self::AbstractManeuver, t::Real) = Tuple(a(t) for a in self.angle)

"""
    `get_RPM(maneuver::AbstractManeuver, i::Int, t::Real)`

Returns the normalized RPM of the i-th rotor system at the non-dimensional time
t.
"""
function get_RPM(self::AbstractManeuver, i::Int, t::Real)
    if i<=0 || i>get_nrtrsys(self)
        error("Invalid rotor system #$i (max is $(get_nrtrsys(self))).")
    end
    if t<0 || t>1
        @warn("Got non-dimensionalized time $(t).")
    end
    return self.RPM[i](t)
end

"""
    `get_RPMs(maneuver::AbstractManeuver, t::Real)`

Returns the normalized RPM of every rotor systems at the non-dimensional time
t.
"""
get_RPMs(self::AbstractManeuver, t::Real) = Tuple(rpm(t) for rpm in self.RPM)

##### COMMON INTERNAL FUNCTIONS  ###############################################

##### END OF ABSTRACT MANEUVER #################################################

################################################################################
# KINEMATIC MANEUVER TYPE
################################################################################
"""
    `KinematicManeuver(angle, RPM, Vvehicle, anglevehicle)`

A vehicle maneuver where the kinematic are prescribed.

# ARGUMENTS
* `angle::NTuple{N, Function}` where `angle[i](t)` returns the angles
        `[Ax, Ay, Az]` (in degrees)of the i-th tilting system at time `t` (t is
        nondimensionalized by the total time of the maneuver, from 0 to 1,
        beginning to end).
* `RPM::NTuple{M, Function}` where `RPM[i](t)` returns the normalized RPM of
        the i-th rotor system at time `t`. This RPM values are normalized by the
        an arbitrary RPM value (usually RPM in hover or cruise).
* `Vvehicle::Function` where `Vvehicle(t)` returns the normalized vehicle
        velocity `[Vx, Vy, Vz]` at the normalized time `t`. Velocity is
        normalized by a reference velocity (typically, cruise velocity).
* `anglevehicle::Function` where `anglevehicle(t)` returns the angles
        `[Ax, Ay, Az]` (in degrees) of the vehicle relative to the global
        coordinate system at the normalized time `t`.
"""
struct KinematicManeuver{N, M} <: AbstractManeuver{N, M}
    angle::NTuple{N, Function}
    RPM::NTuple{M, Function}
    Vvehicle::Function
    anglevehicle::Function
end

# # Implicit N and M constructor
# KinematicManeuver(a::NTuple{N, Function}, b::NTuple{M, Function},
#                     c::Function, d::Function
#                  ) where {N, M} = KinematicManeuver{N, M}(a, b, c, d)

##### FUNCTIONS  ###############################################################
function calc_dV(self::KinematicManeuver, vehicle::AbstractVehicle, t::Real,
                                            dt::Real, ttot::Real, Vref::Real)
    return Vref * (self.Vvehicle((t+dt)/ttot) - self.Vvehicle(t/ttot))
end

function calc_dW(self::KinematicManeuver, vehicle::AbstractVehicle, t::Real,
                                                         dt::Real, ttot::Real)
    prev_W = (self.anglevehicle(t/ttot) - self.anglevehicle((t-dt)/ttot)) / dt
    cur_W = (self.anglevehicle((t+dt)/ttot) - self.anglevehicle(t/ttot)) / dt
    return pi/180 * (cur_W - prev_W)
end

##### INTERNAL FUNCTIONS  ######################################################

##### END OF KINEMATICMANEUVER  ################################################

"""
    MassProp(m, Ixx, Iyy, Izz, Ixz, Ixy, Iyz)

Mass and moments of inertia in the body frame.
Ixx = int(y^2 + z^2, dm)
Ixz = int(xz, dm)

We can model our quadcopter as two thin uniform rods crossed at the origin with a point mass
(motor) at the end of each. With this in mind, it’s clear that the symmetries result in a diagonal
inertia matrix
"""
struct MassProp{TF}
    m::TF
    Ixx::TF
    Iyy::TF
    Izz::TF
end

"""
    ConstantAtmosphere(Wi, Wb, rho, asound, g)

Constant atmospheric properties.
"""
struct Atmosphere{TF,TV}
    Wi::TV
    Wb::TV
    rho::TF
    asound::TF
    g::TF
end

################################################################################
# DYNAMIC MANEUVER TYPE
################################################################################
"""
    `DynamicManeuver(angle, RPM, Vvehicle, anglevehicle)`

A vehicle maneuver where the kinematic are prescribed.

# ARGUMENTS
* `angle::NTuple{N, Function}` where `angle[i](t)` returns the angles
        `[Ax, Ay, Az]` (in degrees)of the i-th tilting system at time `t` (t is
        nondimensionalized by the total time of the maneuver, from 0 to 1,
        beginning to end).
* `RPM::NTuple{M, Function}` where `RPM[i](t)` returns the normalized RPM of
        the i-th rotor system at time `t`. This RPM values are normalized by the
        an arbitrary RPM value (usually RPM in hover or cruise).
* `Vvehicle::Function` where `Vvehicle(t)` returns the normalized vehicle
        velocity `[Vx, Vy, Vz]` at the normalized time `t`. Velocity is
        normalized by a reference velocity (typically, cruise velocity).
* `anglevehicle::Function` where `anglevehicle(t)` returns the angles
        `[Ax, Ay, Az]` (in degrees) of the vehicle relative to the global
        coordinate system at the normalized time `t`.
"""
mutable struct DynamicManeuver{N, M}
    ttot::Real
    angle::NTuple{N, Function}
    RPM::NTuple{M, Function} #maybe change later
    RPMref::Real
    Vvehicle::Function
    anglevehicle::Function
    mp::MassProp{Float64}
    atm::Atmosphere{Float64,Array{Float64,1}}
    rotors::Array{FLOWVLM.Rotor,1} #change to 4 when we use quadrotor?
    k::Real
    b::Real
    kd::Real
    l::Real
    x::Array{Float64,1}
    xprev::Array{Float64,1}
    xdot::Array{Float64,1}
    xdotprev::Array{Float64,1}
    theta::Array{Float64,1}
    thetaprev::Array{Float64,1}
    thetadot::Array{Float64,1}
    omega::Array{Float64,1}
    omegaprev::Array{Float64,1}
    omegadot::Array{Float64,1}
    omegadotprev::Array{Float64,1}
    a::Array{Float64,1}
    tau::Array{Float64,1}
    MovementMatrix::Matrix{Float64}
    nsteps::Int
end

##### INTERNAL FUNCTIONS  ######################################################
function calc_dV(self::DynamicManeuver,
            dt::Real, Vref::Real, t::Real, nt::Int)
self.xdotprev = self.xdot
self.omegaprev = self.omega
self.thetaprev = self.theta

self.a ,self.omegadot, self.omega, self.thetadot,
self.theta, self.xdot, self.x, self.tau, self.MovementMatrix =
equations(dt, self, t, nt)
return (self.xdot - self.xdotprev)
end

"""
  'calc_dW'
"""
function calc_dW(self::DynamicManeuver,
    dt::Real, t::Real, nt::Int)

self.xdotprev = self.xdot
self.omegaprev = self.omega
self.thetaprev = self.theta

self.a ,self.omegadot, self.omega, self.thetadot,
self.theta, self.xdot, self.x, self.tau, self.MovementMatrix = equations(dt, self, t, nt)

return (self.omega- self.omegaprev)
end
##### END OF KINEMATICMANEUVER  ################################################

#


text:
using Canvas
using Markdown

import Canvas: @api, render

desc = md"""
# Hello, World!
You can write any **Markdown** text and place it in the middle of other
tiles.
""" |> pad(1inch)

main(window) = begin
    push!(window.assets, "layout2")
    push!(window.assets, "widgets")

    t, p = wire(tabs([hbox(icon("home"), hskip(1em), "Home"),
                      hbox(icon("info-outline"), hskip(1em),  "Notifications"),
                      hbox(icon("settings"), hskip(1em), "Settings")]),
                pages([desc, "b", "c"]), :lol, :selected)

    vbox(toolbar([iconbutton("face"), "My App", flex(), iconbutton("search")]),
         maxwidth(600px, t),
         pad(1em, p))
end


text:
using UncertaintyQuantification, Random

Random.seed!(2547);
N = 100
a = 2
b = 1

changepoint = floor(Int, N*rand(Uniform(0,1)))

lambda = Gamma(a, 1/b)

lambdas = zeros(N)

for i in 1:N
    if i<changepoint
        lambdas[i] = rand(lambda) *changepoint
    else
        lambdas[i] = rand(lambda) *(N-changepoint)
    end
end

sampledata = [rand(Poisson(i)) for i in lambdas]
Random.seed!()

x = sampledata

lambda1(in) = rand(Gamma(a + sum(x[1:changepoint-1]), 1/(changepoint + b)))
lambda2(in) = rand(Gamma(a + sum(x[changepoint:N]), 1/(N - changepoint + b)))

function multinomialN(in)
    propn = zeros(N)
    for i in 1:N
    propn[i] = sum(x[1:i-1])*log(in[1]) - i*in[1] + sum(x[i:N])*log(in[2]) - (N-changepoint)*in[2]
    end
    propn = exp.(propn .- maximum(propn))
    multiN = rand(Multinomial(1, propn./sum(propn)))
    rv = 0
    for i in 1:N
        if multiN[i]==1 rv = i end
    end
    return rv
end

condist = [lambda1, lambda2, multinomialN]

rsample = gibbssample(condist, [3, 6, 50], 10000, 200)


text:
module OmegaOptim

import OmegaCore: argmax

using Optim

"Simple initialization algorithm"
struct RandInit end

randinit(rng, ::BinaryChoice{Bool}) = rand(rng, Bool)
randinit(rng, ::BinaryChoice{T})  where {T <: Integer} = rand(rng, (zero(T), one(T)))
randinit(rng, ::Unit{T}) = rand(rng, T)

function randindepinit(rng, x)
  ω = defω()
end

function init(x)
end

function OmegaCore.argmax(ℓ, ω0, ::Optim.AbstractOptimizer)
  ## Step 1

end

function test()
  x = 1 ~ Interval(0, 10)
  y = 2 ~ Interval(0, 10)
  ℓ = x +ₚ y
  ω0 = LazyΩ(x => 0.3, y => 0.9)
  optimize(ℓ, ω0, )
end

## What's needed?

## Initialization
  # Like lazy Omega
  # Or execute the program
## Overload optimize or use different name?
## function linearize(x) = xs::AbstractArray -> genω(xs)

end # module

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
################################################################################

@testset "Reset" begin

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

  # load small array into hidden Markov model
  include("data/smallArray.jl")

  # reset
  resetedhmm = HiddenMarkovModelReaders.reset!(HMM([zeros(size(x, 2))], [zeros(size(x, 1) + 1)], Array{Int32}(undef, 0)))

  # testing
  @test resetedhmm == [[0.0; repeat([-1.0], size(x, 1))]]
  @test hmm.data == [zeros(size(x, 2))]
  @test hmm.model == [[0.0; repeat([-1.0], size(x, 1))]]

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

  # read file into hidden Markov model
  include("data/readfile.jl")

  # reset
  resetedhmm = HiddenMarkovModelReaders.reset!(HMM([zeros(size(v, 2))], [zeros(size(v, 1) + 1)], Array{Int32}(undef, 0)))

  # testing
  @test resetedhmm == [[0.0; repeat([-1.0], size(v, 1))]]
  @test hmm.data == [zeros(size(v, 2))]
  @test hmm.model == [[0.0; repeat([-1.0], size(v, 1))]]

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

end

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


text:
using HttpServer
using WebSockets
using Merlin
using JukaiNLP
using JukaiNLP.TokenizationSentence
using JukaiNLP.Tagging
using JLD
using JSON
using PyCall

unshift!(PyVector(pyimport("sys")["path"]), "mitie-v0.2-python-2.7")
unshift!(PyVector(pyimport("sys")["path"]), "entity-disambi")
@pyimport NER
is_linux() && (@pyimport mymodule)

type Token
    form::String
    cat::String
    ne::Vector{Any}
    dep::Int
end

Token(form::String) = Token(form, "", [], 0)

function todict(t::Token)
    Dict("form"=>t.form, "cat"=>t.cat, "ne"=>t.ne)
end

const samples = [
    Token("Pierre", "NNP", [], 2),
    Token("Vinken", "NNP", [], 8),
    Token(",", ",", [], 2),
    Token("61", "CD", [], 5),
    Token("years", "NNS", [], 6),
    Token("old", "JJ", [], 2),
    Token(",", ",", [], 2),
    Token("will", "MD", [], 0),
    Token("join", "VB", [], 8),
    Token("the", "DT", [], 11),
    Token("board", "NN", [], 9),
    Token("Nov.", "NNP", [], 9),
    Token("29", "CD", [], 12),
    Token(".", ".", [], 8)
]
const samples2 = [ "Diamond", "Shamrock", "Offshore", "Partners","said","it",
    "had","discovered","gas","offshore","Louisiana", "."]

function load_tokenizer()
    println("loading tokenizer...")
    path = joinpath(Pkg.dir("JukaiNLP"), "web/models")
    model = h5load("$(path)/tokenizer.h5")
    dict = JLD.load("$(path)/tokenizer.jld", "tokenizer")
    Tokenizer(dict, model)
end
const tokenizer = load_tokenizer()

function load_postagger()
    println("loading postagger...")
    path = joinpath(Pkg.dir("JukaiNLP"), "web/models")
    m = h5load("$(path)/postagger.h5")
    model = Tagging.POSModel(m["wordfun"],m["charfun"],m["sentfun"])
    t = JLD.load("$(path)/postagger.jld", "postagger")
    t.model = model
    t
end
const postagger = load_postagger()

function setpostag!(tokens::Vector{Token}, forms::Vector{String})
    tags = postagger(forms)
    for i = 1:length(tokens)
        tokens[i].cat = tags[i]
    end
end

function setner!(tokens::Vector{Token}, forms::Vector{String})
    ners = NER.predict(forms)
    for (range,tag) in ners
        tokens[start(range)+1].ne = Any[tag,length(range)]
        if is_linux()
            #target = forms[start(range)+1:last(range)+1]
            #mymodule.predict(target, forms)
        end
    end
end

include("conf.jl")

const filepath = dirname(@__FILE__)
const clients = Dict()
const conf = begin
    e = readconf(joinpath(filepath,"conf/visual.conf"))
    d = Dict("entity_types" => e)
    #=
    "relation_types" => [
        Dict(
            "type" => "Anaphora",
            "labels" => ["Anaphora", "Ana"],
            "dashArray" => "3,3",
            "color" => "purple",
            "args" => [
                Dict("role" => "Anaphor", "targets" => ["Person"]),
                Dict("role" => "Entity", "targets" => ["Person"])
            ]
        )
    ]
    =#
    JSON.json(d)
end

wsh = WebSocketHandler() do req, client
    println("Client: $(client.id) is connected.")
    write(client, conf)
    while true
        #println("Request from $(client.id) recieved.")
        msg = bytestring(read(client))
        length(msg) > 1000 && continue
        chars = Vector{Char}(msg)

        doc = tokenizer(chars)
        doc = map(doc) do sent
            tokens = map(w -> Token(w), sent)
            setpostag!(tokens, sent)
            setner!(tokens, sent)
            map(todict, tokens)
        end
        res = JSON.json(doc)
        write(client, res)
    end
end

onepage = readstring(joinpath(dirname(@__FILE__), "index.html"))
httph = HttpHandler() do req::Request, res::Response
    Response(onepage)
end
server = Server(httph, wsh)
run(server, 3000)

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
"""
    import_abundance_table(file::AbstractString; delim::Char='\t')

Given a file path or paths to abundance tables (eg humann2 or metaphlan2),
create abundance table. Table is presumed to have samples in columns and
features in rows. First column is taken as feature IDs.
"""
function import_abundance_table(file::AbstractString; delim::Char='\t')
    @info "Importing abundance table" file
    df = CSV.File(file, delim=delim) |> DataFrame
    rename!(df, names(df)[1] => :col1)
    return df
end

function import_abundance_tables(files::Array{<:AbstractString, 1}; delim::Char='\t')
    @info "Importing abundance tables"
    fulltable = DataFrame(col1=String[])
    for t in files
        df = import_abundance_table(t, delim=delim)
        fulltable = join(fulltable, df, on=:col1, kind=:outer)
    end

    # replace all missing values (from joins) with 0.
    fulltable = map(c -> eltype(c) <: Union{<:Number, Missing} ? collect(Missings.replace(c, 0)) : c, eachcol(fulltable))
    return fulltable
end

function clean_abundance_tables(files::Array{String, 1};
                        delim::Char='\t',
                        col1::Symbol=:taxon,
                        suffix::String="_taxonomic_profile")
    @info "Cleaning"
    t = import_abundance_tables(files, delim=delim)

    rename!(t, :col1 => col1)
    rename!(n-> Symbol(replace(String(n), suffix => "")), t)

    return t
end

function rm_strat!(df::DataFrame; col::Union{Int, Symbol}=1)
    filter!(row->!occursin(r"\|", row[1]), df)
end

"""
    permanova(dm::Array{<:Real,2}, metadata::AbstractVector, nperm::Int=999;
                label=nothing, datafilter=x->true)
    permanova(dm::Array{<:Real,2}, metadata::AbstractDataFrame, nperm=999;
                fields=names(metadata), kwargs...)

Performs PERMANOVA analysis from R's [`vegan`](https://www.rdocumentation.org/packages/vegan/versions/2.4-2) package
using the `adonis` function.

**Positional arguments**:

- `dm`: a symetric distance matrix.
- `metadata`: either a vector of numerical or categorical data to test against,
  or a DataFrame with columns for each variable to test against.
  Any missing data in the vector or rows of the DataFrame with missing data
  will be filtered out.
- `nperm`=999: number of permutations for PERMANOVA.

**Keyword Arguments**:

- `datafilter=x-> true`: a function to filter elements (or rows) of `metadata`.
  Removal of missing values occurs before this function is applied.
- `label=nothing`: If provided, adds a column `label` to the results
  filled with this value.
  Useful if performing multiple runs that will be combined in a single DataFrame.
- `fields`: if passing a DataFrame as `metadata`,
  an array of symbols may be passed to select only certain columns
  and/or determine their order for the resulting PERMANOVA.

Note: this will throw an error if `vegan` is not installed.
To install:

```julia
using RCall

reval("install.packages('vegan')")
```
"""
function permanova(dm::Array{<:Real,2}, metadata::AbstractVector, nperm::Int=999;
                    datafilter=x->true, label=nothing)
    size(dm,1) != size(dm,2) && throw(ArgumentError("dm must be symetrical distance matrix"))
    size(dm,2) != length(metadata) && throw(ArgumentError("Metadata does not match the size of distance matrix"))
    let notmissing = map(!ismissing, metadata)
        metadata = metadata[notmissing]
        dm = dm[notmissing, notmissing]
    end

    filt = map(datafilter, metadata)
    r_dm = dm[filt, filt]
    metadata = metadata[filt]
    @rput r_dm
    @rput metadata
    reval("library(vegan)")

    reval("p <- adonis(r_dm ~ metadata, permutations = $nperm)")

    @rget p

    p = p[:aov_tab]
    if !isnothing(label)
        p[!,:label] = fill(label, size(p, 1))
    end

    return p
end

function permanova(dm::Array{<:Real,2}, metadata::AbstractDataFrame, nperm=999;
            datafilter=x->true,
            label=nothing,
            fields=names(metadata))

    let notmissing = map(row->all(!ismissing, row[fields]), eachrow(metadata))
        metadata = metadata[notmissing, :]
        dm = dm[notmissing, notmissing]
    end

    fields = join(String.(fields), " + ")

    filt = map(datafilter, eachrow(metadata))
    r_meta = metadata[filt, :]
    r_dm = dm[filt,filt]
    @rput r_meta
    @rput r_dm
    reval("library(vegan)")
    reval("p <- adonis(r_dm ~ $fields, data=r_meta, permutations = $nperm)")

    @rget p
    p = p[:aov_tab]
    if !isnothing(label)
        p[!,:label] = fill(label, size(p, 1))
    end

    return p
end
#==============
PanPhlAn Utils
==============#

function panphlan_calcs(df::DataFrame)
    abun = abundancetable(df)
    dm = pairwise(Jaccard(), abun, dims=2)
    rowdm = pairwise(Jaccard(), abun, dims=1)
    col_clust = hclust(dm, :single)
    row_clust = hclust(rowdm, :single)
    optimalorder!(col_clust, dm)
    optimalorder!(row_clust, rowdm)

    mds = fit(MDS, dm, distances=true)

    return abun, dm, col_clust, row_clust, mds
end


text:
__precompile__()

module Convex
using OrderedCollections: OrderedDict
using LinearAlgebra
using SparseArrays
using LDLFactorizations
using AbstractTrees: AbstractTrees

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

# Functions
export conv,
    dotsort, entropy, exp, geomean, hinge_loss, huber, inner_product, invpos
export log_perspective,
    logisticloss,
    logsumexp,
    matrixfrac,
    neg,
    norm2,
    norm_1,
    norm_inf,
    nuclearnorm
export partialtrace,
    partialtranspose, pos, qol_elementwise, quadform, quadoverlin, rationalnorm
export relative_entropy,
    scaledgeomean,
    sigmamax,
    square,
    sumlargest,
    sumlargesteigs,
    sumsmallest,
    sumsquares
export GeomMeanHypoCone,
    GeomMeanEpiCone,
    RelativeEntropyEpiCone,
    quantum_relative_entropy,
    quantum_entropy
export trace_logm, trace_mpower, lieb_ando

# rexports from LinearAlgebra
export diag, diagm, Diagonal, dot, eigmax, eigmin, kron, logdet, norm, tr

# Constraints
export Constraint
export isposdef, ⪰, ⪯ # PSD constraints
export socp
export Constraint # useful for making abstractly-typed vectors via `Constraint[]`

# Variables
export Constant, ComplexVariable, HermitianSemidefinite, Semidefinite, Variable
export curvature, evaluate, fix!, free!, monotonicity, sign, vexity
export BinVar, IntVar, ContVar, vartype, vartype!
export constraints, add_constraint!, set_value!, evaluate

# Signs
export Positive, Negative, ComplexSign, NoSign

# Problems
export add_constraints!, maximize, minimize, Problem, satisfy, solve!

# Module level globals

"""
    emit_dcp_warnings

Controls whether or not warnings are emitted for when an expression fails to be
of disciplined convex form. To turn warnings off, override the method via

    Convex.emit_dcp_warnings() = false

This will cause Julia's method invalidation to recompile any functions emitting
DCP warnings and remove them. This should be run from top-level (not within a function).
"""
emit_dcp_warnings() = true

"""
    MAXDEPTH

Controls depth of tree printing globally for Convex.jl; defaults to 3. Set via

    Convex.MAXDEPTH[] = 5
"""
const MAXDEPTH = Ref(3)

"""
    MAXWIDTH

Controls width of tree printing globally for Convex.jl; defaults to 15. Set via

    Convex.MAXWIDTH[] = 15
"""
const MAXWIDTH = Ref(15)

"""
    MAXDIGITS

When priting IDs of variables, only show the initial and final digits
if the full ID has more than double the number of digits specified
here.  So, with the default setting MAXDIGITS=3, any ID longer than 7
digits would be shortened; for example, ID `14656210999710729289`
would be printed as `146…289`.

This setting controls tree printing globally for Convex.jl; defaults to 3.

Set via:

    Convex.MAXDIGITS[] = 3
"""
const MAXDIGITS = Ref(3)

### modeling framework
include("dcp.jl")
include("expressions.jl")
# need to define `Variable` before `UniqueConicForms`
include("variable.jl")
include("conic_form.jl")
# need to define `conic_form!` for `Variable`s after `UniqueConicForms`
include("variable_conic_form.jl")
include("constant.jl")
include("constraints/constraints.jl")
include("constraints/signs_and_sets.jl")
include("constraints/soc_constraints.jl")
include("constraints/exp_constraints.jl")
include("constraints/sdp_constraints.jl")
include("problems.jl")
include("solution.jl")

### affine atoms
include("atoms/affine/add_subtract.jl")
include("atoms/affine/multiply_divide.jl")
include("atoms/affine/sum.jl")
include("atoms/affine/transpose.jl")
include("atoms/affine/index.jl")
include("atoms/affine/diag.jl")
include("atoms/affine/diagm.jl")
include("atoms/affine/stack.jl")
include("atoms/affine/dot.jl")
include("atoms/affine/reshape.jl")
include("atoms/affine/trace.jl")
include("atoms/affine/partialtrace.jl")
include("atoms/affine/partialtranspose.jl")
include("atoms/affine/conv.jl")
include("atoms/affine/real_imag.jl")
include("atoms/affine/inner_product.jl")
include("atoms/affine/conjugate.jl")
include("atoms/affine/kron.jl")

### lp atoms
include("atoms/lp_cone/abs.jl")
include("atoms/lp_cone/maximum.jl")
include("atoms/lp_cone/minimum.jl")
include("atoms/lp_cone/max.jl")
include("atoms/lp_cone/min.jl")
include("atoms/lp_cone/sumlargest.jl")
include("atoms/lp_cone/dotsort.jl")

### SOC atoms
include("atoms/second_order_cone/norm.jl") # also includes some lp atoms
include("atoms/second_order_cone/norm2.jl")
include("atoms/second_order_cone/quadoverlin.jl")
include("atoms/second_order_cone/qol_elementwise.jl")
include("atoms/second_order_cone/geomean.jl")
include("atoms/second_order_cone/quadform.jl")
include("atoms/second_order_cone/power_to_socp.jl")
include("atoms/second_order_cone/rationalnorm.jl")
include("atoms/second_order_cone/huber.jl")

### SDP atoms
include("atoms/sdp_cone/nuclearnorm.jl")
include("atoms/sdp_cone/operatornorm.jl")
include("atoms/sdp_cone/eig_min_max.jl")
include("atoms/sdp_cone/matrixfrac.jl")
include("atoms/sdp_cone/sumlargesteigs.jl")
include("atoms/sdp_cone/geom_mean_hypocone.jl")
include("atoms/sdp_cone/geom_mean_epicone.jl")
include("atoms/sdp_cone/relative_entropy_epicone.jl")
include("atoms/sdp_cone/quantum_relative_entropy.jl")
include("atoms/sdp_cone/quantum_entropy.jl")
include("atoms/sdp_cone/trace_logm.jl")
include("atoms/sdp_cone/trace_mpower.jl")
include("atoms/sdp_cone/lieb_ando.jl")

### exponential atoms
include("atoms/exp_cone/exp.jl")
include("atoms/exp_cone/log.jl")
include("atoms/exp_cone/logsumexp.jl")
include("atoms/exp_cone/entropy.jl")
include("atoms/exp_cone/relative_entropy.jl")

### exp + sdp atoms
include("atoms/exp_+_sdp_cone/logdet.jl")

### utilities
include("utilities/tree_print.jl")
include("utilities/tree_interface.jl")
include("utilities/show.jl")
include("utilities/iteration.jl")
include("utilities/broadcast.jl")
include("problem_depot/problem_depot.jl")

end


text:
# Number of points used in rand testing.
nrand = 100

@testset "Interval" begin
    @test dim(Interval(rand())) == 1
    @test codim(Interval(rand())) == 0
    @test ambientdim(Interval(rand())) == 1

    @testset "Basics" begin
        for _ in 1:nruns
            len = 5rand() + 1
            t   = rand()
            i1  = Interval(len)
            val = i1(t)
            @test length(val) == 1
            @test val[1] ≈ t * len

            i2  = Interval(x -> x^2)
            val = i2(t, scale = len)
            @test val[1] ≈ t * len^2
        end
    end

    @testset "rand" begin
        for _ in 1:nruns
            len = 5rand() + 1
            i1  = Interval(len)
            @test all(0 .< rand(i1, nrand) .< len)

            i2  = Interval(x -> x^2)
            @test all(0 .< rand(i1, nrand, scale = len) .< len^2)
        end
    end
end

@testset "NSphere" begin
    @testset "Basics" begin
        for d in 1:nruns
            rad = 5rand() + 1
            sp1 = NSphere{d}(rad)
            val = sp1(fill(rand(), d)...)
            @test length(val) == d + 1
            @test norm(val) ≈ rad

            sp2 = NSphere{d}(x -> x^2)
            val = sp2(fill(rand(), d)..., scale = rad)
            @test norm(val) ≈ rad^2

            @test dim(sp1) == d
            @test codim(sp1) == 1
            @test ambientdim(sp1) == d + 1
        end
    end

    @testset "rand" begin
        for d in 1:nruns
            rad = 5rand() + 1
            sp1 = NSphere{d}(rad)
            @test all(mapslices(norm, rand(sp1, nrand), 1) .≈ rad)

            sp2 = NSphere{d}(x -> x^2)
            @test all(mapslices(norm, rand(sp2, nrand, scale = rad), 1) .≈ rad^2)
        end
    end
end

@testset "Knot" begin
    # TODO: How to test this????
    @testset "Basics" begin
        for _ in 1:nruns
        end
    end
end

@testset "UnitSpace, copy" begin
    u = UnitSpace()
    for space in [Interval(),
                  Knot(),
                  NSphere{10}(),
                  Circle() * Circle(),
                  Circle() × Circle(),
                  Interval()^2 * (Knot() * NSphere{3}() × Interval() × Sphere()),
                  UnitSpace()]
        T = typeof(space)
        @test space == u * space == space * u == u * space * u
        @test space == u × space == space × u == u × space × u
        @test u == one(space) == one(T)

        @test copy(space) ≡ space
    end
end

@testset "ProductSpace" begin
    # Do Torus, Cylinder, Interval^n

    @testset "Torus" begin
        torus31 = Circle(3.0) * Circle()
        circ2   = Circle(2.0)
        circ3   = Circle(3.0)
        circ4   = Circle(4.0)

        for _ in 1:nruns
            t = rand()
            @test torus31(t, 0)    ≈ [circ2(t); 0] atol = 1e-16
            @test torus31(t, 0.5)  ≈ [circ4(t); 0] atol = 1e-16
            @test torus31(t, 0.25) ≈ [circ3(t); 1] atol = 1e-16
            @test torus31(t, 0.75) ≈ [circ3(t);-1] atol = 1e-16
        end
    end

    @testset "Cylinder" begin
        cyl  = Interval() * Circle()
        circ = Circle()

        for _ in 1:nruns
            t1, t2 = rand(2)
            @test cyl(t1, t2) ≈ [t1; circ(t2)]
        end
    end

    @testset "Ring" begin
        ring  = Circle(3) * Interval()
        circ2 = Circle(2)
        circ3 = Circle(3)

        for _ in 1:nruns
            t = rand()
            @test ring(t, 1) ≈ circ2(t)
            @test ring(t, 0) ≈ circ3(t)
        end
    end

    @testset "Interval^n" begin
        for d in 1 + (1:nruns)
            lens = 5rand(d) + 1
            intd = prod(Interval(l) for l in lens)

            @test typeof(intd) == ProductSpace{d, 0, d}
            ts = rand(d)
            @test intd(ts...) ≈ ts .* lens
        end
    end

    @testset "Circle^n" begin
        for n in 1 + (1:nruns)
            mf = Circle()^n
            @test mf(fill(.25, n)...) ≈ vcat(0, fill(1, n))  atol = 1e-16
            @test mf(fill(.75, n)...) ≈ vcat(0, fill(-1, n)) atol = 1e-16
        end
    end

    @testset "tnbframe errors" begin
        # No error.
        circxglome = Circle() * NSphere{3}()
        @test circxglome            ≠ nothing
        @test Circle()^5            ≠ nothing
        @test Sphere()^5            ≠ nothing
        @test Interval()^5          ≠ nothing
        @test Circle() * circxglome ≠ nothing
        # Error / tnbframe not defined.
        @test_throws ErrorException NSphere{3}() * Circle()
        @test_throws ErrorException circxglome * Circle()
        @test_throws ErrorException circxglome * circxglome
    end
end

@testset "CartesianSpace" begin

    @testset "Clifford Torus" begin
        r1, r2 = 5rand(2) + 1
        clifft = Circle(r1) × Circle(r2)
        @test typeof(clifft) == CartesianSpace{2, 2, 2}

        for _ in 1:nruns
            t1, t2 = rand(2)
            @test clifft(t1, t2) ≈ [r1*cos(2π*t1), r1*sin(2π*t1),
                                    r2*cos(2π*t2), r2*sin(2π*t2)]
        end
    end

    @testset "Cylinder" begin
        cyl  = Interval() × Circle()
        circ = Circle()

        for _ in 1:nruns
            t1, t2 = rand(2)
            @test cyl(t1, t2) ≈ [t1; circ(t2)]
        end
    end

    @testset "Interval^n" begin
        for d in 1 + (1:nruns)
            lens = 5rand(d) + 1
            intd = reduce(×, Interval(l) for l in lens)

            @test typeof(intd) == CartesianSpace{d, 0, d}
            ts = rand(d)
            @test intd(ts...) ≈ ts .* lens
        end
    end
end

@testset "when * == ×" begin
    for n in 1:nruns
        t, u, v = rand(3)
        @test (Interval() * Interval())(t, u)  ≈ (Interval() × Interval())(t, u)
        @test (Interval() * Circle())(t, u)    ≈ (Interval() × Circle())(t, u)
        @test (Interval() * Sphere())(t, u, v) ≈ (Interval() × Sphere())(t, u, v)

        prod_space = reduce(*, fill(Interval(), n)) * NSphere{n}()
        cart_space = reduce(×, fill(Interval(), n)) × NSphere{n}()
        inputs     = rand(n + n)

        @test prod_space(inputs...) ≈ cart_space(inputs...)
    end
end


text:
module DemultiplexFASTQ

# Write your package code here.

using CSV, CodecZlib
include("demultiplex.jl")
include("run.jl")

export processargs

end


text:
import Base.Sort
immutable BubbleSortAlg <: Sort.Algorithm end
const BubbleSort = BubbleSortAlg()

function Base.sort!(v::AbstractVector, lo::Int, hi::Int, ::BubbleSortAlg, o::Sort.Ordering)
while true
    clean = true
        for i = lo:hi-1
             if Sort.lt(o, v[i+1], v[i])
                v[i+1], v[i] = v[i], v[i+1]
                clean = false
            end
        end
        clean && break
    end
    return v
end

max_stars_repo_name: kwcto/atom-script

max_issues_repo_name: kwcto/atom-script

max_forks_repo_name: kwcto/atom-script


text:
function test_jump_soc(solver)
    @testset "Second-order Cone Programming" begin
        @testset "SOC1" begin
            m = Model(solver=solver)
            @variable(m, x)
            @variable(m, y)
            @variable(m, t >= 0)
            @objective(m, Min, t)
            @constraint(m, x + y >= 1)
            @constraint(m, [t,x,y] in MOI.SecondOrderCone(3))

            JuMP.solve(m)

            @test JuMP.isattached(m)
            @test JuMP.hasvariableresult(m)
            @test JuMP.terminationstatus(m) == MOI.Success
            @test JuMP.primalstatus(m) == MOI.FeasiblePoint

            @test JuMP.objectivevalue(m) ≈ sqrt(1/2) atol=1e-6

            @test JuMP.resultvalue.([x,y,t]) ≈ [0.5,0.5,sqrt(1/2)] atol=1e-3

        end

        @testset "RotatedSOC1" begin
            m = Model(solver=solver)

            @variable(m, x[1:5] >= 0)
            @variable(m, 0 <= u <= 5)
            @variable(m, v)
            @variable(m, t1 == 1)
            @variable(m, t2 == 1)

            @objective(m, Max, v)

            @constraint(m, [t1/sqrt(2),t2/sqrt(2),x...] in MOI.RotatedSecondOrderCone(7))
            @constraint(m, [x[1]/sqrt(2), u/sqrt(2), v] in MOI.RotatedSecondOrderCone(3))

            JuMP.solve(m)

            @test JuMP.isattached(m)
            @test JuMP.hasvariableresult(m)
            @test JuMP.terminationstatus(m) == MOI.Success
            @test JuMP.primalstatus(m) == MOI.FeasiblePoint

            @test JuMP.resultvalue.(x) ≈ [1,0,0,0,0] atol=1e-2
            @test JuMP.resultvalue(u) ≈ 5 atol=1e-4
            @test JuMP.resultvalue(v) ≈ sqrt(5) atol=1e-6
        end
    end
end


text:
using Base.Test

using Perf
@test 1 == 2


text:
abstract type AbstractMechanism{T,Nn,Nb,Ne,Ni} end

"""
$(TYPEDEF)

A `Mechanism` contains the [`Origin`](@ref), [`Body`](@ref)s, and [`EqualityConstraint`](@ref)s of a system and can be used for simulation.
# Important attributes
* `origin`:        The origin of a mechanism.
* `bodies`:        The bodies of a mechanism (Dict).
* `eqconstraints`: The equality constraints (joints) of a mechanism (Dict).
* `Δt`:            The time step of the mechanism.
* `g`:             The gravitational acceleration in z-direction.

# Constuctors
    Mechanism(origin, bodies; Δt, g)
    Mechanism(origin, bodies, eqcs; Δt, g)
    Mechanism(urdf_filename; floating, Δt, g)
"""
mutable struct Mechanism{T,Nn,Nb,Ne,Ni} <: AbstractMechanism{T,Nn,Nb,Ne,Ni}
    origin::Origin{T}
    bodies::UnitDict{Base.OneTo{Int64},Body{T}}
    eqconstraints::UnitDict{UnitRange{Int64},<:EqualityConstraint{T}}
    ineqconstraints::UnitDict{UnitRange{Int64},<:InequalityConstraint{T}}

    graph::Graph{Nn}
    ldu::SparseLDU{T}

    # TODO remove once EqualityConstraint is homogenous
    normf::T
    normΔs::T

    Δt::T
    g::T

    α::T
    μ::T

    function Mechanism(origin::Origin{T},bodies::Vector{Body{T}},
            eqcs::Vector{<:EqualityConstraint{T}}, ineqcs::Vector{<:InequalityConstraint{T}};
            Δt::Real = .01, g::Real = -9.81
        ) where T

        resetGlobalID()
        order = getGlobalOrder()

        for body in bodies
            if norm(body.m) == 0 || norm(body.J) == 0
                @info "Potentially bad inertial properties detected"
            end
        end

        Nb = length(bodies)
        Ne = length(eqcs)
        Ni = length(ineqcs)
        Nn = Nb + Ne

        if Nb < Ne
            @info "More constraints than bodies. Potentially bad behavior."
        end

        currentid = 1

        bdict = Dict{Int64,Int64}()
        for (ind, body) in enumerate(bodies)
            initknotpoints!(body.state, order)

            for eqc in eqcs
                eqc.parentid == body.id && (eqc.parentid = currentid)
                for (ind, bodyid) in enumerate(eqc.childids)
                    if bodyid == body.id
                        eqc.childids = setindex(eqc.childids, currentid, ind)
                        eqc.constraints[ind].childid = currentid
                    end
                end
            end

            for ineqc in ineqcs
                ineqc.parentid == body.id && (ineqc.parentid = currentid)
            end

            body.id = currentid
            currentid += 1

            bdict[body.id] = ind
        end

        eqdict = Dict{Int64,Int64}()
        for (ind, eqc) in enumerate(eqcs)
            eqc.id = currentid
            currentid += 1

            eqdict[eqc.id] = ind
        end

        ineqdict = Dict{Int64,Int64}()
        for (ind, ineqc) in enumerate(ineqcs)
            ineqc.id = currentid
            currentid += 1

            ineqdict[ineqc.id] = ind
        end

        normf = 0
        normΔs = 0

        graph = Graph(origin, bodies, eqcs, ineqcs)
        ldu = SparseLDU(graph, bodies, eqcs, ineqcs, bdict, eqdict, ineqdict)

        # storage = Storage{T}(steps, Nb, Ne)

        bodies = UnitDict(bodies)
        eqcs = UnitDict((eqcs[1].id):(eqcs[Ne].id), eqcs)
        if Ni > 0
            ineqcs = UnitDict((ineqcs[1].id):(ineqcs[Ni].id), ineqcs)
        else
            ineqcs = UnitDict(0:0, ineqcs)
        end

        α = 1
        μ = 1

        new{T,Nn,Nb,Ne,Ni}(origin, bodies, eqcs, ineqcs, graph, ldu, normf, normΔs, Δt, g, α, μ)
    end

    function Mechanism(origin::Origin{T},bodies::Vector{Body{T}},eqcs::Vector{<:EqualityConstraint{T}};
            Δt::Real = .01, g::Real = -9.81
        ) where T

        ineqcs = InequalityConstraint{T}[]
        return Mechanism(origin, bodies, eqcs, ineqcs, Δt = Δt, g = g)
    end

    function Mechanism(origin::Origin{T},bodies::Vector{Body{T}},ineqcs::Vector{<:InequalityConstraint{T}};
            Δt::Real = .01, g::Real = -9.81
        ) where T

        eqc = EqualityConstraint{T}[]
        for body in bodies
            push!(eqc, EqualityConstraint(Floating(origin, body)))
        end
        return Mechanism(origin, bodies, eqc, ineqcs, Δt = Δt, g = g)
    end

    function Mechanism(origin::Origin{T},bodies::Vector{Body{T}};
            Δt::Real = .01, g::Real = -9.81
        ) where T

        eqc = EqualityConstraint{T}[]
        for body in bodies
            push!(eqc, EqualityConstraint(Floating(origin, body)))
        end
        return Mechanism(origin, bodies, eqc, Δt = Δt, g = g)
    end

    function Mechanism(filename::AbstractString; floating::Bool=false, type::Type{T} = Float64, Δt::Real = .01, g::Real = -9.81) where T
        origin, links, joints, loopjoints = parse_urdf(filename, floating, T)
        mechanism = Mechanism(origin, links, [joints;loopjoints], Δt = Δt, g = g)
        set_parsed_values!(mechanism, loopjoints)

        return mechanism
    end
end

mutable struct LinearMechanism{T,Nn,Nb,Ne,Ni} <: AbstractMechanism{T,Nn,Nb,Ne,Ni}
    ## Mechanism attributes
    origin::Origin{T}
    bodies::UnitDict{Base.OneTo{Int64},Body{T}}
    eqconstraints::UnitDict{UnitRange{Int64},<:EqualityConstraint{T}}
    ineqconstraints::UnitDict{UnitRange{Int64},<:InequalityConstraint{T}}

    graph::Graph{Nn}
    ldu::SparseLDU{T}

    # TODO remove once EqualityConstraint is homogenous
    normf::T
    normΔs::T

    Δt::T
    g::T

    α::T
    μ::T

    ## LinearMechanism attributes
    A::AbstractMatrix{T}
    Bu::AbstractMatrix{T}
    Bλ::AbstractMatrix{T}
    G::AbstractMatrix{T}

    xd::Vector{SVector{3,T}}
    vd::Vector{SVector{3,T}}
    qd::Vector{UnitQuaternion{T}}
    ωd::Vector{SVector{3,T}}
    # Fτd::Vector{SVector{3,T}}

    z::Vector{T}
    zsol::Vector{Vector{T}}
    Δz::Vector{T}
    λ::Vector{T}
    λsol::Vector{Vector{T}}
    Δλ::Vector{T}
    u::Vector{T}

    function LinearMechanism(mechanism::Mechanism{T,Nn,Nb,Ne,Ni}, xd, vd, qd, ωd, Fτd, eqcids) where {T,Nn,Nb,Ne,Ni}

        A, Bu, Bλ, G = linearsystem(mechanism, xd, vd, qd, ωd, Fτd, getid.(mechanism.bodies), eqcids)

        z = zeros(T,Nb*12)
        zsol = [zeros(T,Nb*12) for i=1:2]
        Δz = zeros(T,Nb*12)

        nc = 0
        for eqc in mechanism.eqconstraints
            nc += length(eqc)
        end
        λ = zeros(T,nc)
        λsol = [zeros(T,nc) for i=1:2]
        Δλ = zeros(T,nc)

        u = zeros(T,size(Bu)[2])

        new{T,Nn,Nb,Ne,Ni}([getfield(mechanism,i) for i=1:getfieldnumber(mechanism)]..., A, Bu, Bλ, G, xd, vd, qd, ωd, z, zsol, Δz, λ, λsol, Δλ, u)
    end

    function LinearMechanism(mechanism::Mechanism{T,Nn,Nb,Ne,Ni};
        xd = [mechanism.bodies[i].state.xc for i=1:Nb],
        vd = [mechanism.bodies[i].state.vc for i=1:Nb],
        qd = [mechanism.bodies[i].state.qc for i=1:Nb],
        ωd = [mechanism.bodies[i].state.ωc for i=1:Nb],
        eqcids = [],
        Fτd = []
    ) where {T,Nn,Nb,Ne,Ni}

        return LinearMechanism(mechanism, xd, vd, qd, ωd, Fτd, eqcids)
    end
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, mechanism::AbstractMechanism{T,Nn,Nb,Ne,0}) where {T,Nn,Nb,Ne}
    summary(io, mechanism)
    println(io, " with ", Nb, " bodies and ", Ne, " constraints")
    println(io, " Δt: "*string(mechanism.Δt))
    println(io, " g:  "*string(mechanism.g))
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, mechanism::AbstractMechanism{T,Nn,Nb,Ne,Ni}) where {T,Nn,Nb,Ne,Ni}
    summary(io, mechanism)
    println(io, " with ", Nb, " bodies, ", Ne, " equality constraints, and ", Ni, " inequality constraints")
    println(io, " Δt: "*string(mechanism.Δt))
    println(io, " g:  "*string(mechanism.g))
end


text:
# Define delay differential equation for change solute concentrations with
# respect to stream length, assuming steady state with respect to time.
macro stream_model(reactionKinetics, shape)
    if shape == "powerLaw"
        if reactionKinetics == "firstOrder"
        return :(
            function stream_model!(du,u,p,t)
                τ₋, τ₊, α, k, qPerUnitHZ, sHZ, sCH, v = p
                integral = u[1] * quadgk(τ -> E(τ, α, τ₋, τ₊)*firstOrder(τ, k),τ₋, τ₊)[1]
                du[1] = (qPerUnitHZ*(sHZ/sCH)*(integral - u[1]))/v
            end
                )
        elseif reactionKinetics == "firstOrderProductionAndConsumption"
        return :(
            function stream_steady_state_withNitrification!(du,u,p,t)
                τ₋, τ₊, α, qPerUnitHZ, sHZ, sCH, v, kprod, kcons, substrate0 = p
                integral =  quadgk(τ -> E(τ, α, τ₋, τ₊)*firstOrderProductionAndConsumption(τ, kprod, kcons, substrate0, u[1]),τ₋, τ₊)[1]
                du[1] = (qPerUnitHZ*(sHZ/sCH)*(integral - u[1]))/v
            end
                )
        end
    end
  end

max_stars_repo_name: elizabethjmohr/hyporheicBGC

max_issues_repo_name: elizabethjmohr/hyporheicBGC

max_forks_repo_name: elizabethjmohr/hyporheicBGC


text:
module SOMone

export SOM,maxentropy

using Statistics, Printf, Random, StatsBase

include("metrics.jl")
include("neighbourhood.jl")
include("initialisation.jl")
include("selforganise.jl")

end


text:
## =============================================================================
## To be deprecated for consolidation, can revisit after all is one function #675
## =============================================================================

# UP fetch + condition

function getMsgUpChannel(cdat::BayesTreeNodeData)
  Base.depwarn("Deprecated - using Edge Messages for consolidation, #675 - 2", :getMsgUpChannel)
  return cdat.upMsgChannel
end
getMsgUpChannel(cliq::TreeClique) = getMsgUpChannel(getCliqueData(cliq))

function putCliqueMsgUp!(cdat::BayesTreeNodeData, upmsg::LikelihoodMessage)
  Base.depwarn("Deprecated - using take! model", :putCliqueMsgUp!)
  # new replace put! interface
  cdc_ = getMsgUpChannel(cdat)
  if isready(cdc_)
    # first clear an existing value
    take!(cdc_)
  end
  put!(cdc_, upmsg)
  # cdat.upMsg = msg
end

## =============================================================================
## To be deprecated fetch
## =============================================================================

"""
    $SIGNATURES

Notify of new up status and message.

DevNotes
- Major part of #459 consolidation effort.
- FIXME require consolidation
  - Perhaps deprecate putMsgUpInitStatus! as separate function?
"""
function prepPutCliqueStatusMsgUp!( csmc::CliqStateMachineContainer,
                                    status::Symbol=getCliqueStatus(csmc.cliq);
                                    dfg::AbstractDFG=csmc.cliqSubFg,
                                    upmsg=prepCliqueMsgUpConsolidated(dfg, csmc.cliq, status, logger=csmc.logger)  )

  Base.depwarn("Deprecated - using take! model", :prepPutCliqueStatusMsgUp!)
  #
  # TODO replace with msg channels only

  # put the init upmsg
  cd = getCliqueData(csmc.cliq)

  setCliqueStatus!(csmc.cliq, status)

  # NOTE consolidate with upMsgChannel #459
  putCliqueMsgUp!(cd, upmsg)
    # TODO remove as part of putCliqueMsgUp!
    # new replace put! interface
    # cdc_ = getMsgUpChannel(cd)
    # if isready(cdc_)
    #   # first clear an existing value
    #   take!(cdc_)
    # end
    # put!(cdc_, upmsg)

  notify(getSolveCondition(csmc.cliq))
  # took ~40 hours to figure out that a double norification fixes the problem with hex init
  sleep(0.1)
  notify(getSolveCondition(csmc.cliq))

  # also notify parent as part of upinit (PUSH NOTIFY PARENT), received at `slowWhileInit_StateMachine`
  prnt = getParent(csmc.tree, csmc.cliq)
  0 < length(prnt) ? notify(getSolveCondition(prnt[1])) : nothing

  infocsm(csmc, "prepPutCliqueStatusMsgUp! -- notified status=$(upmsg.status) with msg keys $(collect(keys(upmsg.belief)))")

  # return new up messages in case the user wants to see
  return upmsg
end

# DOWN fetch + condition
function getDwnMsgConsolidated(btnd::BayesTreeNodeData)
  Base.depwarn("Deprecated - using Edge Messages for consolidation, #675 - 2", :getDwnMsgConsolidated)
  btnd.dwnMsgChannel
end
getDwnMsgConsolidated(cliq::TreeClique) = getDwnMsgConsolidated(getCliqueData(cliq))

function fetchDwnMsgConsolidated(btnd::BayesTreeNodeData)
  Base.depwarn("Deprecated - using take! model", :fetchDwnMsgConsolidated)
  fetch(getDwnMsgConsolidated(btnd))
end
fetchDwnMsgConsolidated(cliq::TreeClique) = fetchDwnMsgConsolidated(getCliqueData(cliq))

function putDwnMsgConsolidated!(btnd::BayesTreeNodeData, msg::LikelihoodMessage)
  Base.depwarn("Deprecated - using take! model", :putDwnMsgConsolidated!)
  # need to get the current solvableDims
  dmc = getDwnMsgConsolidated(btnd)
  if isready(dmc)
    take!(dmc)
  end
  put!(dmc, msg)
end
putDwnMsgConsolidated!(cliq::TreeClique, msg::LikelihoodMessage) = putDwnMsgConsolidated!(getCliqueData(cliq), msg)

"""
    $SIGNATURES

Consolidated downward messages generator and Channel sender.

Notes
- Post #459
"""
function prepPutCliqueStatusMsgDwn!(csmc::CliqStateMachineContainer,
                                    status::Symbol=getCliqueStatus(csmc.cliq);
                                    dfg::AbstractDFG=csmc.cliqSubFg,
                                    childSolvDims::Dict{Int,Float64} = Dict{Int,Float64}(),
                                    dwnmsg=prepSetCliqueMsgDownConsolidated!(dfg, csmc.cliq, LikelihoodMessage(status=status, childSolvDims=childSolvDims), csmc.logger, status=status )  )
  #
  Base.depwarn("Deprecated - using take! model", :prepPutCliqueStatusMsgDwn!)

  cd = getCliqueData(csmc.cliq)

  setCliqueStatus!(csmc.cliq, status)

  # NOTE consolidate with upMsgChannel #459
  putDwnMsgConsolidated!(cd, dwnmsg)

  notify(getSolveCondition(csmc.cliq))
  # double notification fixes the problem with hex init (likely due to old lock or something, remove with care)
  sleep(0.1)
  notify(getSolveCondition(csmc.cliq))

  infocsm(csmc, "prepPutCliqueStatusMsgDwn! -- notified status=$(dwnmsg.status), msgs $(collect(keys(dwnmsg.belief))), childSolvDims=$childSolvDims")

  status
end

"""
    $(SIGNATURES)

Return the last up message stored in This `cliq` of the Bayes (Junction) tree.
"""
fetchMsgUpThis(cdat::BayesTreeNodeData) = fetch(getMsgUpChannel(cdat))  # cdat.upMsg    # TODO rename to fetchMsgUp
fetchMsgUpThis(cliql::TreeClique) = fetchMsgUpThis(getCliqueData(cliql))
fetchMsgUpThis(btl::AbstractBayesTree, frontal::Symbol) = fetchMsgUpThis(getClique(btl, frontal))

function fetchMsgsUpChildrenDict( treel::AbstractBayesTree,
                                  cliq::TreeClique )
  #
  msgs = Dict{Int, LikelihoodMessage}()
  for chld in getChildren(treel, cliq)
    msgs[chld.index] = fetchMsgUpThis(chld)
  end

  return msgs
end
fetchMsgsUpChildrenDict( csmc::CliqStateMachineContainer ) = fetchMsgsUpChildrenDict( csmc.tree, csmc.cliq )

"""
    $SIGNATURES

Get and return upward belief messages as stored in child cliques from `treel::AbstractBayesTree`.

Notes
- Use last parameter to select the return format.
- Pull model #674

DevNotes
- Consolidate fetchChildrenStatusUp, getMsgsUpInitChildren
- FIXME update refactor to fetch or take, #855

Related

fetchMsgsUpChildrenDict
"""
function fetchMsgsUpChildren( treel::AbstractBayesTree,
                              cliq::TreeClique,
                              ::Type{TreeBelief} )
  #
  chld = getChildren(treel, cliq)
  retmsgs = Vector{LikelihoodMessage}(undef, length(chld))
  for i in 1:length(chld)
    retmsgs[i] = fetchMsgUpThis(chld[i])
  end
  return retmsgs
end

function fetchMsgsUpChildren( csmc::CliqStateMachineContainer,
                              ::Type{TreeBelief}=TreeBelief )
  #
  # TODO, replace with single channel stored in csmcs or cliques
  fetchMsgsUpChildren(csmc.tree, csmc.cliq, TreeBelief)
end

## TODO Consolidate/Deprecate below

"""
    $SIGNATURES

Fetch (block) caller until child cliques of `cliq::TreeClique` have valid csm status.

Notes:
- Return `::Dict{Symbol}` indicating whether next action that should be taken for each child clique.
- See status options at `getCliqueStatus(..)`.
- Can be called multiple times
"""
function fetchChildrenStatusUp( tree::AbstractBayesTree,
                                cliq::TreeClique,
                                logger=ConsoleLogger() )
  #
  ret = Dict{Int, Symbol}()
  chlr = getChildren(tree, cliq)
  for ch in chlr
      # # FIXME, why are there two steps getting cliq status????
      # chst = getCliqueStatus(ch)  # TODO, remove this
    with_logger(logger) do
      @info "cliq $(cliq.index), child $(ch.index) isready(initUpCh)=$(isready(getMsgUpChannel(ch)))."
    end
    flush(logger.stream)
    # either wait to fetch new result, or report or result
    ret[ch.index] = (fetch(getMsgUpChannel(ch))).status
    with_logger(logger) do
      @info "ret[$(ch.index)]=$(ret[ch.index])."
    end
  end

  return ret
end

# FIXME TEMPORARY CONSOLIDATION FUNCTIONS
# this method adds children and own up msg info to the returning Dict.
# own information is added to capture information from cousins during down init.
function getMsgsUpInitChildren( treel::AbstractBayesTree,
                                cliq::TreeClique,
                                ::Type{TreeBelief};
                                skip::Vector{Int}=Int[])
  #
  chld = getChildren(treel, cliq)
  retmsgs = Dict{Int, LikelihoodMessage}()
  # add possible information that may have come via grandparents from elsewhere in the tree
  if !(cliq.index in skip)
    thismsg = fetchMsgUpThis(cliq)
    retmsgs[cliq.index] = thismsg
  end

  # now add information from each of the child cliques (no longer all stored in prnt i.e. old push #674)
  for ch in chld
    chmsg = fetchMsgUpThis(ch)
    if !(ch.index in skip)
      retmsgs[ch.index] = chmsg
    end
  end
  return retmsgs
end

function getMsgsUpInitChildren( csmc::CliqStateMachineContainer,
                                ::Type{TreeBelief}=TreeBelief;
                                skip::Vector{Int}=Int[] )
  #
  # TODO, replace with single channel stored in csmcs or cliques
  getMsgsUpInitChildren(csmc.tree, csmc.cliq, TreeBelief, skip=skip)
end


text:
export RTMHandler,
       rtm_connect,
       send_event,
       RTMClient,
       on_reply,
       on_event,
       on_error,
       on_disconnect,
       on_connect,
       attach

using DandelionWebSockets
import DandelionWebSockets: on_text, on_binary,
                            state_connecting, state_open, state_closing, state_closed,
                            AbstractWSClient
import JSON

#
# RTMHandler defines an interface for handling RTM events.
#

abstract RTMHandler

on_reply(h::RTMHandler, id::Int, event::Event) =
    error("on_reply not implemented for $(h) and/or $(event)")

on_event(h::RTMHandler, event::Event) =
    error("on_event not implemented for $(h) and/or $(event)")

# on error is called when there is a problem with receiving the message, such as invalid JSON.
# This is not an error sent by Slack, but an error caught in this code.
on_error(h::RTMHandler, e::EventError) =
    error("on_error not implemented for $(h) or $e")

on_disconnect(h::RTMHandler) = error("on_disconnect not implemented for $(h)")
on_connect(h::RTMHandler) = error("on_connect not implemented for $(h)")

#
#
#
immutable UnsetRTMHandler <: RTMHandler end

on_reply(h::UnsetRTMHandler, id::Int, event::Event) =
    error("No RTM handler set! Call attach(::RTMClient, ::RTMHandler) before connect")

on_event(h::UnsetRTMHandler, event::Event) =
    error("No RTM handler set! Call attach(::RTMClient, ::RTMHandler) before connect")

# on error is called when there is a problem with receiving the message, such as invalid JSON.
# This is not an error sent by Slack, but an error caught in this code.
on_error(h::UnsetRTMHandler, e::EventError) =
    error("No RTM handler set! Call attach(::RTMClient, ::RTMHandler) before connect")

on_disconnect(h::UnsetRTMHandler) =
    error("No RTM handler set! Call attach(::RTMClient, ::RTMHandler) before connect")

on_connect(h::UnsetRTMHandler) =
    error("No RTM handler set! Call attach(::RTMClient, ::RTMHandler) before connect")

#
# RTMWebSocketHandler takes events from a WebSocket connection and converts to RTM events.
#

type RTMWebSocket <: WebSocketHandler
    handler::RTMHandler
    connection_retry::AbstractRetry

    RTMWebSocket(retry::AbstractRetry) = new(UnsetRTMHandler(), retry)
end

function on_text(rtm::RTMWebSocket, text::UTF8String)
    dict::Dict{AbstractString, Any} = Dict()
    try
        dict = JSON.parse(text)
    catch ex
        if isa(ex, ErrorException)
            on_error(rtm.handler, InvalidJSONError(text))
            return
        end
        rethrow(ex)
    end

    if !haskey(dict, "type")
        # Special case: A message ack does not have a type property, but is instead identified by
        # having a text and a reply_to property.
        # We fake a type "message_ack" for it, so it's handled by the same code as everything else.
        if haskey(dict, "text") && haskey(dict, "reply_to")
            dict["type"] = utf8("message_ack")
        else
            on_error(rtm.handler, MissingTypeError(text))
            return
        end
    end

    event_type = find_event(dict["type"])
    if event_type == nothing
        on_error(rtm.handler, UnknownEventTypeError(text, dict["type"]))
        return
    end

    event = nothing
    try
        event = deserialize(event_type, dict)
    catch ex
        if isa(ex, KeyError)
            on_error(rtm.handler, DeserializationError(ex.key, text, event_type))
            return
        else
            rethrow(ex)
        end
    end

    if haskey(dict, "reply_to")
        message_id = dict["reply_to"]
        on_reply(rtm.handler, message_id, event)
    else
        on_event(rtm.handler, event)
    end
end

function state_open(rtm::RTMWebSocket)
    on_connect(rtm.handler)
    reset(rtm.connection_retry)
end

function state_closed(rtm::RTMWebSocket)
    on_disconnect(rtm.handler)
    retry(rtm.connection_retry)
end

# Implement a warning here, as Slack shouldn't send binary messages.
on_binary(::RTMWebSocket, ::Vector{UInt8}) = nothing
state_connecting(::RTMWebSocket) = nothing
state_closing(t::RTMWebSocket) = nothing

attach(t::RTMWebSocket, handler::RTMHandler) = t.handler = handler

#
# RTMClient is an object for sending events to Slack.
#

abstract AbstractRTMClient

default_backoff = RandomizedBackoff(Backoff(5.0, 200.0), MersenneTwister(), 3.0)

throttled_client_factory = handler -> ThrottledWSClient(WSClient(), 1.0)

type RTMClient <: AbstractRTMClient
    ws_client::AbstractWSClient
    next_id::Int64
    rtm_ws::RTMWebSocket
    token::Token

    function RTMClient(token::Token;
                       connection_retry::AbstractRetry=Retry(default_backoff, x -> nothing),
                       ws_client_factory=throttled_client_factory)
        rtm_ws = RTMWebSocket(connection_retry)
        ws_client = ws_client_factory(rtm_ws)
        c = new(ws_client, 1, rtm_ws, token)
        set_function(connection_retry, () -> rtm_connect(c))
        c
    end
end

show(io::IO, c::RTMClient) = show(io, "RTMClient($(c.next_id), $(c.token), $(c.rtm_ws))")

function send_event(c::RTMClient, event::OutgoingEvent)
    this_id = c.next_id
    c.next_id += 1

    dict = serialize(event)
    dict["id"] = this_id
    text = utf8(JSON.json(dict))
    send_text(c.ws_client, text)

    this_id
end

close(c::RTMClient) = stop(c.ws_client)
attach(c::RTMClient, handler::RTMHandler) = attach(c.rtm_ws, handler)

function rtm_connect(client::RTMClient;
                     requests=real_requests)
    try
        status, response = makerequest(
            RtmStart(client.token, Nullable(), Nullable(), Nullable()), requests)
        wsconnect(client.ws_client, Requests.URI(response.url), client.rtm_ws)
    catch ex
        println("Exception when connecting: $ex")
        state_closed(client.rtm_ws)
    end
end


text:
using Plots
using DifferentialEquations
using LaTeXStrings
using ProgressMeter
plotlyjs()

# Rossler space equation
g = @ode_def Rossler begin
    dx = -y - z
    dy = x + a*y
    dz = b + z*(x-c)
end a=>0.2 b=>0.2 c=>5.7

# Simulate
time_span = (0.00, 3000.0)
x0 = [1.0, 1.0, 2.0]
prob = ODEProblem(g, x0, time_span)
sol = solve(prob, RK4(), dt=0.1)

plot(sol, vars=(1, 2))

# Plot Poincare section
xconst = 0
plot!(zeros(100), linspace(0.0, 10.0, 100), lw=4, color="black")

# Get Poincare section
init = 600
poincare_y = []
poincare_z = []
p = zeros(3)
for i = init:length(sol)-1
    if sol[i][1] > 0 && sol[i+1][1] < 0
        p = (sol[i+1] - sol[i])/(sol[i+1][1] - sol[i][1])*(xconst - sol[i][1]) + sol[i]
        push!(poincare_y, p[2])
        push!(poincare_z, p[3])
    end
end

scatter(poincare_y, poincare_z)

# Bifurcation diagram
list_y = []
list_params = []
@showprogress 1 "Computing..." for param = 0.01:0.01:2.00
    # Solve Rossler
    h = Rossler(b=param)
    prob = ODEProblem(h, x0, time_span)
    sol = solve(prob, RK4(), dt=0.1)
    # Find variables
    for i = init:length(sol)-1
        if sol[i][1] > 0 && sol[i+1][1] < 0
            p = (sol[i+1] - sol[i])/(sol[i+1][1] - sol[i][1])*(xconst - sol[i][1]) + sol[i]
            push!(list_y, p[2])
            push!(list_params, param)
        end
    end
end

pyplot()
scatter(list_params, list_y, color="black", markersize=0.01, markerstrokewidth=0.5)

max_stars_repo_name: MACSIN-UFMG/julia-rossler-example

max_issues_repo_name: MACSIN-UFMG/julia-rossler-example

max_forks_repo_name: MACSIN-UFMG/julia-rossler-example