text:
#=

=#

################################################################################
#
# Create tiny layers and tiny network
#
################################################################################

########################################
# Input Layer
########################################

# Slices of Arrays
tinyslice{R<:Real}(X::Array{R,1}, ld_out, hd_dim, num_channels, num_obs) = X[1:num_obs]
tinyslice{R<:Real}(X::Array{R,2}, ld_out, hd_dim, num_channels, num_obs) = X[1:ld_out, 1:num_obs]
tinyslice{R<:Real}(X::Array{R,3}, ld_out, hd_dim, num_channels, num_obs) = X[1:hd_dim, 1:hd_dim, 1:num_obs]
tinyslice{R<:Real}(X::Array{R,4}, ld_out, hd_dim, num_channels, num_obs) = X[1:hd_dim, 1:hd_dim, 1:num_channels, 1:num_obs]

# InputLayer
function tinylayer(ℓ::InputLayer, ld_out, hd_dim,
                   num_channels, batch_size)
    X = tinyslice(ℓ.A, ld_out, hd_dim, num_channels, batch_size)
    rows = prod(size(X)[1:end-1])
    cols = size(X)[end]
    InputLayer(
    X, reshape(X, rows, cols)
    )
end

########################################
# Hidden Layers
########################################

# FullyConnectedLayer
function tinylayer(ℓ::FullyConnectedLayer, ld_in, ld_out, hd_dim, num_channels,
                                         patch_dim, num_patches, pooling_stride,
                                         num_outputs, num_classes, batch_size)

    # Datatype
    F = eltype(ℓ.A)

    # Take small slices out of ℓ
    b = ℓ.b[1:ld_out]
    W = ℓ.W[1:ld_out, 1:ld_in]

    # Create layer
    FullyConnectedLayer(
    b, W,
    zeros(F, ld_out, batch_size),
    zeros(F, ld_out, batch_size),
    zeros(F, ld_out, batch_size),
    zeros(b), zeros(W), zeros(b), zeros(W),
    ℓ.f!, ℓ.∇f!
    )
end

# WordEmbeddingLayer
function tinylayer(ℓ::WordEmbeddingLayer, ld_in, ld_out, hd_dim, num_channels,
                                         patch_dim, num_patches, pooling_stride,
                                         num_outputs, num_classes, batch_size)

    # Datatype
    F = eltype(ℓ.A)

    # Small constants
    vocab_size = 20
    num_words = 2
    num_feats = ld_out / num_words

    # Take small slices out of ℓ
    W = ℓ.W[1:ld_out, 1:ld_in] # TODO

    # Create layer
    _, W = init_uniform(F, num_feats, vocab_size)
    A = zeros(F, num_words * num_feats, batch_size)
    Δ = zeros(F, num_words * num_feats, batch_size)
    WordEmbeddingLayer(
    W,
    A, Δ,
    zeros(W), zeros(W)
    )
end

# SparseEncoderLayer
function tinylayer(ℓ::SparseEncoderLayer, ld_in, ld_out, hd_dim, num_channels,
                                         patch_dim, num_patches, pooling_stride,
                                         num_outputs, num_classes, batch_size)

    # Datatype
    F = eltype(ℓ.A)

    # Take small slices out of ℓ
    b = ℓ.b[1:ld_out]
    W = ℓ.W[1:ld_out, 1:ld_in]

    # Create layer
    SparseEncoderLayer(
    b, W,
    zeros(F, ld_out, batch_size),
    zeros(F, ld_out, batch_size),
    zeros(F, ld_out, batch_size),
    zeros(b), zeros(W), zeros(b), zeros(W),
    zeros(b), zeros(b), ℓ.p, ℓ.β,
    ℓ.f!, ℓ.∇f!, ℓ.δ!, F(batch_size)
    )
end

# ConvolutionLayer
function tinylayer(ℓ::ConvolutionLayer, ld_in, ld_out, hd_dim, num_channels,
                                        patch_dim, num_patches, pooling_stride,
                                        num_outputs, num_classes, batch_size)

    # Datatype
    F = eltype(ℓ.A)

    # Patch slice indices
    num_patch_dims = length(size(ℓ.W))
    if num_patch_dims == 3
        patch_ids = (1:patch_dim, 1:patch_dim, 1:num_patches)
    elseif num_patch_dims == 4
        patch_ids = (1:patch_dim, 1:patch_dim, 1:num_channels, 1:num_patches)
    end

    # Take small slices out of b and W
    b = ℓ.b[patch_ids[end]]
    W = ℓ.W[patch_ids...]

    # Convolution dims
    num_conv_dims = length(size(ℓ.A))
    if num_conv_dims == 4
        conv_dims = (hd_dim, hd_dim, num_patches, batch_size)
    elseif num_conv_dims == 5
        conv_dims = (hd_dim, hd_dim, num_channels, num_patches, batch_size)
    end

    # Create Z, A and Δ
    Z = zeros(F, conv_dims)
    A = zeros(F, conv_dims)
    Δ = zeros(F, conv_dims)

    # For reshaping
    rows = prod(conv_dims[1:end-1])

    # Create layer
    ConvolutionLayer(
    b, W,
    Z, reshape(Z, rows, batch_size),
    A, reshape(A, rows, batch_size),
    Δ, reshape(Δ, rows, batch_size),
    zeros(b), zeros(W), zeros(b), zeros(W),
    ℓ.f!, ℓ.∇f!
    )
end

# MeanPoolLayer
function tinylayer(ℓ::MeanPoolLayer, ld_in, ld_out, hd_dim, num_channels,
                                     patch_dim, num_patches, pooling_stride,
                                     num_outputs, num_classes, batch_size)

    # Datatype
    F = eltype(ℓ.A)

    # Pool dims
    num_dims = length(size(ℓ.A))
    # TODO - Can there be 5 dims??
    if num_dims == 4
        dims = (hd_dim, hd_dim, num_patches, batch_size)
    elseif num_dims == 5
        dims = (hd_dim, hd_dim, num_channels, num_patches, batch_size)
    end

    # Create empty A and Δ
    A = zeros(dims)
    Δ = zeros(dims)

    # For reshaping
    rows = prod(dims[1:end-1])

    # Create layer
    MeanPoolLayer(
    A, reshape(A, rows, batch_size),
    Δ, reshape(Δ, rows, batch_size)
    )
end

# MaxPoolLayer
function tinylayer(ℓ::MaxPoolLayer, ld_in, ld_out, hd_dim, num_channels,
                                    patch_dim, num_patches, pooling_stride,
                                    num_outputs, num_classes, batch_size)

    # Datatype
    F = eltype(ℓ.A)

    # Pool dims
    num_dims = length(size(ℓ.A))
    prev_dim = hd_dim * pooling_stride
    # TODO - Can there be 5 dims??
    if num_dims == 4
        dims = (hd_dim, hd_dim, num_patches, batch_size)
        dims_in = (prev_dim, prev_dim, num_patches, batch_size)
    elseif num_dims == 5
        dims = (hd_dim, hd_dim, num_channels, num_patches, batch_size)
        dims_in = (prev_dim, prev_dim, num_channels, num_patches, batch_size)
    end

    # Create empty A and Δ
    A = zeros(dims)
    Δ = zeros(dims)

    # For reshaping
    rows = prod(dims[1:end-1])

    # Create layer
    MaxPoolLayer(
    A, reshape(A, rows, batch_size),
    Δ, reshape(Δ, rows, batch_size),
    zeros(Int, dims_in...)
    )
end

########################################
# Output Layers
########################################

# LinearOutputLayer
function tinylayer(ℓ::LinearOutputLayer, ld_in, ld_out, hd_dim, num_channels,
                                         patch_dim, num_patches, pooling_stride,
                                         num_outputs, num_classes, batch_size)

    # Datatypes
    F = eltype(ℓ.A)

    # Take small slices out of ℓ weights and target
    b = ℓ.b[:]
    W = ℓ.W[:, 1:ld_in]
    y = ℓ.y[1:batch_size]

    # Create layer
    LinearOutputLayer(
    b, W,
    zeros(F, 1, batch_size), zeros(F, 1, batch_size),
    zeros(b), zeros(W), zeros(b), zeros(W),
    ℓ.δ!, y, ℓ.cost, F(batch_size)
    )
end

# MultiLinearOutputLayer
function tinylayer(ℓ::MultiLinearOutputLayer, ld_in, ld_out, hd_dim, num_channels,
                                         patch_dim, num_patches, pooling_stride,
                                         num_outputs, num_classes, batch_size)

    # Datatype
    F = eltype(ℓ.A)

    # Take small slices out of ℓ weights and target
    b = ℓ.b[1:num_outputs]
    W = ℓ.W[1:num_outputs, 1:ld_in]
    y = ℓ.y[1:num_outputs, 1:batch_size]

    # Create layer
    MultiLinearOutputLayer(
    b, W,
    zeros(F, num_outputs, batch_size), zeros(F, num_outputs, batch_size),
    zeros(b), zeros(W), zeros(b), zeros(W),
    ℓ.δ!, y, ℓ.cost, F(batch_size)
    )
end

# LogisticOutputLayer
function tinylayer(ℓ::LogisticOutputLayer, ld_in, ld_out, hd_dim, num_channels,
                                         patch_dim, num_patches, pooling_stride,
                                         num_outputs, num_classes, batch_size)

    # Datatype
    F = eltype(ℓ.A)

    # Take small slices out of ℓ weights and target
    b = ℓ.b[:]
    W = ℓ.W[:, 1:ld_in]
    y = rand(0:1, batch_size)   #ℓ.y[1:batch_size]

    # Create layer
    LogisticOutputLayer(
    b, W,
    zeros(F, 1, batch_size), zeros(F, 1, batch_size), zeros(F, 1, batch_size),
    zeros(b), zeros(W), zeros(b), zeros(W),
    ℓ.f!, ℓ.δ!, y, ℓ.cost, F(batch_size)
    )
end

# SoftmaxOutputLayer
function tinylayer(ℓ::SoftmaxOutputLayer, ld_in, ld_out, hd_dim, num_channels,
                                         patch_dim, num_patches, pooling_stride,
                                         num_outputs, num_classes, batch_size)

    # Datatype
    F = eltype(ℓ.A)

    # Take small slices out of ℓ weights and target
    b = ℓ.b[1:num_classes]
    W = ℓ.W[1:num_classes, 1:ld_in]
    y = rand(1:num_classes, batch_size)   #ℓ.y[1:batch_size]

    # Create layer
    SoftmaxOutputLayer(
    b, W,
    zeros(F, num_classes, batch_size), zeros(F, 1, batch_size),
    zeros(F, num_classes, batch_size), zeros(F, num_classes, batch_size),
    zeros(b), zeros(W), zeros(b), zeros(W),
    ℓ.f!, ℓ.δ!, y, ℓ.cost, F(batch_size)
    )
end

########################################
# Network
########################################

# Network
function tinynet(net::NeuralNet)

    # Start values that gett changed as layers are built
    ld_in = 9
    ld_out = 9
    hd_dim = 3

    # Constant values - not changed
    num_channels = 3
    patch_dim = 3
    num_patches = 5
    pooling_stride = 2
    num_outputs = min(9, size(net.data.X_train, 1))
    num_classes = min(5, length(unique(net.data.y_train)))
    batch_size = min(40, net.batch_size)

    # Create tiny hidden layers from last to first
    tinylayers = NetLayer[]
    for (ℓ, prev_ℓ) in zip(net.layers[end:-1:2], net.layers[(end-1):-1:1])

        # Update ld_in var if chaining previous high dim with this low dim layer
        if typeof(ℓ) <: LDL && typeof(prev_ℓ) <: HDL
            prev_ℓ_dims = length(size(prev_ℓ.A))
            if prev_ℓ_dims == 4
                ld_in = hd_dim * hd_dim * num_patches
            elseif prev_ℓ_dims == 5
                ld_in = hd_dim * hd_dim * num_channels * num_patches
            end
        end

        # Create tiny layer and add to vector
        push!(tinylayers, tinylayer(ℓ, ld_in, ld_out, hd_dim, num_channels,
                                    patch_dim, num_patches, pooling_stride,
                                    num_outputs, num_classes, batch_size))

        # Update high dimension var for use in previous high dimension layer
        if typeof(ℓ) <: PoolLayer
            hd_dim = hd_dim * pooling_stride
        elseif typeof(ℓ) <: ConvolutionLayer
            hd_dim = hd_dim + patch_dim - 1
        end
    end

    # Create tiny input layer
    ℓ = net.layers[1]
    push!(tinylayers, tinylayer(ℓ, ld_out, hd_dim, num_channels, batch_size))

    # Reverse - i.e. flip to ordering from input o output
    tinylayers = reverse(tinylayers)

    # Check that the layer dimensions allow chaining
    checkdims(tinylayers)

    # Create tiny data container
    tdata = Data(tinylayers[1].A, tinylayers[end].y)

    # Create tiny network
    F = net.dtype
    NeuralNet(F, tdata, tinylayers, net.reg_cost, net.reg_grad!, net.λ,
              batch_size, Array(F, 0), Array(F, 0), Array(F, 0))
end

################################################################################
#
# Display functions
#
################################################################################

# Display detail - show difference for every parameter
function display_grad_check(numeric_grad, grad)
    # Display the difference between numeric and analytic gradient for each parameter
    println()
    @printf "  Numeric Gradient     Analytic Gradient         Difference\n"
    for i=1:length(numeric_grad)
        @printf "%18e %21e %18e \n" numeric_grad[i] grad[i] numeric_grad[i] - grad[i]
    end
    println()
end

# Display summary - show relative difference
function display_grad_reldiff(numeric_grad, grad)
    # Display the norm of the difference between two solutions.
    # The diff should be less than 1e-7
    diff = norm(view(grad, :) - view(numeric_grad, :))
    denom = max(norm(view(grad, :)), norm(view(numeric_grad, :)))
    rel_diff = diff == 0.0 ? 0.0 : diff / denom
    println("  $rel_diff")
    println()
end

################################################################################
#
# Check the gradients
#
################################################################################

function check_gradients(bignet::NeuralNet; detail=false)

    # Craete tinynet
    net = tinynet(bignet)

    # Calculate gradients analytically
    fwd_prop!(net)
    bwd_prop!(net)

    # Size of perturbation
    ϵ = 1e-4

    # Output Layer
    L = net.layers[end]

    # Print
    println("Check gradients of a tiny network with similar architecture:\n")
    ℓ = net.layers[1]
    str = "Layer 1, " * string(typeof(ℓ))[7:end] * ", Dimensions - $(size(ℓ.A)) :\n\n"
    print_with_color(:blue, str)

    for (n, ℓ) in enumerate(net.layers[2:end])

        # Print results
        print_with_color(:blue, "Layer $(n+1), " * displaystr(ℓ) * " :\n\n")

        # Skip pooling layers
        if typeof(ℓ) <: PoolLayer
            println("  No parameters in pooling layer.\n")
            continue
        end

        # Initialise arrays
        numeric_grad_b = zeros(ℓ.b)
        numeric_grad_W = zeros(ℓ.W)

        # Estimate gradient for each bias parameter separately
        for i in eachindex(ℓ.b)

            # Store value of parameter
            tmp = ℓ.b[i]

            # Perturb the parameter and calculate the cost
            ℓ.b[i] += ϵ
            fwd_prop!(net)
            cost1 = L.cost(L.A, L.y) / L.m +
                    sparse_cost(net)

            ℓ.b[i] -= 2ϵ
            fwd_prop!(net)
            cost2 = L.cost(L.A, L.y) / L.m +
                    sparse_cost(net)

            # Calculate numeric gradient using central difference
            numeric_grad_b[i] = (cost1 - cost2) ./ (2ϵ)

            # Reset the parameter
            ℓ.b[i] = tmp
        end

        # Estimate gradient for each weight parameter separately
        for i in eachindex(ℓ.W)

            # Store value of parameter
            tmp = ℓ.W[i]

            # Perturb the parameter and calculate the cost
            ℓ.W[i] += ϵ
            fwd_prop!(net)
            cost1 = L.cost(L.A, L.y) / L.m +
                    net.reg_cost(net) +
                    sparse_cost(net)

            ℓ.W[i] -= 2ϵ
            fwd_prop!(net)
            cost2 = L.cost(L.A, L.y) / L.m +
                    net.reg_cost(net) +
                    sparse_cost(net)

            # Calculate numeric gradient using central difference
            numeric_grad_W[i] = (cost1 - cost2) ./ (2ϵ)

            # Reset the parameter
            ℓ.W[i] = tmp
        end

        # Detail
        if detail
            println("  Bias parameters :")
            display_grad_check(numeric_grad_b, ℓ.∇b)
            println("  Weight parameters :")
            display_grad_check(numeric_grad_W, ℓ.∇W)
        end

        # Summary
        println("  Bias parameters relative difference :")
        display_grad_reldiff(numeric_grad_b, ℓ.∇b)
        println("  Weight parameters relative difference :")
        display_grad_reldiff(numeric_grad_W, ℓ.∇W)

    end
end


text:
@node Categorical Stochastic [ out, p ]

@average_energy Categorical (q_out::Categorical, q_p::Any) = -sum(probvec(q_out) .* logmean(q_p))

ext: jl


text:
"Computes the joint state and control dimensions of the combined system."
function xu_dims(subsystems::NTuple)
    nx, nu = 0, 0;
    xids, uids = [], []
    for sub in subsystems
        nxᵢ = n_states(sub); nuᵢ = n_controls(sub)
        push!(xids, SVector{nxᵢ}((nx+1):(nx+nxᵢ)))
        push!(uids, SVector{nuᵢ}((nu+1):(nu+nuᵢ)))
        nx += nxᵢ; nu += nuᵢ
    end
    xyids = (xid[xyindex(sub)] for (xid, sub) in zip(xids, subsystems))
    return nx, nu, Tuple(xids), Tuple(uids), Tuple(xyids)
end

function ξ_dims(subsystems::NTuple)
    nξ = 0
    ξids = []
    for sub in subsystems
        nξᵢ = n_linstates(sub)
        push!(ξids, SVector{nξᵢ}((nξ+1):(nξ+nξᵢ)))
        nξ += nξᵢ
    end
    ξxyids = (ξid[ξxyindex(sub)] for (ξid, sub) in zip(ξids, subsystems))
    return nξ, Tuple(ξids), Tuple(ξxyids)
end

Base.@pure function common_linstyle(subsystems::NTuple)
    linstyles = unique(LinearizationStyle.(subsystems))
    return length(linstyles) == 1 ? linstyles[1] : return JacobianLinearization()
end

struct ProductSystem{ΔT,nx,nu,np,nξ,TS<:NTuple{np,<:ControlSystem{ΔT}},TXI<:NTuple{np},
                     TUI<:NTuple{np},TXY<:NTuple{np},TXII<:NTuple{np},TXIXY<:NTuple{np}}<:ControlSystem{ΔT,nx,nu}
    subsystems::TS
    xids::TXI
    uids::TUI
    xyids::TXY
    ξids::TXII
    ξxyids::TXIXY
    function ProductSystem(subsystems::TS) where {ΔT,np,TS<:NTuple{np,<:ControlSystem{ΔT}}}
        nx, nu, xids, uids, xyids = xu_dims(subsystems)
        # Notice: this very local type instability is acceptable because the
        # ProductSystem will only be construted once
        if common_linstyle(subsystems) isa FeedbackLinearization
            nξ, ξids, ξxyids = ξ_dims(subsystems)
        else
            nξ, ξids, ξxyids = 0, Tuple(0 for i in 1:np), Tuple(0 for i in 1:np)
        end
        new{ΔT,nx,nu,np,nξ,TS,typeof(xids),typeof(uids),
            typeof(xyids),typeof(ξids),typeof(ξxyids)}(subsystems, xids, uids,
                                                       xyids, ξids, ξxyids)
    end
end
n_players(cs::ProductSystem{ΔT,nx,nu,np}) where {ΔT,nx,nu,np} = np
subsystems(cs::ProductSystem) = cs.subsystems
xindex(cs::ProductSystem) = cs.xids
uindex(cs::ProductSystem) = cs.uids
xyindex(cs::ProductSystem) = cs.xyids

function dx(cs::ProductSystem, x::SVector, u::SVector, t::AbstractFloat)
    dxs = map(subsystems(cs), xindex(cs), uindex(cs)) do sub, xid, uid
        xᵢ = x[xid]
        uᵢ = u[uid]
        return dx(sub, xᵢ, uᵢ, t)
    end
    return vcat(dxs...)
end

# TODO: maybe make this a generated funtions. Depending on the dimensionality,
# it might be benificial to do the sparse or dense `next_x`.
function next_x(cs::ProductSystem, x::SVector, u::SVector, t::Float64)
    xs = map(subsystems(cs), xindex(cs), uindex(cs)) do sub, xid, uid
        xᵢ = x[xid]
        uᵢ = u[uid]
        return next_x(sub, xᵢ, uᵢ, t)
    end
    return vcat(xs...)
end

# computing large matrix exponentials is expensive. Therefore, we exploit the
# sparsity of the dynamics and compute the linearization for each subsystems
@inline function linearize_discrete(cs::ProductSystem, x::SVector, u::SVector, t::AbstractFloat)
    nx = n_states(cs)
    nu = n_controls(cs)
    # the full matrices
    A = @MMatrix zeros(nx, nx)
    B = @MMatrix zeros(nx, nu)
    for (i, sub) in enumerate(subsystems(cs))
        xid = xindex(cs)[i]
        uid = uindex(cs)[i]

        lin_sub = linearize_discrete(sub, x[xid], u[uid], t)
        A[xid, xid] = lin_sub.A
        B[xid, uid] = lin_sub.B
    end
    return LinearSystem{samplingtime(cs)}(SMatrix(A), SMatrix(B))
end

"------------------- Implement Feedback Linearization Interface -------------------"

# For now, the product system is only feedback linearizable if all subsystems are
LinearizationStyle(cs::ProductSystem) = common_linstyle(subsystems(cs))

n_linstates(cs::ProductSystem{ΔT,nx,nu,np,nξ}) where {ΔT,nx,nu,np,nξ} = nξ
ξxyindex(cs::ProductSystem) = cs.ξxyids
ξindex(cs::ProductSystem) = cs.ξids

function feedbacklin(cs::ProductSystem)
    nx = n_states(cs)
    nu = n_controls(cs)
    # the full matrices
    A = @MMatrix zeros(nx, nx)
    B = @MMatrix zeros(nx, nu)
    for (i, sub) in enumerate(subsystems(cs))
        xid = xindex(cs)[i]
        uid = uindex(cs)[i]

        lin_sub = feedbacklin(sub)
        A[xid, xid] = lin_sub.dyn.A
        B[xid, uid] = lin_sub.dyn.B
    end
    dyn = LinearSystem{samplingtime(cs)}(SMatrix(A), SMatrix(B))
    return LTISystem(dyn, ξxyindex(cs), ξindex(cs))
end

"A sparse version of the `feedbacklin` method that may be faster if the product
system contains of many subsystems (exploits block diagonal structure of linearized
system matrix)."
function sparse_feedbacklin(cs::ProductSystem)
    feedbacklin_subs = map(subsystems(cs)) do sub
        feedbacklin(sub)
    end
    return ProductSystem(feedbacklin_subs)
end

function x_from(cs::ProductSystem, ξ::SVector)
    subs = subsystems(cs)
    ξids = ξindex(cs)
    xs = map(pindex(cs)) do i
        x_from(subs[i], ξ[ξids[i]])
    end
    # Notice: we can safely vcat because in a product system, xids will always be
    # ordered in an ascending order.
    return vcat(Tuple(xs)...)
end

function ξ_from(cs::ProductSystem, x::SVector)
    subs = subsystems(cs)
    xids = xindex(cs)
    ξs = map(pindex(cs)) do i
        ξ_from(subs[i], x[xids[i]])
    end
    # Notice: we can safely vcat because in a product system, ξids will always be
    # ordered in an ascending order.
    return vcat(Tuple(ξs)...)
end

function λ_issingular(cs::ProductSystem, ξ::SVector)
    return any(zip(subsystems(cs), ξindex(cs))) do (subᵢ, ξis)
        λ_issingular(subᵢ, ξ[ξis])
    end
end

# TODO: why does this allocate?
function inverse_decoupling_matrix(cs::ProductSystem, x::SVector)
    nu = n_controls(cs)
    Minv = @MMatrix(zeros(nu, nu))
    for (sub, xid, uid) in zip(subsystems(cs), xindex(cs), uindex(cs))
        xᵢ = x[xid]
        Minv[uid, uid] = inverse_decoupling_matrix(sub, xᵢ)
    end
    return SMatrix(Minv)
end

function decoupling_drift_term(cs::ProductSystem, x)
    nu = n_controls(cs)
    subs = subsystems(cs)
    xids = xindex(cs)
    uids = uindex(cs)

    ms = map(pindex(cs)) do i
        decoupling_drift_term(subs[i], x[xids[i]])
    end
    # Notice: we can safely vcat because in a product system, uids will always be
    # ordered in an ascending order.
    return vcat(Tuple(ms)...)
end

function transformed_cost(cs::ProductSystem, c::PlayerCost, np::Int)
    # forward to the corresponding subsystem
    return transformed_cost(subsystems(cs)[player_id(c)], c, np)
end


text:
using testpkg2
using Test

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


text:
import DistributedFactorGraphs: getVariableType

"""
    CliqStatus
Clique status message enumerated type with status.
"""
@enum CliqStatus NULL NO_INIT INITIALIZED UPSOLVED MARGINALIZED DOWNSOLVED UPRECYCLED ERROR_STATUS

# Used for UPWARD_DIFFERENTIAL, UPWARD_COMMON, DOWNWARD_COMMON marginalized types
abstract type MessagePassDirection end
struct UpwardPass <: MessagePassDirection end
struct DownwardPass <: MessagePassDirection end

abstract type MessageType end
struct NonparametricMessage <: MessageType end
struct ParametricMessage <: MessageType end

abstract type PackedSamplableBelief end

const SamplableBelief = Union{Distributions.Distribution, KDE.BallTreeDensity, AMP.ManifoldKernelDensity, AliasingScalarSampler, FluxModelsDistribution, HeatmapGridDensity, LevelSetGridNormal}

#Supported types for parametric
const ParametricTypes = Union{Normal, MvNormal}

"""
    $TYPEDEF

INTERMEDIATE DATA STRUCTURE DURING REFACTORING.

Representation of the belief of a single variable.

Notes:
- we want to send the joint, this is just to resolve consolidation #459 first.
- Long term objective is single joint definition, likely called `LikelihoodMessage`.
"""
struct TreeBelief{T <: InferenceVariable, P, M <: MB.AbstractManifold}
  val::Vector{P}
  bw::Array{Float64,2}
  infoPerCoord::Vector{Float64}
  # see DFG #603, variableType defines the domain and manifold as well as group operations for a variable in the factor graph
  variableType::T
  # TODO -- DEPRECATE
  manifold::M # Tuple{Vararg{Symbol}} # NOTE added during #459 effort
  # only populated during up as solvableDims for each variable in clique, #910
  solvableDim::Float64
end

TreeBelief( p::ManifoldKernelDensity,
            ipc::AbstractVector{<:Real}=[0.0;],
            variableType::T=ContinuousScalar(),
            manifold=getManifold(variableType),
            solvableDim::Real=0) where {T <: InferenceVariable} = TreeBelief(getPoints(p), getBW(p), ipc, variableType, manifold, solvableDim)

TreeBelief( val::AbstractVector{P},
            bw::Array{Float64,2},
            ipc::AbstractVector{<:Real}=[0.0;],
            variableType::T=ContinuousScalar(),
            manifold::M=getManifold(variableType),
            solvableDim::Real=0) where {P, T <: InferenceVariable, M <:MB.AbstractManifold} = TreeBelief{T,P,M}(val, bw, ipc, variableType, manifold, solvableDim)

function TreeBelief(vnd::VariableNodeData{T}, solvDim::Real=0) where T
  TreeBelief( vnd.val, vnd.bw, vnd.infoPerCoord, getVariableType(vnd), getManifold(T), solvDim )
end

TreeBelief( vari::DFGVariable,
            solveKey::Symbol=:default;
            solvableDim::Real=0 ) = TreeBelief( getSolverData(vari, solveKey) , solvableDim)
#

getVariableType(tb::TreeBelief) = tb.variableType

getManifold(treeb::TreeBelief) = getManifold(treeb.variableType)

function compare(t1::TreeBelief, t2::TreeBelief)
  TP = true
  TP = TP && norm(t1.val - t2.val) < 1e-5
  TP = TP && norm(t1.bw - t2.bw) < 1e-5
  TP = TP && isapprox(t1.infoPerCoord, t2.infoPerCoord, atol= 1e-4)
  TP = TP && t1.variableType == t2.variableType
  TP = TP && abs(t1.solvableDim - t2.solvableDim) < 1e-5
  return TP
end

#


text:
# batch run converted into matrix operations,
# batch_size preallocation for nn layers

using LinearAlgebra
using MLDatasets
using Random

struct network_v3
    num_layers::Int64
    sizearr::Array{Int64,1}
    biases::Array{Array{Float64,1},1}
    weights::Array{Array{Float64,2},1}
    zs::Array{Array{Float64,2},1}
    activations::Array{Array{Float64,2},1}
    ∇_b::Array{Array{Float64,1},1}
    ∇_w::Array{Array{Float64,2},1}
    δs::Array{Array{Float64,2},1}
end

σ(z) = 1/(1+exp(-z)) # @inline ?
σ_grad(z) = σ(z)*(1-σ(z))

function (net::network_v3)(a)
    for (w, b) in zip(net.weights, net.biases)
        a = σ.(w*a .+ b)
    end

    return a
end

function network_v3(sizes, batch_size)
    num_layers = length(sizes)
    sizearr = sizes
    biases = [randn(y) for y in sizes[2:end]]
    weights = [randn(y, x) for (x, y) in zip(sizes[1:end-1], sizes[2:end])]
    zs = [zeros(y, batch_size) for y in sizes[2:end]]
    activations = [zeros(y, batch_size) for y in sizes[2:end]]
    ∇_b = [zeros(y) for y in sizes[2:end]]
    ∇_w = [zeros(y, x) for (x, y) in zip(sizes[1:end-1], sizes[2:end])]
    δs = [zeros(y, batch_size) for y in sizes[2:end]]
    network_v3(num_layers, sizearr, biases, weights, zs, activations,∇_b,∇_w,δs)
end

function update_batch(net::network_v3, x, y, η)
    ∇_b = net.∇_b
    ∇_w = net.∇_w

    # not need to fill ∇_b ∇_w with zero!

    backprop!(net, x, y)

    coef = η/size(x,2)

    for i in 1:length(∇_b)
        net.biases[i] .-= coef.*∇_b[i]
    end

    for i in 1:length(∇_w)
        net.weights[i] .-= coef.*∇_w[i]
    end
end

function backprop!(net::network_v3, x, y)
    ∇_b = net.∇_b
    ∇_w = net.∇_w

    len = net.num_layers - 1
    activations = net.activations
    zs = net.zs
    δs = net.δs

    input = x
    for i in 1:len
        b, w, z = net.biases[i], net.weights[i], zs[i]
        mul!(z, w, input) # z = w * input
        z .+= b
        activations[i] .= σ.(z)
        input = activations[i]
    end

    δ = δs[end]
    δ .= (activations[end] .- y) .* σ_grad.(zs[end])
    sum!(∇_b[end], δ)

    for l in 1:len-1
        mul!(∇_w[end-l+1], δ, activations[end-l]') # ∇_w[end-l+1] = δ * activations[end-l]'
        z = zs[end-l]
        mul!(δs[end-l], net.weights[end-l+1]', δ) # δs[end-l] = net.weights[end-l+1]' * δ
        δ = δs[end-l]
        δ .*= σ_grad.(z)
        sum!(∇_b[end-l], δ)
    end
    mul!(∇_w[1], δ, x') # ∇_w[1] = δ * x'

    return nothing
end

function SGDtrain(net::network_v3, traindata, epochs, batch_size, η, testdata=nothing)
    n_test = testdata != nothing ? size(testdata[1], 2) : nothing
    n = size(traindata[1], 2)

    idx = randperm(n) # one time shuffle for performance, then only take random batch index
    train_x = traindata[1][:,idx]
    train_y = traindata[2][:,idx]

    # reorganize data in batches
    batch = [(train_x[:, k-batch_size+1 : k], train_y[:, k-batch_size+1 : k]) for k in batch_size:batch_size:n]

    println("========")
    for j in 1:epochs
        idx = randperm(length(batch))

        @time for k in idx
            update_batch(net, batch[k]..., η)
        end

        if testdata != nothing
            println("Epoch ", j,": ", evaluate(net, testdata...), "/", n_test)
        else
            println("Epoch ", j," complete.")
        end
    end
end

function evaluate(net::network_v3, x, y)
    hits = 0
    for i = 1:size(x, 2)
        out = net(x[:,i])
        if (findmax(out)[2] - 1) == y[i]
            hits += 1
        end
    end
    hits
end

function loaddata(rng = 1:60000)
    train_x, train_y = MNIST.traindata(Float64, Vector(rng))
    train_x = reshape(train_x, size(train_x,1)*size(train_x,2), :) # 28 x 28 x N -> 28*28 x N
    train_y = vectorize(train_y)
    test_x, test_y = MNIST.testdata(Float64)
    test_x = reshape(test_x, size(test_x,1)*size(test_x,2), :) # 28 x 28 x N -> 28*28 x N
    return (train_x, train_y), (test_x, test_y)
end

function vectorize(vec)
    N = 10
    len = length(vec)
    mtx = zeros(N, len)
    for i = 1:len
        mtx[vec[i]+1, i] = 1
    end
    return mtx
end

function main_v3()
    epochs = 10
    batch_size = 10
    η = 1.25
    net = network_v3([784, 30, 10], batch_size)
    traindata, testdata = loaddata()
    SGDtrain(net, traindata, epochs, batch_size, η, testdata)
	# @profiler SGDtrain(net, traindata, 1, batch_size, η, testdata)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""Represents storage that is managed by an external CSI volume driver (Beta feature)

    IoK8sApiCoreV1CSIPersistentVolumeSource(;
        controllerExpandSecretRef=nothing,
        controllerPublishSecretRef=nothing,
        driver=nothing,
        fsType=nothing,
        nodePublishSecretRef=nothing,
        nodeStageSecretRef=nothing,
        readOnly=nothing,
        volumeAttributes=nothing,
        volumeHandle=nothing,
    )

    - controllerExpandSecretRef::IoK8sApiCoreV1SecretReference : ControllerExpandSecretRef is a reference to the secret object containing sensitive information to pass to the CSI driver to complete the CSI ControllerExpandVolume call. This is an alpha field and requires enabling ExpandCSIVolumes feature gate. This field is optional, and may be empty if no secret is required. If the secret object contains more than one secret, all secrets are passed.
    - controllerPublishSecretRef::IoK8sApiCoreV1SecretReference : ControllerPublishSecretRef is a reference to the secret object containing sensitive information to pass to the CSI driver to complete the CSI ControllerPublishVolume and ControllerUnpublishVolume calls. This field is optional, and may be empty if no secret is required. If the secret object contains more than one secret, all secrets are passed.
    - driver::String : Driver is the name of the driver to use for this volume. Required.
    - fsType::String : Filesystem type to mount. Must be a filesystem type supported by the host operating system. Ex. \"ext4\", \"xfs\", \"ntfs\".
    - nodePublishSecretRef::IoK8sApiCoreV1SecretReference : NodePublishSecretRef is a reference to the secret object containing sensitive information to pass to the CSI driver to complete the CSI NodePublishVolume and NodeUnpublishVolume calls. This field is optional, and may be empty if no secret is required. If the secret object contains more than one secret, all secrets are passed.
    - nodeStageSecretRef::IoK8sApiCoreV1SecretReference : NodeStageSecretRef is a reference to the secret object containing sensitive information to pass to the CSI driver to complete the CSI NodeStageVolume and NodeStageVolume and NodeUnstageVolume calls. This field is optional, and may be empty if no secret is required. If the secret object contains more than one secret, all secrets are passed.
    - readOnly::Bool : Optional: The value to pass to ControllerPublishVolumeRequest. Defaults to false (read/write).
    - volumeAttributes::Dict{String, String} : Attributes of the volume to publish.
    - volumeHandle::String : VolumeHandle is the unique volume name returned by the CSI volume plugin’s CreateVolume to refer to the volume on all subsequent calls. Required.
"""
mutable struct IoK8sApiCoreV1CSIPersistentVolumeSource <: SwaggerModel
    controllerExpandSecretRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1SecretReference } # spec name: controllerExpandSecretRef
    controllerPublishSecretRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1SecretReference } # spec name: controllerPublishSecretRef
    driver::Any # spec type: Union{ Nothing, String } # spec name: driver
    fsType::Any # spec type: Union{ Nothing, String } # spec name: fsType
    nodePublishSecretRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1SecretReference } # spec name: nodePublishSecretRef
    nodeStageSecretRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1SecretReference } # spec name: nodeStageSecretRef
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly
    volumeAttributes::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: volumeAttributes
    volumeHandle::Any # spec type: Union{ Nothing, String } # spec name: volumeHandle

    function IoK8sApiCoreV1CSIPersistentVolumeSource(;controllerExpandSecretRef=nothing, controllerPublishSecretRef=nothing, driver=nothing, fsType=nothing, nodePublishSecretRef=nothing, nodeStageSecretRef=nothing, readOnly=nothing, volumeAttributes=nothing, volumeHandle=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1CSIPersistentVolumeSource, Symbol("controllerExpandSecretRef"), controllerExpandSecretRef)
        setfield!(o, Symbol("controllerExpandSecretRef"), controllerExpandSecretRef)
        validate_property(IoK8sApiCoreV1CSIPersistentVolumeSource, Symbol("controllerPublishSecretRef"), controllerPublishSecretRef)
        setfield!(o, Symbol("controllerPublishSecretRef"), controllerPublishSecretRef)
        validate_property(IoK8sApiCoreV1CSIPersistentVolumeSource, Symbol("driver"), driver)
        setfield!(o, Symbol("driver"), driver)
        validate_property(IoK8sApiCoreV1CSIPersistentVolumeSource, Symbol("fsType"), fsType)
        setfield!(o, Symbol("fsType"), fsType)
        validate_property(IoK8sApiCoreV1CSIPersistentVolumeSource, Symbol("nodePublishSecretRef"), nodePublishSecretRef)
        setfield!(o, Symbol("nodePublishSecretRef"), nodePublishSecretRef)
        validate_property(IoK8sApiCoreV1CSIPersistentVolumeSource, Symbol("nodeStageSecretRef"), nodeStageSecretRef)
        setfield!(o, Symbol("nodeStageSecretRef"), nodeStageSecretRef)
        validate_property(IoK8sApiCoreV1CSIPersistentVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        validate_property(IoK8sApiCoreV1CSIPersistentVolumeSource, Symbol("volumeAttributes"), volumeAttributes)
        setfield!(o, Symbol("volumeAttributes"), volumeAttributes)
        validate_property(IoK8sApiCoreV1CSIPersistentVolumeSource, Symbol("volumeHandle"), volumeHandle)
        setfield!(o, Symbol("volumeHandle"), volumeHandle)
        o
    end
end # type IoK8sApiCoreV1CSIPersistentVolumeSource

const _property_map_IoK8sApiCoreV1CSIPersistentVolumeSource = Dict{Symbol,Symbol}(Symbol("controllerExpandSecretRef")=>Symbol("controllerExpandSecretRef"), Symbol("controllerPublishSecretRef")=>Symbol("controllerPublishSecretRef"), Symbol("driver")=>Symbol("driver"), Symbol("fsType")=>Symbol("fsType"), Symbol("nodePublishSecretRef")=>Symbol("nodePublishSecretRef"), Symbol("nodeStageSecretRef")=>Symbol("nodeStageSecretRef"), Symbol("readOnly")=>Symbol("readOnly"), Symbol("volumeAttributes")=>Symbol("volumeAttributes"), Symbol("volumeHandle")=>Symbol("volumeHandle"))
const _property_types_IoK8sApiCoreV1CSIPersistentVolumeSource = Dict{Symbol,String}(Symbol("controllerExpandSecretRef")=>"IoK8sApiCoreV1SecretReference", Symbol("controllerPublishSecretRef")=>"IoK8sApiCoreV1SecretReference", Symbol("driver")=>"String", Symbol("fsType")=>"String", Symbol("nodePublishSecretRef")=>"IoK8sApiCoreV1SecretReference", Symbol("nodeStageSecretRef")=>"IoK8sApiCoreV1SecretReference", Symbol("readOnly")=>"Bool", Symbol("volumeAttributes")=>"Dict{String, String}", Symbol("volumeHandle")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1CSIPersistentVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1CSIPersistentVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1CSIPersistentVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1CSIPersistentVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1CSIPersistentVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1CSIPersistentVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1CSIPersistentVolumeSource)
    (getproperty(o, Symbol("driver")) === nothing) && (return false)
    (getproperty(o, Symbol("volumeHandle")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1CSIPersistentVolumeSource }, name::Symbol, val)
end


text:
#######################################
#
# Filtering
#
#######################################

"""
    filter_highpass(a::EEG; cutOff::Real=2, fs::Real=samplingrate(a), order::Int=0, kwargs...)

Apply a high pass filter.

If order is zero then its computed magically.

# Examples
```julia
a = read_EEG(fname)
b = filter_highpass(a)
c = filter_highpass(a, cutOff = 1)
```
"""
function filter_highpass(
    a::EEG;
    cutOff::Union{typeof(1.0u"Hz"),typeof(1u"Hz")} = 2u"Hz",
    order::Int = 0,
    phase::AbstractString = "zero",
    kwargs...,
)
    Wn = (cutOff |> u"Hz" |> ustrip) / (samplingrate(Float64, a) / 2)
    if order == 0
        order = default_fir_filterorder(Highpass(Wn), samplingrate(Float64, a)) + 1
    end
    f = digitalfilter(Highpass(Wn), FIRWindow(DSP.hamming(order)))
    a = filter(a, f, phase)
end

"""
    filter_lowpass(a::EEG; cutOff::Real=150, fs::Real=samplingrate(a), order::Int=0, kwargs...)

Apply a low pass filter.

If order is zero then its computed magically.

# Examples
```julia
a = read_EEG(fname)
b = filter_lowpass(a)
c = filter_lowpass(a, cutOff = 1)
```
"""
function filter_lowpass(
    a::EEG;
    cutOff::Union{typeof(1.0u"Hz"),typeof(1u"Hz")} = 150u"Hz",
    order::Int = 0,
    phase::AbstractString = "zero",
    kwargs...,
)
    Wn = (cutOff |> u"Hz" |> ustrip) / (samplingrate(Float64, a) / 2)
    if order == 0
        order = default_fir_filterorder(Lowpass(Wn), samplingrate(Float64, a)) + 1
    end
    f = digitalfilter(Lowpass(Wn), FIRWindow(DSP.hamming(order)))
    a = filter(a, f, phase)
end

"""
    filter_bandpass(a::EEG, lower, upper; n::Int=24, rp::Number = 0.0001, kwargs...)

Apply a band pass filter.

# Examples
```julia
a = read_EEG(fname)
a = filter_bandpass(a)
```
"""
function filter_bandpass(
    a::EEG,
    lower::Union{typeof(1.0u"Hz"),typeof(1u"Hz")},
    upper::Union{typeof(1.0u"Hz"),typeof(1u"Hz")};
    n::Int = 24,
    rp::Number = 0.0001,
    phase::AbstractString = "zero-double",
    kwargs...,
)

    Wn_lower = (lower |> u"Hz" |> ustrip) / (samplingrate(Float64, a) / 2)
    Wn_upper = (upper |> u"Hz" |> ustrip) / (samplingrate(Float64, a) / 2)

    f = digitalfilter(Bandpass(Wn_lower, Wn_upper), Chebyshev1(n, rp))
    a = filter(a, f, phase)
end


text:
## Illustration in 2D
# note: behavior is generally less consistent in lower dimensions but 2D is convenient for illustration/visualization

using HePPCAT
using CairoMakie, StableRNGs

## Setup
d = 2
n = [400,100]
v = [4,0.04]
F = ones(d,1)/sqrt(d)

## Generate data
rng = StableRNG(0)
X = [F*randn(rng,nl)' + sqrt(vl)*randn(rng,d,nl) for (nl,vl) in zip(n,v)]
# X = reshape.(collect.(eachcol(reduce(hcat,X))),:,1) # flatten to give each sample its own noise variance estimate

## Compute estimates
models = [
    "(Homoscedastic) PPCA"           => heppcat(X,1,0),    # init is homoscedastic PPCA
    "Heteroscedastic PPCA (HePPCAT)" => heppcat(X,1,1000)
]

## Plot
scene, layout = layoutscene(20,resolution=(800,600))
datacolors = [:dodgerblue1,:green3]
for (idx,(title,M)) in enumerate(models)
    # Data and axes
    ax = layout[1,idx] = LAxis(scene,title=title)

    # Latent axis
    text!(ax,"true component",position=(5.7,5.7),textsize=0.55,align=(:right,:top),rotation=pi/4)
    lines!(ax,-6.25..6.25,identity,linewidth=3,color=:black)

    # Data points
    for (Xl,color) in zip(X,datacolors)
        scatter!(ax,Xl[1,:],Xl[2,:],strokewidth=0.2,markersize=6,color=color)
    end

    # Estimate
    text!(ax,"estimate",position=(M.U[1],M.U[2]).*([4.5,5.7][idx]*sqrt(2)),textsize=0.55,align=(:right,:bottom),rotation=atan(M.U[2]/M.U[1]),color=:darkorange1)
    lines!(ax,-6.25..6.25,x->x*M.U[2]/M.U[1],linewidth=3,color=:darkorange1)

    # Formatting
    ax.autolimitaspect = 1
    ax.xticks = 0:0
    ax.yticks = 0:0
    limits!(ax,(-6.25,6.25),(-6.25,6.25))
    hidedecorations!(ax,grid=false)

    # Noise variance estimates
    ax = layout[2,idx] = LAxis(scene)
    scatter!(ax,reduce(vcat,fill.(M.v,n)),color=reduce(vcat,fill.(datacolors,n)),strokecolor=:transparent,markersize=2)
    # scatter!(ax,fill(-1,sum(n)),color=reduce(vcat,fill.(datacolors,n)),strokecolor=:transparent,markersize=5)
    # scatter!(ax,reduce(vcat,fill.(v,n)),color=:black,strokecolor=:transparent,markersize=5)
    # scatter!(ax,reduce(vcat,fill.(M.v,n)),color=:darkorange1,strokecolor=:transparent,markersize=3)
    ax.xticks = [0,n[1],sum(n)]
    ax.yticks = v
    limits!(ax,(0,sum(n)),(-1,5))
    hidedecorations!(ax,grid=false)
end

# Legend
leg = layout[end+1,:] = LLegend(scene,
    [MarkerElement(marker=:circle,color=c,strokecolor=:transparent) for c in datacolors],
    ["Noisier samples","Cleaner samples"]
)
leg.labelsize = 15
leg.tellheight = true
leg.orientation = :horizontal
leg.framevisible = false

# Format and save
layout[1,1,Left()] = LText(scene,"data and estimated components",textsize=14,rotation=pi/2)
layout[2,1,Left()] = LText(scene,"est. noise\nvariances",textsize=14,rotation=pi/2)
rowsize!(layout, 2, Relative(1/6))

CairoMakie.save("illustration-2D.png",scene)
scene


text:
using Test
using QuantumOpticsBase
using LinearAlgebra, Random

@testset "operators-lazyproduct" begin

Random.seed!(0)

D(op1::AbstractOperator, op2::AbstractOperator) = abs(tracedistance_nh(dense(op1), dense(op2)))
D(x1::StateVector, x2::StateVector) = norm(x2-x1)

b1a = GenericBasis(2)
b1b = GenericBasis(3)
b2a = GenericBasis(1)
b2b = GenericBasis(4)
b3a = GenericBasis(1)
b3b = GenericBasis(5)

b_l = b1a⊗b2a⊗b3a
b_r = b1b⊗b2b⊗b3b

# Test creation
@test_throws ArgumentError LazyProduct()
@test_throws QuantumOpticsBase.IncompatibleBases LazyProduct(randoperator(b_l, b_r), randoperator(b_l, b_r))
@test_throws QuantumOpticsBase.IncompatibleBases LazyProduct(randoperator(b_l, b_r), sparse(randoperator(b_l, b_r)))

# Test copy
op1 = 2*LazyProduct(randoperator(b_l, b_r), sparse(randoperator(b_r, b_l)))
op2 = copy(op1)
@test op1 == op2
@test !(op1 === op2)
op2.operators[1].data[1,1] = complex(10.)
@test op1.operators[1].data[1,1] != op2.operators[1].data[1,1]
op2.factor = 3.
@test op2.factor != op1.factor

# Test dense & sparse
op1 = randoperator(b_l, b_r)
op2 = randoperator(b_r, b_l)
@test 0.1*(op1*op2) == dense(LazyProduct([sparse(op1), sparse(op2)], 0.1))
@test 0.1*(sparse(op1)*sparse(op2)) == sparse(LazyProduct([op1, op2], 0.1))

# Arithmetic operations
# =====================
op1a = randoperator(b_l, b_r)
op1b = randoperator(b_r, b_l)
op2a = randoperator(b_l, b_r)
op2b = randoperator(b_r, b_l)
op3a = randoperator(b_l, b_l)
op1 = LazyProduct([op1a, sparse(op1b)])*0.1
op1_ = 0.1*(op1a*op1b)
op2 = LazyProduct([sparse(op2a), op2b], 0.3)
op2_ = 0.3*(op2a*op2b)
op3 = LazyProduct(op3a)
op3_ = op3a

x1 = Ket(b_l, rand(ComplexF64, length(b_l)))
x2 = Ket(b_l, rand(ComplexF64, length(b_l)))
xbra1 = Bra(b_l, rand(ComplexF64, length(b_l)))
xbra2 = Bra(b_l, rand(ComplexF64, length(b_l)))

# Addition
@test_throws ArgumentError op1 + op2
@test 1e-14 > D(-op1_, -op1)

# Test multiplication
@test_throws DimensionMismatch op1a*op1a
@test 1e-11 > D(op1*(x1 + 0.3*x2), op1_*(x1 + 0.3*x2))
@test 1e-11 > D((xbra1 + 0.3*xbra2)*op1, (xbra1 + 0.3*xbra2)*op1_)
@test 1e-11 > D(op1*x1 + 0.3*op1*x2, op1_*x1 + 0.3*op1_*x2)
@test 1e-12 > D(dagger(x1)*dagger(0.3*op2), dagger(x1)*dagger(0.3*op2_))
@test 0.3*LazyProduct(op1, sparse(op2)) == LazyProduct([op1, sparse(op2)], 0.3)
@test 0.3*LazyProduct(op1)*LazyProduct(sparse(op2)) == LazyProduct([op1, sparse(op2)], 0.3)

# Test division
@test 1e-14 > D(op1/7, op1_/7)

# Test identityoperator
Idense = identityoperator(DenseOpType, b_l)
id = identityoperator(LazyProduct, b_l)
@test isa(id, LazyProduct)
@test dense(id) == Idense
@test 1e-11 > D(id*x1, x1)
@test 1e-11 > D(xbra1*id, xbra1)

# Test tr and normalize
op1 = randoperator(b_l)
op2 = randoperator(b_l)
op = LazyProduct(op1, op2)
@test_throws ArgumentError tr(op)
@test_throws ArgumentError ptrace(op, [1, 2])
@test_throws ArgumentError normalize(op)
@test_throws ArgumentError normalize!(op)

# Test expect
op1 = randoperator(b_l)
op2 = randoperator(b_l)
op = 0.3*LazyProduct(op1, sparse(op2))
op_ = 0.3*op1*op2

state = Ket(b_l, rand(ComplexF64, length(b_l)))
@test expect(op, state) ≈ expect(op_, state)

state = DenseOperator(b_l, b_l, rand(ComplexF64, length(b_l), length(b_l)))
@test expect(op, state) ≈ expect(op_, state)

# Permute systems
op1 = randoperator(b_l)
op2 = randoperator(b_l)
op3 = randoperator(b_l)
op = 0.3*LazyProduct(op1, op2, sparse(op3))
op_ = 0.3*op1*op2*op3

@test 1e-14 > D(permutesystems(op, [1, 3, 2]), permutesystems(op_, [1, 3, 2]))
@test 1e-14 > D(permutesystems(op, [2, 1, 3]), permutesystems(op_, [2, 1, 3]))
@test 1e-14 > D(permutesystems(op, [2, 3, 1]), permutesystems(op_, [2, 3, 1]))
@test 1e-14 > D(permutesystems(op, [3, 1, 2]), permutesystems(op_, [3, 1, 2]))
@test 1e-14 > D(permutesystems(op, [3, 2, 1]), permutesystems(op_, [3, 2, 1]))

# Test gemv
op1 = randoperator(b_l, b_r)
op2 = randoperator(b_r, b_l)
op3 = randoperator(b_l, b_r)
op = LazyProduct([op1, sparse(op2), op3], 0.2)
op_ = 0.2*op1*op2*op3

state = Ket(b_r, rand(ComplexF64, length(b_r)))
result_ = Ket(b_l, rand(ComplexF64, length(b_l)))
result = copy(result_)
QuantumOpticsBase.mul!(result,op,state,complex(1.),complex(0.))
@test 1e-11 > D(result, op_*state)

result = copy(result_)
alpha = complex(1.5)
beta = complex(2.1)
QuantumOpticsBase.mul!(result,op,state,alpha,beta)
@test 1e-11 > D(result, alpha*op_*state + beta*result_)

state = Bra(b_l, rand(ComplexF64, length(b_l)))
result_ = Bra(b_r, rand(ComplexF64, length(b_r)))
result = copy(result_)
QuantumOpticsBase.mul!(result,state,op,complex(1.),complex(0.))
@test 1e-11 > D(result, state*op_)

result = copy(result_)
alpha = complex(1.5)
beta = complex(2.1)
QuantumOpticsBase.mul!(result,state,op,alpha,beta)
@test 1e-11 > D(result, alpha*state*op_ + beta*result_)

# Test gemm
op1 = randoperator(b_l, b_r)
op2 = randoperator(b_r, b_l)
op3 = randoperator(b_l, b_r)
op = LazyProduct([op1, sparse(op2), op3], 0.2)
op_ = 0.2*op1*op2*op3

state = randoperator(b_r, b_r)
result_ = randoperator(b_l, b_r)
result = copy(result_)
QuantumOpticsBase.mul!(result,op,state,complex(1.),complex(0.))
@test 1e-11 > D(result, op_*state)

result = copy(result_)
alpha = complex(1.5)
beta = complex(2.1)
QuantumOpticsBase.mul!(result,op,state,alpha,beta)
@test 1e-11 > D(result, alpha*op_*state + beta*result_)

state = randoperator(b_l, b_l)
result_ = randoperator(b_l, b_r)
result = copy(result_)
QuantumOpticsBase.mul!(result,state,op,complex(1.),complex(0.))
@test 1e-11 > D(result, state*op_)

result = copy(result_)
alpha = complex(1.5)
beta = complex(2.1)
QuantumOpticsBase.mul!(result,state,op,alpha,beta)
@test 1e-11 > D(result, alpha*state*op_ + beta*result_)

end # testset


text:
export gumbel_softmax

"""
    gumbel_softmax(x::Var, tau)

* x: logarithm of probability distribution, e.g., output of `logsoftmax`.
* tau: temprature
"""
function gumbel_softmax(x::Var, tau)
    tau > 0.0 || throw("tau: $(tau) must be positive.")
    tau = eltype(x.data)(tau)
    y = Var(nothing, gumbel_softmax, (x,tau))
    g = gumbel_softmax(x.data, tau)
    y.data = vec(argmax(g,1))
    y.grad = zeros(g)
    y.df! = () -> begin
        isvoid(x.grad) && return
        ∇gumbel_softmax!(g, y.grad, x.grad, tau)
    end
    y
end

gumbel{T}(u::T) = -log(-log(u+1e-20)+1e-20)

function gumbel_softmax{T}(x::Matrix{T}, tau::T)
    y = similar(x)
    u = rand(T, size(x))
    tau = T(1) / tau
    @inbounds for i = 1:length(x)
        y[i] = (x[i] + gumbel(u[i])) * tau
    end
    softmax(y)
end

function ∇gumbel_softmax!{T}(y::Vector{Int}, gy::Matrix{T}, work::Matrix{T}, gx::Matrix{T}, tau::T)
    for j = 1:size(y,2)
        sum = T(0)
        for i = 1:size(y,1)
            sum += gy[i,j] * y[i,j]
        end
        y[j]
        for i = 1:size(y,1)
            gx[i,j] += y[i,j] * (gy[i,j] - sum)
        end
    end
end


text:
# ==============================================================================
# Five-dimensional system
#
# system type: continuous LTI system
# state dimension: 5
# input dimension: 1
#
# This five-dimensional system is taken from Example 4.1 in [1].
#
# [1] Colas Le Guernic: Reachability analysis of hybrid systems with linear
#     continuous dynamics. PhD thesis, 2009.
# ==============================================================================
using MathematicalSystems, LazySets

function fiveDimSys_model()
    # system matrix
    D = [-1.0 -4.0  0.0  0.0  0.0;
          4.0 -1.0  0.0  0.0  0.0;
          0.0  0.0 -3.0  1.0  0.0;
          0.0  0.0 -1.0 -3.0  0.0;
          0.0  0.0  0.0  0.0 -2.0]
    P = [ 0.6 -0.1  0.1  0.7 -0.2;
         -0.5  0.7 -0.1 -0.8  0.0;
          0.9 -0.5  0.3 -0.6  0.1;
          0.5 -0.7  0.5  0.6  0.3;
          0.8  0.7  0.6 -0.3  0.2]
    A = P * D * inv(P)

    # state domain
    X = Universe(5)

    # input domain
    U = Ball2(zeros(5), 0.01)

    # continuous LTI system
    S = @system(x' = Ax + u, x ∈ X, u ∈ U)

    return S
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
using MiniPETSc

using Test

if MPI.Comm_size(MPI.COMM_WORLD) == 1
    include("test_Mat.jl")
    include("test_Vec.jl")
    include("test_KSP.jl")
else
    include("test_GhostedVec.jl")
    include("test_parallel_Vec.jl")
end


text:
# for the DCHSBM representation
function linear_density(d_n::Tuple{AbstractDict{<:Integer, <:AbstractDict}, UnitRange{<:Integer}})
    d = d_n[1]
    edge_count = sum([sum(values(x)) for (k,x) in d])
    n = length(d_n[2])
    kmax = maximum(keys(d))
    possible_edges = sum([binomial(n,k) for k in 1:kmax])
    edge_count/possible_edges
end
function log_density(d_n::Tuple{AbstractDict{<:Integer, <:AbstractDict}, UnitRange{<:Integer}})
    size = hypergraphsize(d_n)
    size == 0 && return NaN
    n = length(d_n[2])
    log(size)/log(n)
end
function FBD.hypergraphsize(d_n::Tuple{AbstractDict{<:Integer, <:AbstractDict}, UnitRange{<:Integer}})
    sum([k*sum(values(x)) for (k,x) in d_n[1]])
end

# for the Kronecker representation
function linear_density(hedges::NTuple{3, <:AbstractVector{<:Integer}})
    edge_count = length(hedges[1])
    edge_count == 0 && return NaN
    n = maximum(maximum.(hedges))
    kmax = 3
    possible_edges = sum([binomial(n,k) for k in 1:kmax])
    edge_count/possible_edges
end
function log_density(hedges::NTuple{3, <:AbstractVector{<:Integer}})
    size = hypergraphsize(hedges)
    size == 0 && return NaN
    n = maximum(maximum.(hedges))
    log(size)/log(n)
end

# for the FBD representation
function linear_density(g::AbstractVector)
    edge_count = length(g)
    edge_count == 0 && return NaN
    n = maximum(maximum.(g))
    kmax = maximum(length.(g))
    possible_edges = sum([binomial(n,k) for k in 1:kmax])
    edge_count/possible_edges
end
function log_density(g::AbstractVector)
    size = hypergraphsize(g)
    size == 0 && return NaN
    n = maximum(maximum.(g))
    log(size)/log(n)
end


text:
convert(::Type{Tile}, md::Markdown.MD) = blocktile(md)

blocktile(x) = convert(Tile, x)
inlinetile(x) = convert(Tile, x)
blocktile(x::Tile) = x
inlinetile(x::Tile) = x

blocktile(md::Markdown.MD) = vbox(map(blocktile, md.content))
inlinetile(x::AbstractString) = x

blocktile{n}(md::Markdown.Header{n}) =
    class("md-heading", heading(n, map(inlinetile, md.text)))

inlinetile(md::Markdown.Code) = code(md.code, language=md.language)
blocktile(md::Markdown.Code)  =
    class("md-codeblock", codemirror(md.code, language=md.language, linenumbers=false, readonly=true))

blocktile(md::Markdown.BlockQuote) =
    class("md-blockquote", blockquote(map(blocktile, md.content)))
blocktile(md::Markdown.List) =
    class("md-list",
        list(map(item -> class("md-list-item", map(inlinetile, item)),
            md.items), ordered=md.ordered))
blocktile(md::Markdown.Paragraph) =
    class("md-paragraph", map(inlinetile, md.content))

inlinetile(md::Markdown.Italic) = emph(map(inlinetile, md.text))
inlinetile(md::Markdown.Bold) = fontweight(bold, map(inlinetile, md.text))
inlinetile(md::Markdown.Link) = link(md.url, map(inlinetile, md.text))

inlinetile(md::Markdown.Image) = image(md.url, alt=md.alt)
blocktile(md::Markdown.Image) = class("md-img", image(md.url, alt=md.alt))

inlinetile(md::Markdown.LaTeX) = tex(md.formula, block=false)
blocktile(md::Markdown.LaTeX) = class("md-tex", tex(md.formula, block=true))


text:
# A quadratic segment in passes through three points: 𝐱₁, 𝐱₂, and 𝐱₃.
# The assumed relation of the points may be seen in the diagram below.
#                 ___𝐱₃___
#            ____/        \____
#        ___/                  \
#     __/                       𝐱₂
#   _/
#  /
# 𝐱₁
#
# NOTE: 𝐱₃ is not necessarily the midpoint, or even between 𝐱₁ and 𝐱₂.
# q(r) = (2r-1)(r-1)𝐱₁ + r(2r-1)𝐱₂ + 4r(1-r)𝐱₃
# See The Visualization Toolkit: An Object-Oriented Approach to 3D Graphics, 4th Edition
# Chapter 8, Advanced Data Representation, in the interpolation functions section
struct QuadraticSegment{F <: AbstractFloat} <: Edge{F}
    points::SVector{3, Point{F}}
end

# # Constructors
# # -------------------------------------------------------------------------------------------------
# QuadraticSegment(p₁::Point, p₂::Point, p₃::Point) = QuadraticSegment(SVector(p₁, p₂, p₃))
#
# # Methods
# # -------------------------------------------------------------------------------------------------
# # Interpolation
# # q(0) = q[1], q(1) = q[2], q(1//2) = q[3]
# function (q::QuadraticSegment{F})(r::Real) where {F <: AbstractFloat}
#     # See The Visualization Toolkit: An Object-Oriented Approach to 3D Graphics, 4th Edition
#     # Chapter 8, Advanced Data Representation, in the interpolation functions section
#     rₜ = F(r)
#     return (2rₜ-1)*( rₜ-1)q[1] + rₜ*(2rₜ-1)q[2] + 4rₜ*( 1-rₜ)q[3]
# end
#
# arclength(q::QuadraticSegment) = arclength(q, Val(15))
# function arclength(q::QuadraticSegment{F}, ::Val{N}) where {F <: AbstractFloat, N}
#     # Numerical integration is used.
#     # (Gauss-Legengre quadrature)
#     #     1                  N
#     # L = ∫ ||∇ q⃗(r)||dr  ≈  ∑ wᵢ||∇ q⃗(rᵢ)||
#     #     0                 i=1
#     #
#     w, r = gauss_legendre_quadrature(F, Val(N))
#     return sum(@. w * norm(∇(q, r)))
# end
#
# # Find the axis-aligned bounding box of the segment.
# function boundingbox(q::QuadraticSegment)
#     # Find the r coordinates where ∂x/∂r = 0, ∂y/∂r = 0
#     # We know ∇ q, so we can directly compute these values
#     r_x = (3q[1].x + q[2].x - 4q[3].x)/(4(q[1].x + q[2].x - 2q[3].x))
#     if 0 < r_x < 1
#         x_extreme = (2r_x-1)*(r_x-1)q[1].x + r_x*(2r_x-1)q[2].x + 4r_x*(1-r_x)q[3].x
#         xmin = min(q[1].x, q[2].x, x_extreme)
#         xmax = max(q[1].x, q[2].x, x_extreme)
#     else
#         xmin = min(q[1].x, q[2].x)
#         xmax = max(q[1].x, q[2].x)
#     end
#
#     r_y = (3q[1].y + q[2].y - 4q[3].y)/(4(q[1].y + q[2].y - 2q[3].y))
#     if 0 < r_y < 1
#         y_extreme = (2r_y-1)*(r_y-1)q[1].y + r_y*(2r_y-1)q[2].y + 4r_y*(1-r_y)q[3].y
#         ymin = min(q[1].y, q[2].y, y_extreme)
#         ymax = max(q[1].y, q[2].y, y_extreme)
#     else
#         ymin = min(q[1].y, q[2].y)
#         ymax = max(q[1].y, q[2].y)
#     end
#     return Rectangle(Point(xmin, ymin), Point(xmax, ymax))
# end
#
# # Return the gradient of q, evalutated at r
# function gradient(q::QuadraticSegment{F}, r::Real) where {F <: AbstractFloat}
#     #∇q =
#     rₜ = F(r)
#     return (4rₜ - 3)*(q[1] - q[3]) + (4rₜ - 1)*(q[2] - q[3])
# end
#
# # Return if the point is left of the quadratic segment
# #   p    ^
# #   ^   /
# # v⃗ |  / u⃗
# #   | /
# #   o
# # function isleft(p::Point, q::QuadraticSegment)
# #     if isstraight(q) || p ∉  boundingbox(q)
# #         # We don't need to account for the curve if q is straight or p is outside
# #         # q's bounding box
# #         u⃗ = q[2] - q[1]
# #         v⃗ = p - q[1]
# #         return u⃗ × v⃗ > 0
# #     else
# #         # Get the nearest point on q to p.
# #         # Construct vectors from a point on q, close to p_near, to p_near and p.
# #         # Use the cross product of these vectors to determine if p isleft.
# #         r, p_near = nearest_point(p, q)
# #
# #         if r < 1e-6 || 1 < r # If r is small or beyond the valid range, just use q[2]
# #             u⃗ = q[2] - q[1]
# #             v⃗ = p - q[1]
# #         else # otherwise use a point on q, close to p_near
# #             q_base = q(0.95r)
# #             u⃗ = p_near - q_base
# #             v⃗ = p - q_base
# #         end
# #         return u⃗ × v⃗ > 0
# #     end
# # end
#
# # If the quadratic segment is effectively linear
# @inline function isstraight(q::QuadraticSegment)
#     # u⃗ × v⃗ = |u⃗||v⃗|sinθ
#     return abs((q[3] - q[1]) × (q[2] - q[1])) < 1e-8
# end
#
# # function intersect(l::LineSegment, q::QuadraticSegment)
# #     ϵ = parametric_coordinate_ϵ
# #     if isstraight(q) # Use line segment intersection.
# #         # See LineSegment for the math behind this.
# #         v⃗ = l[2] - l[1]
# #         u⃗ = q[2] - q[1]
# #         vxu = v⃗ × u⃗
# #         # Parallel or collinear lines, return.
# #         1e-8 < abs(vxu) || return (0x00000000, SVector(Point(), Point()))
# #         w⃗ = q[1] - l[1]
# #         # Delay division until r,s are verified
# #         if 0 <= vxu
# #             lowerbound = (-ϵ)vxu
# #             upperbound = (1 + ϵ)vxu
# #         else
# #             upperbound = (-ϵ)vxu
# #             lowerbound = (1 + ϵ)vxu
# #         end
# #         r_numerator = w⃗ × u⃗
# #         s_numerator = w⃗ × v⃗
# #         if (lowerbound ≤ r_numerator ≤ upperbound) && (lowerbound ≤ s_numerator ≤ upperbound)
# #             return (0x00000001, SVector(l(s_numerator/vxu), Point()))
# #         else
# #             return (0x00000000, SVector(Point(), Point()))
# #         end
# #     else
# #         # q(r) = (2r-1)(r-1)𝐱₁ + r(2r-1)𝐱₂ + 4r(1-r)𝐱₃
# #         # q(r) = 2r²(𝐱₁ + 𝐱₂ - 2𝐱₃) + r(-3𝐱₁ - 𝐱₂ + 4𝐱₃) + 𝐱₁
# #         # Let D⃗ = 2(𝐱₁ + 𝐱₂ - 2𝐱₃), E⃗ = (-3𝐱₁ - 𝐱₂ + 4𝐱₃), F⃗ = x₁
# #         # q(r) = r²D⃗ + rE⃗ + F⃗
# #         # l(s) = 𝐱₄ + sw⃗
# #         # If D⃗ × w⃗ ≠ 0
# #         #   𝐱₄ + sw⃗ = r²D⃗ + rE⃗ + F⃗
# #         #   sw⃗ = r²D⃗ + rE⃗ + (F⃗ - 𝐱₄)
# #         #   0 = r²(D⃗ × w⃗) + r(E⃗ × w⃗) + (F⃗ - 𝐱₄) × w⃗
# #         #   Let A = (D⃗ × w⃗), B = (E⃗ × w⃗), C = (F⃗ - 𝐱₄) × w⃗
# #         #   0 = Ar² + Br + C
# #         #   r = (-B - √(B²-4AC))/2A, -B + √(B²-4AC))/2A)
# #         #   s = ((q(r) - p₄)⋅w⃗/(w⃗ ⋅ w⃗)
# #         #   r is invalid if:
# #         #     1) A = 0
# #         #     2) B² < 4AC
# #         #     3) r < 0 or 1 < r   (Curve intersects, segment doesn't)
# #         #   s is invalid if:
# #         #     1) s < 0 or 1 < s   (Line intersects, segment doesn't)
# #         # If D⃗ × w⃗ = 0, there is only one intersection and the equation reduces to line
# #         # intersection.
# #         npoints = 0x00000000
# #         p₁ = Point()
# #         p₂ = Point()
# #         D⃗ = 2(q[1] +  q[2] - 2q[3])
# #         E⃗ =  4q[3] - 3q[1] -  q[2]
# #         w⃗ = l[2] - l[1]
# #         A = D⃗ × w⃗
# #         B = E⃗ × w⃗
# #         C = (q[1] - l[1]) × w⃗
# #         w = w⃗ ⋅ w⃗
# #         if abs(A) < 1e-8
# #             # Line intersection
# #             # Can B = 0 if A = 0 for non-trivial 𝐱?
# #             r = -C/B
# #             (-ϵ ≤ r ≤ 1 + ϵ) || return 0x00000000, SVector(p₁, p₂)
# #             p₁ = q(r)
# #             s = (p₁ - l[1]) ⋅ w⃗/w
# #             if (-ϵ ≤ s ≤ 1 + ϵ)
# #                 npoints = 0x00000001
# #             end
# #         elseif B^2 ≥ 4A*C
# #             # Quadratic intersection
# #             # The compiler seem seems to catch the √(B^2 - 4A*C), for common subexpression
# #             # elimination, so leaving for readability
# #             r₁ = (-B - √(B^2 - 4A*C))/2A
# #             r₂ = (-B + √(B^2 - 4A*C))/2A
# #             if (-ϵ ≤ r₁ ≤ 1 + ϵ)
# #                 p = q(r₁)
# #                 s₁ = (p - l[1]) ⋅ w⃗/w
# #                 if (-ϵ ≤ s₁ ≤ 1 + ϵ)
# #                     p₁ = p
# #                     npoints += 0x00000001
# #                 end
# #             end
# #             if (-ϵ ≤ r₂ ≤ 1 + ϵ)
# #                 p = q(r₂)
# #                 s₂ = (p - l[1]) ⋅ w⃗/w
# #                 if (-ϵ ≤ s₂ ≤ 1 + ϵ)
# #                     p₂ = p
# #                     npoints += 0x00000001
# #                 end
# #             end
# #             if npoints === 0x00000001 && p₁ === Point()
# #                 p₁ = p₂
# #             end
# #         end
# #         return npoints, SVector(p₁, p₂)
# #     end
# # end
#
# # Return the Laplacian of q, evalutated at r
# function laplacian(q::QuadraticSegment, r::Real)
#     return 4(q[1] + q[2] - 2q[3])
# end
#
# nearest_point(p::Point, q::QuadraticSegment) = nearest_point(p, q, 15)
# # Return the closest point on the curve to point p and the value of r such that q(r) = p_nearest
# # Uses at most N iterations of Newton-Raphson
# function nearest_point(p::Point, q::QuadraticSegment, N::Int64)
#     r = 0.5
#     Δr = 0.0
#     for i = 1:N
#         err = p - q(r)
#         grad = ∇(q, r)
#         if abs(grad[1]) > abs(grad[2])
#             Δr = err[1]/grad[1]
#         else
#             Δr = err[2]/grad[2]
#         end
#         r += Δr
#         if abs(Δr) < 1e-7
#             break
#         end
#     end
#     return r, q(r)
# end
#
# #function newton_raphson(f::Function, J⁻¹::Function, x)
# #    xₙ₊₁ = xₙ - J⁻¹(xₙ)f
# #end
#
#
#
#
#
#
#
#
#
#
#
#
# # Plot
# # -------------------------------------------------------------------------------------------------
# if enable_visualization
#     function convert_arguments(LS::Type{<:LineSegments}, q::QuadraticSegment)
#         rr = LinRange(0, 1, 15)
#         points = q.(rr)
#         coords = reduce(vcat, [[points[i], points[i+1]] for i = 1:length(points)-1])
#         return convert_arguments(LS, coords)
#     end
#
#     function convert_arguments(LS::Type{<:LineSegments}, Q::Vector{<:QuadraticSegment})
#         point_sets = [convert_arguments(LS, q) for q in Q]
#         return convert_arguments(LS, reduce(vcat, [pset[1] for pset in point_sets]))
#     end
# end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
"""
    FlipSignBridge{T, S1, S2, F, G}

Bridge a `G`-in-`S1` constraint into an `F`-in-`S2` constraint by multiplying
the function by `-1` and taking the point reflection of the set across the
origin. The flipped `F`-in-`S` constraint is stored in the `constraint`
field by convention.
"""
abstract type FlipSignBridge{
    T, S1<:MOI.AbstractSet, S2<:MOI.AbstractSet,
    F<:MOI.AbstractFunction, G<:MOI.AbstractFunction} <: SetMapBridge{T, S2, S1, F, G} end

map_function(::Type{<:FlipSignBridge{T}}, func) where {T} = MOIU.operate(-, T, func)
# The map is an involution
inverse_map_function(BT::Type{<:FlipSignBridge}, func) = map_function(BT, func)
# The map is symmetric
adjoint_map_function(BT::Type{<:FlipSignBridge}, func) = map_function(BT, func)
# The map is a symmetric involution
inverse_adjoint_map_function(BT::Type{<:FlipSignBridge}, func) = map_function(BT, func)

function MOI.delete(model::MOI.ModelLike, bridge::FlipSignBridge, i::IndexInVector)
    func = MOI.get(model, MOI.ConstraintFunction(), bridge.constraint)
    idx = setdiff(1:MOI.output_dimension(func), i.value)
    new_func = MOIU.eachscalar(func)[idx]
    set = MOI.get(model, MOI.ConstraintSet(), bridge.constraint)
    new_set = MOI.update_dimension(set, MOI.dimension(set) - 1)
    MOI.delete(model, bridge.constraint)
    bridge.constraint = MOI.add_constraint(model, new_func, new_set)
end

function MOI.modify(model::MOI.ModelLike, bridge::FlipSignBridge,
                    change::MOI.ScalarCoefficientChange)
    MOI.modify(
        model, bridge.constraint,
        MOI.ScalarCoefficientChange(change.variable, -change.new_coefficient))
end
function MOI.modify(model::MOI.ModelLike, bridge::FlipSignBridge,
                    change::MOI.MultirowChange{T}) where T
    new_coefficients = Tuple{Int64, T}[
        (index, -coef) for (index, coef) in change.new_coefficients]
    MOI.modify(model, bridge.constraint,
               MOI.MultirowChange(change.variable,
                                  new_coefficients))
end

"""
    GreaterToLessBridge{T, F<:MOI.AbstractScalarFunction, G<:MOI.AbstractScalarFunction} <:
        FlipSignBridge{T, MOI.GreaterThan{T}, MOI.LessThan{T}, F, G}

Transforms a `G`-in-`GreaterThan{T}` constraint into an `F`-in-`LessThan{T}`
constraint.
"""
struct GreaterToLessBridge{T, F<:MOI.AbstractScalarFunction, G<:MOI.AbstractScalarFunction} <:
    FlipSignBridge{T, MOI.GreaterThan{T}, MOI.LessThan{T}, F, G}
    constraint::CI{F, MOI.LessThan{T}}
end
map_set(::Type{<:GreaterToLessBridge}, set::MOI.GreaterThan) = MOI.LessThan(-set.lower)
inverse_map_set(::Type{<:GreaterToLessBridge}, set::MOI.LessThan) = MOI.GreaterThan(-set.upper)
function concrete_bridge_type(::Type{<:GreaterToLessBridge{T}},
                              G::Type{<:MOI.AbstractScalarFunction},
                              ::Type{MOI.GreaterThan{T}}) where T
    F = MOIU.promote_operation(-, T, G)
    return GreaterToLessBridge{T, F, G}
end

"""
    LessToGreaterBridge{T, F<:MOI.AbstractScalarFunction, G<:MOI.AbstractScalarFunction} <:
        FlipSignBridge{T, MOI.LessThan{T}, MOI.GreaterThan{T}, F, G}

Transforms a `G`-in-`LessThan{T}` constraint into an `F`-in-`GreaterThan{T}`
constraint.
"""
struct LessToGreaterBridge{T, F<:MOI.AbstractScalarFunction, G<:MOI.AbstractScalarFunction} <:
    FlipSignBridge{T, MOI.LessThan{T}, MOI.GreaterThan{T}, F, G}
    constraint::CI{F, MOI.GreaterThan{T}}
end
map_set(::Type{<:LessToGreaterBridge}, set::MOI.LessThan) = MOI.GreaterThan(-set.upper)
inverse_map_set(::Type{<:LessToGreaterBridge}, set::MOI.GreaterThan) = MOI.LessThan(-set.lower)
function concrete_bridge_type(::Type{<:LessToGreaterBridge{T}},
                              G::Type{<:MOI.AbstractScalarFunction},
                              ::Type{MOI.LessThan{T}}) where T
    F = MOIU.promote_operation(-, T, G)
    return LessToGreaterBridge{T, F, G}
end

"""
    NonnegToNonposBridge{T, F<:MOI.AbstractVectorFunction, G<:MOI.AbstractVectorFunction} <:
        FlipSignBridge{T, MOI.Nonnegatives, MOI.Nonpositives, F, G}

Transforms a `G`-in-`Nonnegatives` constraint into a `F`-in-`Nonpositives`
constraint.
"""
mutable struct NonnegToNonposBridge{T, F<:MOI.AbstractVectorFunction, G<:MOI.AbstractVectorFunction} <:
    FlipSignBridge{T, MOI.Nonnegatives, MOI.Nonpositives, F, G}
    constraint::CI{F, MOI.Nonpositives}
end
map_set(::Type{<:NonnegToNonposBridge}, set::MOI.Nonnegatives) = MOI.Nonpositives(set.dimension)
inverse_map_set(::Type{<:NonnegToNonposBridge}, set::MOI.Nonpositives) = MOI.Nonnegatives(set.dimension)
function concrete_bridge_type(::Type{<:NonnegToNonposBridge{T}},
                              G::Type{<:MOI.AbstractVectorFunction},
                              ::Type{MOI.Nonnegatives}) where T
    F = MOIU.promote_operation(-, T, G)
    return NonnegToNonposBridge{T, F, G}
end

"""
    NonposToNonnegBridge{T, F<:MOI.AbstractVectorFunction, G<:MOI.AbstractVectorFunction} <:
        FlipSignBridge{T, MOI.Nonpositives, MOI.Nonnegatives, F, G}

Transforms a `G`-in-`Nonpositives` constraint into a `F`-in-`Nonnegatives`
constraint.
"""
mutable struct NonposToNonnegBridge{T, F<:MOI.AbstractVectorFunction, G<:MOI.AbstractVectorFunction} <:
    FlipSignBridge{T, MOI.Nonpositives, MOI.Nonnegatives, F, G}
    constraint::CI{F, MOI.Nonnegatives}
end
map_set(::Type{<:NonposToNonnegBridge}, set::MOI.Nonpositives) = MOI.Nonnegatives(set.dimension)
inverse_map_set(::Type{<:NonposToNonnegBridge}, set::MOI.Nonnegatives) = MOI.Nonpositives(set.dimension)
function concrete_bridge_type(::Type{<:NonposToNonnegBridge{T}},
                              G::Type{<:MOI.AbstractVectorFunction},
                              ::Type{MOI.Nonpositives}) where T
    F = MOIU.promote_operation(-, T, G)
    return NonposToNonnegBridge{T, F, G}
end


text:
module test_05_ArraySignalsWithUnit

using ModiaPlot

include("test_base/test_05_ArraySignalsWithUnit.jl")

end


text:
using Test, LinearAlgebra, StaticArrays

using DomainSets
using CompositeTypes.Indexing

include("test_common.jl")
include("test_maps.jl")
include("test_generic_domain.jl")
include("test_specific_domains.jl")
include("test_canonical.jl")
include("test_setoperations.jl")
include("test_applications.jl")
include("test_readme.jl")


text:
include("zygote_patch.jl")

import YaoExtensions, Random

c = YaoExtensions.variational_circuit(5)
h = YaoExtensions.heisenberg(5)

function loss(h, c, θ) where N
    # the assign is nessesary!
    c = dispatch!(c, fill(θ, nparameters(c)))
    reg = apply!(zero_state(nqubits(c)), c)
    real(expect(h, reg))
end

reg0 = zero_state(5)
zygote_grad = Zygote.gradient(θ->loss(h, c, θ), 0.5)[1]

# check gradients
using Test
dispatch!(c, fill(0.5, nparameters(c)))
greg, gparams = expect'(h, zero_state(5)=>c)
true_grad = sum(gparams)

@test true_grad ≈ true_grad


text:
# Rubrene Ordejon
# Rubrene electron-phonon matrix elements

# Ordejón, P., Boskovic, D., Panhans, M., Ortmann, F., 2017. Ab initio study of
# electron-phonon coupling in rubrene. Phys. Rev. B 96, 035202.
# https://doi.org/10.1103/PhysRevB.96.035202

# Table II

# Units: cm^-1 meV
EffectiveH=(1208.9, 106.8)
EffectiveP=(117.9, 21.9)

EffectiveHGirlando=(1277, 99)
EffectivePGirlando=(77,   20)

# Table IV
# ωp      ωpg0  ωpgi
# cm^-1    meV   meV
PerModeOrdejon=
[57.8    -1.7   0.85
 59.6     1.4  -0.83
 89.0     1.6  -4.8
 107.3   -0.14  2.8
 139.1   -2.3  -3.7
 639.1   -7.5   1.0
 1011.2  -3.6  -0.04
 1344.7  19.8   0.04
 1593.3 -42.0  -0.12
]

PerModeGirlando =
[37.4    -0.9   3.4
 66.6     1.6  -6.6
 86.7    -0.6  -9.3
 106.3    0    -4.4
 125.1    1.4  -4.7
 631.2  -10.8   1.3
 1002.3  24.6   0
 1348.6  49.9   0
 1593.8 -45.6   1.6
]


text:
# This file was generated, do not modify it. # hide
X2 = hcat(X, X.LStat .* X.Age);


text:
abstract type FindBestAlg end
struct BruteForce <: FindBestAlg end
struct MaybeFaster <: FindBestAlg end

best_acronym(af::AcronymFamily, params::ScoreParams) = best_acronym(BruteForce(), af, params)

best_acronym(::BruteForce, af::AcronymFamily, params::ScoreParams) = maximum_by(last, (a, score(a, params)) for a in acronyms(af))

best_acronym(::MaybeFaster, af::AcronymFamily, params::ScoreParams) = maximum_by(last, _best_faster(af, trace, params) for trace in _family_traces(af))

function _best_faster(af::AcronymFamily, trace::AbstractVector{NTuple{2, Int}}, params::ScoreParams)
    ac_words = Union{String, Nothing}[]
    s = 0.
    group_seen = falses(length(params.word_groups))

    for (i, n) in trace
        words = af.components[n, i]

        if isempty(words)
            push!(ac_words, nothing)
            s -= params.gap_penalty

        else
            best_word = nothing
            best_score = -Inf

            for word in words
                g = params.word_to_group[word]
                ws = params.group_scores[g]

                if g in params.required_groups && !group_seen[g]
                    best_word = word
                    best_score = 0.
                    break
                elseif ws > best_score
                    best_word = word
                    best_score = ws
                end
            end

            push!(ac_words, best_word)
            s += best_score
            group_seen[params.word_to_group[best_word]] = true
        end
    end

    all(group_seen[g] for g in params.required_groups) || (s = -Inf)

    a = Acronym(af.word, ac_words, last.(trace))
    return a, s
end


text:
##### Beginning of file

module GitHubDstHost # Begin submodule RepoSnapshots.Hosts.GitHubDstHost

__precompile__(true)

import ..Types
import ..Utils

import Dates
import GitHub
import TimeZones

import ..delayederror

function new_github_session(
        ;
        github_organization::AbstractString,
        github_bot_username::AbstractString,
        github_bot_personal_access_token::AbstractString,
        )::Function

    _github_organization::String = strip(
        convert(String, github_organization)
        )
    _provided_github_bot_username::String = strip(
        convert(String, github_bot_username,)
        )
    _github_bot_personal_access_token::String = strip(
        convert(String, github_bot_personal_access_token)
        )
    function _get_github_username(auth::GitHub.Authorization)::String
        user_information::AbstractDict = GitHub.gh_get_json(
            GitHub.DEFAULT_API,
            "/user";
            auth = auth,
            )
        username::String = user_information["name"]
        username_stripped::String = strip(username)
        return username_stripped
    end

    @info("Attempting to authenticate to GitHub...")
    auth::GitHub.Authorization = GitHub.authenticate(
        _github_bot_personal_access_token
        )
    _github_username::String = _get_github_username(auth)
    if lowercase(strip(_github_username)) !=
            lowercase(strip(_provided_github_bot_username))
        delayederror(
            string(
                "Provided GitHub username ",
                "(\"$(_provided_github_bot_username)\") ",
                "does not match ",
                "actual GitHub username ",
                "(\"$(_github_username)\").",
                )
            )
    else
        @info(
            string(
                "Provided GitHub username matches ",
                "actual GitHub username.",
                ),
            _provided_github_bot_username,
            _github_username,
            )
    end
    @info("Successfully authenticated to GitHub :)")

    @info(
        string(
            "GitHub username: ",
            "$(_get_github_username(auth))",
            )
        )
    @info(
        string(
            "GitHub organization: ",
            "$(_github_organization)",
            )
        )

    repository_owner = GitHub.owner(
        _github_organization,
        true;
        auth = auth,
        )

    function _create_gist(params::AbstractDict)::Nothing
        gist_description::String = strip(params[:gist_description])
        gist_content::String = strip(params[:gist_content])
        @info("Attempting to create gist on GitHub...")
        create_gist_function = () ->
            GitHub.create_gist(
                ;
                auth = auth,
                params = Dict(
                    :public => true,
                    :description => gist_description,
                    :files => Dict(
                        "list.txt" => Dict("content" => gist_content,),
                        ),
                    ),
                )
        Utils.retry_function_until_success(
            create_gist_function;
            max_attempts = 10,
            seconds_to_wait_between_attempts = 180,
            )
        @info("Successfully created gist on GitHub.")
        return nothing
    end

    function _get_all_gists()::Vector{GitHub.Gist}
        @info("Loading the list of all of my GitHub gists")
        full_gist_list::Vector{GitHub.Gist} = GitHub.Gist[]
        need_to_continue::Bool = true
        current_page_number::Int = 1
        while need_to_continue
            gists, page_data = GitHub.gists(
                _github_username;
                params = Dict(
                    "per_page" => 100,
                    "page" => current_page_number,
                    ),
                auth = auth,
                )
            if length(gists) == 0
                need_to_continue = false
            else
                for x in gists
                    if x in full_gist_list
                    else
                        push!(full_gist_list, x)
                    end
                end
                need_to_continue = true
                current_page_number += 1
            end
        end
        unique_gist_list::Vector{GitHub.Gist} = unique(full_gist_list)
        return unique_gist_list
    end

    function _retrieve_gist(params::AbstractDict)::String
        gist_description_to_match::String = params[:gist_description]
        correct_gist_id::String = ""
        all_my_gists = _get_all_gists()
        for gist in all_my_gists
            if gist.description == gist_description_to_match
                correct_gist_id = gist.id
            end
        end
        result::String = ""
        if length(correct_gist_id) > 0
            @info("Downloading the correct GitHub gist")
            correct_gist::GitHub.Gist = GitHub.gist(
                correct_gist_id;
                auth = auth,
                )
            correct_gist_content::String = correct_gist.files[
                "list.txt"]["content"]
            result = correct_gist_content
        else
            result = ""
        end
        if length(result) == 0
            delayederror("Could not find the matching Gist")
        end
        return result
    end

    function _delete_gists(params::AbstractDict)::Nothing
        gist_description_to_match::String = params[:gist_description]
        list_of_gist_ids_to_delete::Vector{String} = String[]
        all_my_gists::Vector{GitHub.Gist} = _get_all_gists()
        for gist in all_my_gists
            if gist.description == gist_description_to_match
                push!(list_of_gist_ids_to_delete, strip(gist.id),)
            end
        end
        for gist_id_to_delete in list_of_gist_ids_to_delete
            GitHub.delete_gist(gist_id_to_delete;auth = auth,)
            @info(string("Deleted GitHub gist id $(gist_id_to_delete)"))
        end
        return nothing
    end

    function _delete_gists_older_than_minutes(params::AbstractDict)::Nothing
        time::TimeZones.ZonedDateTime =
            params[:time]
        delete_gists_older_than_minutes::Int =
            params[:delete_gists_older_than_minutes]
        max_gist_age_milliseconds::Int =
            delete_gists_older_than_minutes*60*1000
        list_of_gist_ids_to_delete::Vector{String} = String[]
        all_my_gists::Vector{GitHub.Gist} = _get_all_gists()
        for gist in all_my_gists
            gist_updated_at = gist.updated_at
            gist_updated_at_zoned = TimeZones.ZonedDateTime(
                gist_updated_at,
                TimeZones.localzone(),
                )
            gist_age = time - gist_updated_at_zoned
            if gist_age.value > max_gist_age_milliseconds
                push!(list_of_gist_ids_to_delete, strip(gist.id),)
            end
        end
        for gist_id_to_delete in list_of_gist_ids_to_delete
            GitHub.delete_gist(gist_id_to_delete;auth = auth,)
            @info(string("Deleted GitHub gist id $(gist_id_to_delete)"))
        end
        return nothing
    end

    function _repo_name_with_org(
            ;
            repo::AbstractString,
            org::AbstractString,
            )::String
        repo_name_without_org::String = _repo_name_without_org(
            ;
            repo = repo,
            org = org,
            )
        org_stripped::String = strip(
            strip(strip(strip(strip(convert(String, org)), '/')), '/')
            )
        result::String = string(
            org_stripped,
            "/",
            repo_name_without_org,
            )
        return result
    end

    function _repo_name_without_org(
            ;
            repo::AbstractString,
            org::AbstractString,
            )::String
        repo_stripped::String = strip(
            strip(strip(strip(strip(convert(String, repo)), '/')), '/')
            )
        org_stripped::String = strip(
            strip(strip(strip(strip(convert(String, org)), '/')), '/')
            )
        if length(org_stripped) == 0
            result = repo_stripped
        else
            repo_stripped_lowercase::String = lowercase(repo_stripped)
            org_stripped_lowercase::String = lowercase(org_stripped)
            org_stripped_lowercase_withtrailingslash::String = string(
                org_stripped_lowercase,
                "/",
                )
            if startswith(repo_stripped_lowercase,
                    org_stripped_lowercase_withtrailingslash)
                index_start =
                    length(org_stripped_lowercase_withtrailingslash) + 1
                result = repo_stripped[index_start:end]
            else
                result = repo_stripped
            end
        end
        return result
    end

    function _get_destination_url(
            ;
            repo_name::String,
            credentials::Symbol,
            )::String
        repo_name_without_org::String = _repo_name_without_org(
            ;
            repo = repo_name,
            org = _github_organization,
            )
        result::String = ""
        if credentials == :with_auth
            result = string(
                "https://",
                _github_username,
                ":",
                _github_bot_personal_access_token,
                "@",
                "github.com/",
                _github_organization,
                "/",
                repo_name_without_org,
                )
        elseif credentials == :with_redacted_auth
            result = string(
                "https://",
                _github_username,
                ":",
                "*****",
                "@",
                "github.com/",
                _github_organization,
                "/",
                repo_name_without_org,
                )
        elseif credentials == :without_auth
            result =string(
                "https://",
                "github.com/",
                _github_organization,
                "/",
                repo_name_without_org,
                )
        else
            delayederror("$(credentials) is not a supported value for credentials")
        end
        return result
    end

    function _github_repo_exists(
            ;
            repo_name::String,
            )::Bool
        repo_name_with_org = _repo_name_with_org(
            ;
            repo = repo_name,
            org = _github_organization,
            )
        result::Bool = try
            repo = GitHub.repo(
                repo_name_with_org;
                auth = auth,
                )
            true
        catch
            false
        end
        return result
    end

    function _create_repo(params::AbstractDict)::Nothing
        repo_name::String = params[:repo_name]
        repo_name_with_org::String = _repo_name_with_org(
            ;
            repo = repo_name,
            org = _github_organization,
            )
        repo_name_without_org::String = _repo_name_without_org(
            ;
            repo = repo_name,
            org = _github_organization,
            )
        repo_destination_url_without_auth = _get_destination_url(
            ;
            repo_name = repo_name_without_org,
            credentials = :without_auth,
            )
        if Utils._url_exists(repo_destination_url_without_auth)
            @info("According to HTTP GET request, the repo exists.")
        else
            if _github_repo_exists(; repo_name = repo_name_with_org)
                @info("According to the GitHub API, the repo exists.")
            else
                @info(
                    string("Creating new repo on GitHub"),
                    repo_destination_url_without_auth,
                    )
                repo = GitHub.create_repo(
                    repository_owner,
                    repo_name_without_org,
                    Dict{String, Any}(
                        "has_issues" => "false",
                        "has_wiki" => "false",
                        );
                    auth = auth,
                    )
            end
        end
        return nothing
    end

    function _generate_new_repo_description(
            params::AbstractDict,
            )::String
        source_url::String = params[:source_url]
        when::TimeZones.ZonedDateTime = params[:when]
        time_zone::TimeZones.TimeZone = params[:time_zone]
        by::String = strip(string("@", _github_username))

        new_description::String = Utils.default_repo_description(
            ;
            from = source_url,
            when = when,
            time_zone = time_zone,
            by = by,
            )

        return new_description
    end

    function _update_repo_description(params::AbstractDict)::Nothing
        repo_name::String = params[:repo_name]
        new_repo_description = params[:new_repo_description]
        _create_repo(
            Dict(
                :repo_name => repo_name,
                ),
            )
        repo_name_with_org::String = _repo_name_with_org(
            ;
            repo = repo_name,
            org = _github_organization,
            )
        github_repo_function = () -> GitHub.repo(
            repo_name_with_org;
            auth = auth,
            )
        repo = Utils.retry_function_until_success(
            github_repo_function;
            max_attempts = 10,
            seconds_to_wait_between_attempts = 180,
            )
        @info("Attempting to update repo description on GitHub...")
        github_update_description_function = () ->
            GitHub.gh_patch_json(
                GitHub.DEFAULT_API,
                "/repos/$(GitHub.name(repo.owner))/$(GitHub.name(repo.name))";
                auth = auth,
                params = Dict(
                    "name" => GitHub.name(repo.name),
                    "description" => new_repo_description,
                    ),
                )
        result = Utils.retry_function_until_success(
            github_update_description_function;
            max_attempts = 10,
            seconds_to_wait_between_attempts = 180,
            )
        @info("Successfully updated repo description on GitHub")
        return nothing
    end

    function _github_provider(task::Symbol)::Function
        if task == :create_gist
            return _create_gist
        elseif task == :retrieve_gist
            return _retrieve_gist
        elseif task == :delete_gists
            return _delete_gists
        elseif task == :create_repo
            return _create_repo
        elseif task == :generate_new_repo_description
            return _generate_new_repo_description
        elseif task == :update_repo_description
            return _update_repo_description
        elseif task == :delete_gists_older_than_minutes
            return _delete_gists_older_than_minutes
        elseif task == :get_destination_url
            return _get_destination_url
        else
            delayederror("$(task) is not a valid task")
        end
    end

    return _github_provider
end

end # End submodule RepoSnapshots.Hosts.GitHubDstHost

##### End of file


text:
module Interactive
using  ..ShellInterop
using  ..MacroExpander
using  ..MacroDatabase
using  ..Tokenize

mutable struct InteractConfig
  expand::Bool
  tokens::Bool
  macrocall::Bool
  quit::Bool
  stacktrace::StackTrace

  InteractConfig() = new(false,false,false,false,[])
end

mutable struct InputState
  top::Bool
  topmacro::Bool
  incomplete_string::Bool
  incomplete_parens::Bool
  escaped_newline::Bool
  nested_macros::Int
  parens_balance::Vector{String}
  quotes_incomplete::Symbol
  found_macro::Bool

  InputState() = new(true,true,false,false,false,0,[],:none,false)
end

macro with(change, actions)
  changed  = change.args[1]
  newvalue = change.args[2]
  quote
    oldvalue = $(esc(changed))
    $(esc(changed)) = $newvalue
    $(esc(actions))
    $(esc(changed)) = oldvalue
  end
end

macro toggle(x)
  esc(:($x=!$x))
end

function main_loop(shell)
  config = InteractConfig()

  while !config.quit
    try
      input = obtain_input(config, InputState())
      local output
      if !config.quit
        if config.expand
          output = input
        else
          output = read_sh(shell, input)
        end
        if output != ""
          pre="| "
          output = "|/\n$pre"*join(split(output,"\n"),"\n$pre")*"\n|\\"
          println(output)
        end
      end
    catch err
      print_with_color(:red, sprint(showerror, err), "\n")
      config.stacktrace = stacktrace()
      println()
    end
  end
end

function obtain_input(config::InteractConfig,is::InputState)::String
  print_prompt(config, is)
  input = readline()

  process_input(input, config, is)
end

function process_input(input, config::InteractConfig, is::InputState)::String
  next_line() = obtain_input(config,is)

  if is.incomplete_string
    tokens = tokenize(input, incomplete=is.quotes_incomplete)
    if isempty(tokens)
      return "\n$(next_line())"
    elseif length(tokens) == 1
      t = tokens[1]
      is.incomplete_string = t.incomplete
      is.quotes_incomplete = t.incomplete ? t.typ : :none
      return "$input\n$(t.incomplete ? next_line() : "")"
    else
      is.incomplete_string = false
      is.quotes_incomplete = :none

      t = tokens[1]
      i = t.index+length(t)+2

      incput = input[1:i-1]
      comput = input[i:end]

      return "$incput$(process_input(comput,config,is))"
    end

  elseif is.top && is_magic_command(input)
    exec_magic(input, config)
  else
    ma = macro_analysis = analyze_line(input)
    print_extra_info(config, input, macro_analysis)

    if ismacro(ma) && ismultiline(ma)
      is.nested_macros += 1
      is.found_macro = true
    elseif isend(ma) && is.nested_macros > 0
      is.nested_macros -= 1
    end

    if is.nested_macros > 0
       @with is.topmacro = false begin
         result = "$input\n$(next_line())"
       end

       return is.topmacro ? compile_str(result) : result
    else
      if is.found_macro
        if isend(macro_analysis)
          return input
        else
          return "$input\n$(next_line())"
        end

      else
        if ismacro(macro_analysis) # single line macro
          return compile_str(input)

        elseif endswith(input, "\\")
          is.escaped_newline = true
          return "$input\n$(next_line())"

        else
          check_balance(input, is)
          if is.incomplete_parens || is.incomplete_string
            return "$input\n$(next_line())"
          else
            return input
          end
        end
      end
    end
  end
end

cute_invert_lparen(p) = (s="-()-{}-[]-")[search(s,p)+1] # swoooon
cute_invert_rparen(p) = (s="-()-{}-[]-")[search(s,p)-1] # this code is so cute !!

invert_parens(x) =
  contains("([{", x) ? cute_invert_lparen(x) : cute_invert_rparen(x)

isparens(t::Token) = t.typ == :special && contains("(){}[]", t.token)

function check_balance(input::String, is::InputState)
  input_tokens = tokenize(input, incomplete=is.quotes_incomplete)
  if !isempty(input_tokens) && last(input_tokens).incomplete
    is.quotes_incomplete = last(input_tokens).typ
  end

  for t in filter(isparens, input_tokens)
    current_parens = t.token

    if contains("([{", current_parens)
      push!(is.parens_balance, current_parens)

    else
      if isempty(is.parens_balance)
        error("Unmatched parens '$p'")

      elseif last(is.parens_balance) == invert_parens(current_parens)
        pop!(is.parens_balance)

      else
        error("Wrong matching parens: '$(last(is.parens_balance))$current_parens'")
      end
    end
  end
  is.incomplete_parens = !isempty(is.parens_balance)

  if !isempty(input_tokens)
    lt = last(input_tokens)
    if lt.incomplete
      is.incomplete_string = true
      is.quotes_incomplete = lt.typ
    end
  end
end

function print_prompt(config::InteractConfig, is::InputState)
  print_with_color(
    prompt_color(config),
    string(prompt_container_thingy(config), prompt_arrow(is)," "),
    bold=true
  )
end

function prompt_container_thingy(config::InteractConfig)
  e = config.expand    ? 'e' : '-'
  t = config.tokens    ? 't' : '-'
  m = config.macrocall ? 'm' : '-'

  return "[$e$t$m]"
end

function prompt_arrow(is::InputState)
  is.incomplete_string  ? " '>" :
  is.incomplete_parens  ? "()>" :
  is.escaped_newline    ? "//>" :
  is.nested_macros != 0 ? " >>" :
  true && true || true  ? "==>" : ""
end

function prompt_color(config::InteractConfig)
  config.expand ? :green : :red
end

function is_magic_command(input)
  startswith(input, ".")
end

function exec_magic(input, config::InteractConfig)::String
  iscommand(cmd) = input in (".$cmd", ".$(cmd[1])")

  if iscommand("quit")
    config.quit = true

  elseif iscommand("expand")
    @toggle config.expand

  elseif iscommand("call")
    @toggle config.macrocall

  elseif iscommand("token")
    @toggle config.tokens

  elseif iscommand("list")
    for (m,f) in Macros
      println("-- $m")
    end

  elseif iscommand("stacktrace")
    display(config.stacktrace)
    println()

  elseif iscommand("help")
    println("""
      .q | .quit       == end session
      .e | .expand     == toggles if macros are expanded or evaluated
      .t | .token      == shows tokens of input line
      .c | .call       == shows macrocall analysis
      .l | .list       == list macros
      .s | .stacktrace == show last stacktrace from a julia error
      .h | .help       == this help text
    """)

  else
    print_with_color(:red, "invalid fxcrd_i command $input. See '.help'\n")

  end
  return ""
end

function print_extra_info(config,input,mm)
  if config.tokens
    print_with_color(:yellow, "[t~~] ")
    for t in tokenize(input)
      print("('$(t.token)'::$(t.typ)) ")
    end
    println()
  end
  if !ismultiline(mm) && ismacro(mm) && config.macrocall
    function pm(s)
      print_with_color(:red, "[~c~]")
      println(s)
    end

    pm("name=$(mm.name)")
    pm("args=$(join(mm.args, " "))")
    pm("input ----------")
    pm("----------------")
    println(mm.input)
    pm("---------------")
  end
end

end

max_stars_repo_name: fcard/Fxcrd

max_issues_repo_name: fcard/Fxcrd

max_forks_repo_name: fcard/Fxcrd


text:
using DiffEqUncertainty, Distributions
import ModelingToolkit: get_defaults, varmap_to_vars
#stab_error(sol,index) = 1e3*sqrt(mean(abs2,cumsum(sol[index,:])))
running_mean(sol,index) = 1e3*sqrt.(map(/,cumsum(cumsum(sol[index,:]).^2),1:length(sol)))
function stab_error(sol, index)
    θ = zero(eltype(sol))
    err2 = zero(eltype(sol))
    N=length(sol)
    @inbounds for ω in sol[index,:]
        θ+=ω
        err2+= abs2(θ)
    end
    1e3*sqrt(err2/N)
end

#utilities
indexof(sym::Num,syms) = findfirst(isequal(sym),syms)
indexof(sym::Symbol,syms) = findfirst(isequal(sym),[s.name for s in syms])

function loss(θ;full=true)
    Z1,Z2,P2,Z3,P3,Z4,K_dc = θ
    prob, p , sys = gimbal_conntroller(;full,
                                       Z1, Z2, P2, Z3, P3, Z4,K_dc)

    pnew = get_defaults(sys)
    pnew[p._k_s] = Uniform(0.2-0.04,0.2+0.04)
    pnew[p._k_v] = Uniform(0.545-0.109,0.545+0.109)
    pnew[p._ω_brk] = Uniform(0.074-0.0148,0.074+0.0148)
    pnew[p._T_c] = Uniform(0.106-0.0212,0.106+0.0212)
    pnew[p._T_brk] = Uniform(0.177-0.0354,0.177+0.0354)
    p_dist = varmap_to_vars(pnew,parameters(sys))

    index = indexof(p.ω,states(sys))
    g(x) = stab_error(x,index)
    expectation(g, prob, prob.u0, p_dist,
                Koopman(), Rodas5(); ireltol=1e-1)[1]
    #expectation(g, prob, prob.u0, prob.p, MonteCarlo(), Rodas5(); trajectories = 100)
end

export loss


text:
# MLJIteration

using MLJ

for control in MLJ.MLJIteration.CONTROLS
    eval(:($control))
end

IterationControl.skip(Step(2))
IterationControl.debug(Step(2))

IteratedModel

# MLJSerialization

Save()
@test MLJ.save isa Function

# MLJOpenML

@test OpenML.load isa Function

true