text:
"""
This is a time-aware distance. It contains an inner distance (currently only [`OptimalTransportRootDistance`](@ref) and [`KernelWassersteinRootDistance`](@ref)supported), and some parameters that are specific to the time dimension, example:
```
dist = TimeDistance(
    inner = OptimalTransportRootDistance(
        domain = Continuous(),
        p      = 2,
        weight = simplex_residueweight,
    ),
    tp = 2,
    c  = 0.1,
)
```
`tp` is the same as `p` but for the time dimension, and `c` trades off the distance along the time axis with the distance along the frequency axis. A smaller `c` makes it cheaper to transport mass across time. The frequency axis spans `[-π,π]` and the time axis is the non-negative integers, which should give you an idea for how to make this trade-off.
"""
TimeDistance
@kwdef struct TimeDistance{D,T} <: AbstractDistance
    inner::D
    tp::Int = 2
    c::T    = 0.1
end

"""
    TimeVaryingAR{T <: AbstractRoots} <: AbstractModel

This model represents a rational spectrogram, i.e., a rational spectrum that changes with time. See [`TimeWindow`](@ref) for the corresponding fit method and [`TimeDistance`](@ref) for a time-aware distance.

Internally, this model stores a vector of [`AR`](@ref).
"""
struct TimeVaryingAR{MT<:AR, VT <: AbstractVector{MT}} <: AbstractModel
    models::VT
end

@inline Base.length(m::TimeVaryingAR) = sum(length, m.models)-length(m.models) # the length of a model is 1 more than the number of poles due to the numerator
@inline Base.size(m::TimeVaryingAR) = (length(m.models[1])-1, length(m.models))
@inline Base.eachindex(m::TimeVaryingAR) = 1:length(m)

@inline function Base.getindex(m::TimeVaryingAR, i::Int, j::Int)
    m.models[j].pc[i]
end

@inline function divmod(a,b)
    c = a÷b
    d = a - c*b
    c,d
end

@inline function mi2ij(m,i)
    r,c = size(m)
    j, i = divmod((i-1), r)
    i+1, j+1
end

@inline function Base.getindex(m::TimeVaryingAR, i::Int)
    i,j = mi2ij(m,i)
    m[i,j]
end

PolynomialRoots.roots(::Continuous, m::TimeVaryingAR) = roots.(Continuous(), m.models)

function Base.isapprox(m1::TimeVaryingAR, m2::TimeVaryingAR, args...; kwargs...)
    all(isapprox(m1, m2, args...; kwargs...) for m1 in m1.models, m2 in m2.models)
end

change_precision(F, m::TimeVaryingAR) = TimeVaryingAR(change_precision.(F, m.models))

"""
We define a custom fit method for fitting time varying spectra, [`TimeWindow`](@ref). It takes as arguments an inner fitmethod, the number of points that form a time window, and the number of points that overlap between two consecutive time windows:
# Example
```
fitmethod = TimeWindow(TLS(na=2), 1000, 500)
y = sin.(0:0.1:100);
model = fitmethod(y)
```
"""
TimeWindow
@kwdef struct TimeWindow{T<:FitMethod} <: FitMethod
    inner::T
    n::Int
    noverlap::Int = n÷2
end

function fitmodel(fm::TimeWindow, x; showprogress = true)
    n,na,noverlap = fm.n, fm.inner.na, fm.noverlap
    if fm.inner.λ > 0
        AS = similar(x, n, na+1)
    else
        AS = similar(x, n-na, na+1)
    end
    if showprogress
        models = @showprogress 1 "Model estimation" map(arraysplit(x,n,noverlap)) do slice
            fitmodel!(AS, fm.inner, slice)
        end
    else
        models = map(arraysplit(x,n,noverlap)) do slice
            fitmodel!(AS, fm.inner, slice)
        end
    end
    TimeVaryingAR(models)
end

preprocess_roots(d, e::Vector{<:AbstractRoots}) = e

distmat_euclidean(m1::AbstractModel,m2::AbstractModel,p, tp, c) = distmat_euclidean!(zeros(length(m1),length(m2)), m1, m2, p, tp, c)

function distmat_euclidean!(D, m1::TimeVaryingAR, m2::TimeVaryingAR, p, tp, c)
    @assert size(D) == (length(m1), length(m2))
    _distmat_kernel!( D, m1, m2, c,
        p == 1 ? Base.FastMath.abs_fast :
        p == 2 ? Base.FastMath.abs2_fast : error("p must be 1 or 2"),
        tp == 1 ? Base.FastMath.abs_fast :
        tp == 2 ? Base.FastMath.abs2_fast : error("tp must be 1 or 2"),
    )
end

function _distmat_kernel!(D,m1,m2,c,f1::F1,f2::F2) where {F1,F2}
    @fastmath @inbounds for j in 1:length(m2)
        for i in 1:length(m1)
            _,t1 = mi2ij(m1,i)
            _,t2 = mi2ij(m2,j)
            D[i,j] = f1(m1[i]-m2[j]) + c*f2(t1-t2)
        end
    end
    D
end

function evaluate(od::TimeDistance{<:OptimalTransportRootDistance}, m1::TimeVaryingAR,m2::TimeVaryingAR; solver=sinkhorn_log!, kwargs...)
    d     = od.inner
    @assert d.domain isa Continuous "TimeDistance currently only works in continuous domain, open an issue with a motivation for why you require support for discrete domain and I might be able to add it."
    D     = distmat_euclidean(m1, m2, d.p, od.tp, od.c)
    w1    = s1(reduce(vcat,map(d.weight, m1.models)))
    w2    = s1(reduce(vcat,map(d.weight, m2.models)))
    C     = solver(D,w1,w2; β=d.β, kwargs...)[1]
    if any(isnan, C)
        @info("Nan in OptimalTransportRootDistance, increasing precision")
        C     = solver(big.(D),big.(w1),big.(w2); β=d.β, kwargs...)[1]
        any(isnan, C) && error("No solution found by solver $(solver), check your input and consider increasing β ($(d.β)).")
        eltype(D).(C)
    end
    dot(C, D)
end

@inline function Base.resize!(D::AbstractMatrix, s1::Int, s2::Int)
    l = s1 * s2
    Dv = resize!(vec(D), l)
    reshape(Dv, s1, s2)
end

function evaluate(
    od::TimeDistance{<:KernelWassersteinRootDistance},
    m1::TimeVaryingAR,
    m2::TimeVaryingAR;
    D = zeros(max(length(m1), length(m2)), max(length(m1), length(m2))),
    kwargs...,
)
    d = od.inner
    λ = d.λ
    @assert d.domain isa Continuous "TimeDistance currently only works in continuous domain, open an issue with a motivation for why you require support for discrete domain and I might be able to add it."

    D = resize!(D, length(m2), length(m2))
    distmat_euclidean!(D, m2, m2, 2, od.tp, od.c)
    @avx @. D = exp(-λ * D)
    c = mean(D)

    D = resize!(D, length(m1), length(m1))
    distmat_euclidean!(D, m1, m1, 2, od.tp, od.c)
    @avx @. D = exp(-λ * D)
    c += mean(D)

    D = resize!(D, length(m1), length(m2))
    distmat_euclidean!(D, m1, m2, 2, od.tp, od.c)
    @avx @. D = exp(-λ * D)
    c -= 2 * mean(D)

    c
end

"""
    distance_profile(od::TimeDistance, q::TimeVaryingAR, y::TimeVaryingAR; normalize_each_timestep = false, kwargs...)

Optimized method to compute the distance profile corresponding to sliding the short query `q` over the longer `y`.
"""
function SlidingDistancesBase.distance_profile(od::TimeDistance{<:OptimalTransportRootDistance}, q::TimeVaryingAR, y::TimeVaryingAR; normalize_each_timestep = false, kwargs...)
    d     = od.inner
    @assert d.domain isa Continuous "TimeDistance currently only works in continuous domain, open an issue with a motivation for why you require support for discrete domain and I might be able to add it."
    any(methods(d.weight).ms) do m
        m.nargs == 3
    end || throw(ArgumentError("distance_profile requires the weight function of the distance to support in-place update on the form `weightfun(weights, roots)`. See simplex_residueweight, residueweight, unitweight"))

    T  = eltype(q.models[1].a)
    N  = length(q)
    na = length(q.models[1].pc)
    nq = length(q.models)
    ny = length(y.models)
    w1 = s1(reduce(vcat,map(d.weight, q.models)))
    w2 = similar(w1)
    C  = Matrix{T}(undef, N, N)

    workspace = SinkhornLogWorkspace(T,N,N)

    @showprogress 1 "Distance profile" map(1:ny-nq) do i
        @views Y = TimeVaryingAR(y.models[i:i+nq-1])
        distmat_euclidean!(C, q, Y, d.p, od.tp, od.c)
        inds = (1:na)
        for i in eachindex(Y.models)
            d.weight(@view(w2[inds]), Y.models[i])
            if normalize_each_timestep
                w2[inds] ./= sum(@view(w2[inds]))
            end
            inds = inds .+ na
        end
        w2  ./= sum(w2)
        Γ = sinkhorn_log!(workspace,C,w1,w2; β=d.β, kwargs...)[1]
        dot(Γ,C)
    end
end

"""
    SlidingDistancesBase.distance_profile(d::ConvOptimalTransportDistance, q::AbstractMatrix, y::AbstractMatrix; stride=1, kwargs...)

Optimized method for [`ConvOptimalTransportDistance`](@ref). To get smooth distance profiles, a slightly higher β than for barycenters is recommended.  β around 0.01 should do fine.
- `stride`: allows you to reduce computations by setting `stride > 1`.

See also [`SlidingConvOptimalTransportDistance`](@ref)
"""

function SlidingDistancesBase.distance_profile(d::ConvOptimalTransportDistance, Q::AbstractMatrix, Y::AbstractMatrix; stride=1, kwargs...)

    T   = eltype(Q)
    m,n = size(Q)
    stride <= n || throw(ArgumentError("The stride can not be longer than the length of `q`"))
    N   = lastlength(Y)
    A   = copy(Q)
    B   = copy(getwindow(Y, n, 1))
    A ./= sum(A)
    workspace = SCWorkspace(A,B,d.β)
    U,V = workspace.U, workspace.V

    D = similar(Q, (N-n)÷stride+1)
    # Di = similar(Q, (N-n)÷stride+1)
    sB = sum(B) - sum(B[:,n-stride+1:n])
    iD = 0
    @views for i = 1:stride:N-n+1
        B .= getwindow(Y, n, i) # TODO: this is wasteful, only update one column and shift the rest
        sB += sum(B[:,n-stride+1:n])
        sB1 = sum(B[:,1:stride])
        @avx B ./= sB
        D[iD += 1] = sinkhorn_convolutional(workspace, A, B; β = d.β, initUV = false, kwargs...)[1]
        if d.invariant_axis != 0
            ia = d.invariant_axis
            sum_axis = ia == 1 ? 2 : 1
            v,u = vec(sum(A, dims=sum_axis)), vec(sum(B, dims=sum_axis))
            # v ./= sum(v); u ./= sum(u)
            invariant_cost = discrete_grid_transportcost(v, u; inplace=true)
            # Di[iD] = invariant_cost
            D[iD] -= invariant_cost
        end
        @avx U[:,1:end-stride] .= exp.(U[:,1+stride:end]) # Warm-start gives 3x imrovement in simple benchmark
        @avx V[:,1:end-stride] .= exp.(V[:,1+stride:end])
        @avx U[:,end-stride+1:end] .= 1#exp.(U[:,end]) # These make it slower
        @avx V[:,end-stride+1:end] .= 1#exp.(V[:,end])
        sB -= sB1
    end
    D#, Di
end

function SlidingDistancesBase.distance_profile(d::ConvOptimalTransportDistance, q::DSP.Periodograms.TFR, y::DSP.Periodograms.TFR; kwargs...)
    df  = d.dynamic_floor
    Q   = normalize_spectrogram(q, df)
    Y   = normalize_spectrogram(y, df)
    distance_profile(d, Q, Y; kwargs...)
end

# function distance_profile(d::AbstractDistance, q::TimeVaryingAR, y::TimeVaryingAR; normalize_each_timestep = false, kwargs...) where F
#     T  = eltype(q.models[1].a)
#     N  = length(q)
#     na = length(q.models[1].pc)
#     nq = length(q.models)
#     ny = length(y.models)
#     w1 = s1(reduce(vcat,map(d.weight, q.models)))
#     w2 = s1(reduce(vcat,map(d.weight, y.models)))
#     C  = Matrix{T}(undef, N, N)
#
#     any(methods(d.weight).ms) do m
#         m.nargs == 3
#     end || throw(ArgumentError("distance_profile requires the weight function of the distance to support in-place update on the form `weightfun(weights, roots)`. See simplex_residueweight, residueweight, unitweight"))
#
#     # workspace = SinkhornLogWorkspace(T,N,N)
#
#     @showprogress 1 "Distance profile"  map(1:ny-nq) do i
#         @views Y = TimeVaryingAR(y.models[i:i+nq-1])
#         inds = (1:na)
#         for i in eachindex(Y.models)
#             if normalize_each_timestep
#                 w2[inds] ./= sum(@view(w2[inds]))
#             end
#             inds = inds .+ na
#         end
#         w2  ./= sum(w2)
#         evaluate(d,e1,e2,w1,w2)
#     end
# end

struct PrimitiveModel{R <: AbstractVector,W <: AbstractVector} <: AbstractModel
    roots::R
    weights::W
end

function PrimitiveModel(d::AbstractDistance, m::AbstractModel)
    r = roots(domain(d), m)
    w = d.weight(r)
    PrimitiveModel(r.r,w)
end

function distmat(od::TimeDistance{<:OptimalTransportRootDistance}, models::Vector{<:TimeVaryingAR}; normalize=false, kwargs...)
    n = length(models)
    d = od.inner
    T = eltype(models[1].models[1].a)
    N = length(models[1])
    C = [Matrix{T}(undef, N, N) for _ = 1:nthreads()]

    weights = map(models) do q
        s1(reduce(vcat,map(d.weight, q.models)))
    end

    workspace = [SinkhornLogWorkspace(T, N, N) for _ = 1:nthreads()]
    dists = [deepcopy(d) for _ = 1:nthreads()]
    D = zeros(T, n, n)
    @showprogress 1 "distmat" for i = 1:n
        m1 = models[i]
        distmat_euclidean!(C[1], m1, m1, d.p, od.tp, od.c)
        Γ = sinkhorn_log!(workspace[1], C[1], weights[i], weights[i]; β = d.β, kwargs...)[1]
        D[i, i] = dot(Γ, C[1])
        Threads.@threads for j = i+1:n
            m2 = models[j]
            id = threadid()
            distmat_euclidean!(C[id], m1, m2, d.p, od.tp, od.c)
            Γ = sinkhorn_log!(
                workspace[id],
                C[id],
                weights[i],
                weights[j];
                β = d.β,
                kwargs...,
            )[1]
            D[i, j] = dot(Γ, C[id])
            (D[j, i] = D[i, j])
        end
    end
    normalize && symmetrize!(D)
    Symmetric(D)

end

function distmat(od::TimeDistance{<:KernelWassersteinRootDistance}, models::Vector{<:TimeVaryingAR}; normalize=false, kwargs...)
    n = length(models)
    d = od.inner
    λ = d.λ
    T = eltype(models[1].models[1].a)
    N = length(models[1])
    C = [Matrix{T}(undef, N, N) for _ = 1:nthreads()]
    dists = [deepcopy(d) for _ = 1:nthreads()]
    D = zeros(T, n, n)
    # error("need to take care of the size of C")
    @showprogress 1 "distmat" for i = 1:n
        m1 = models[i]
        Threads.@threads for j = i:n
            m2 = models[j]
            id = threadid()

            C[id] = resize!(C[id], length(m2), length(m2))
            distmat_euclidean!(C[id], m2, m2, 2, od.tp, od.c)
            @avx @. C[id] = exp(-λ * C[id])
            c = mean(C[id])

            C[id] = resize!(C[id], length(m1), length(m1))
            distmat_euclidean!(C[id], m1, m1, 2, od.tp, od.c)
            @avx @. C[id] = exp(-λ * C[id])
            c += mean(C[id])

            C[id] = resize!(C[id], length(m1), length(m2))
            distmat_euclidean!(C[id], m1, m2, 2, od.tp, od.c)
            @avx @. C[id] = exp(-λ * C[id])
            c -= 2 * mean(C[id])

            D[i, j] = c
            (D[j, i] = D[i, j])
        end
    end
    normalize && symmetrize!(D)
    Symmetric(D)

end


text:
module ICATT

include("elements.jl")
include("kepler.jl")
include("lambert.jl")
include("dopri.jl")

end

ext: jl

max_issues_repo_name: OpenAstrodynamics/benchmarks

max_forks_repo_name: OpenAstrodynamics/benchmarks


text:
using Plots
using ArgParse
using DrWatson
@quickactivate
using CSV
using BSON
using Random
using FileIO
using DataFrames

# s = datadir("evaluation/zkouska_graf_rel_params.bson")

plotly()
function plot_npar_rank(s)
    f = load(s)
    df = f[:df]
    pl = plot()
    for n in unique(df[:modelname])
        n_df = df[df[:modelname] .== n, :]
        plot!(n_df[:relnpars], n_df[:rank], label = n)
    end
    display(pl)
    return pl
end


text:
function igyp03()
#
#  generalization of Igarash and Ypma
#
    p = poly([10.0*(1+im)*[1;1;1];1;im;2;2im;3;4im;5])
    z = [10*(1+im) 3; 1 1; im 1; 2 1; 2im 1; 3 1; 4im 1; 5 1];
    p, PolyZeros(z)
end


text:
using StableRNGs

include("mlj_functions.jl")

#%%

filename = "./data/df.data"
# N_rows = 1_000_000
N_rows = 10_000

X, y = get_Xy(filename, N_rows);
train, test = partition(eachindex(y), 0.75; shuffle = true, rng = StableRNG(123));

X_train = selectrows(X, train);
y_train = selectrows(y, train);
X_test = selectrows(X, test);
y_test = selectrows(y, test);

hot = OneHotEncoder();
mach = fit!(machine(hot, X))
Xt = MLJ.transform(mach, X)
Xt |> schema

using MLJ
using Flux
using MLDatasets

# helper function
function flatten(x::AbstractArray)
    return reshape(x, :, size(x)[end])
end

import MLJFlux
mutable struct MyConvBuilder
    filter_size::Int
    channels1::Int
    channels2::Int
    channels3::Int
end

function MLJFlux.build(b::MyConvBuilder, rng, n_in, n_out, n_channels)

    k, c1, c2, c3 = b.filter_size, b.channels1, b.channels2, b.channels3

    mod(k, 2) == 1 || error("`filter_size` must be odd. ")

    # padding to preserve image size on convolution:
    p = div(k - 1, 2)

    front = Chain(
        Conv((k, k), n_channels => c1, pad = (p, p), relu),
        MaxPool((2, 2)),
        Conv((k, k), c1 => c2, pad = (p, p), relu),
        MaxPool((2, 2)),
        Conv((k, k), c2 => c3, pad = (p, p), relu),
        MaxPool((2, 2)),
        flatten,
    )
    d = Flux.outputsize(front, (n_in..., n_channels, 1)) |> first
    return Chain(front, Dense(d, n_out))
end

N = 500
Xraw, yraw = MNIST.traindata();
Xraw = Xraw[:, :, 1:N];
yraw = yraw[1:N];

scitype(Xraw)
scitype(yraw)

X = coerce(Xraw, GrayImage);
y = coerce(yraw, Multiclass);

ImageClassifier = @load ImageClassifier
clf = ImageClassifier(
    builder = MyConvBuilder(3, 16, 32, 32),
    epochs = 10,
    loss = Flux.crossentropy,
)

machine(clf, X, y)

evaluate!(
    mach,
    resampling = Holdout(rng = 123, fraction_train = 0.7),
    operation = predict_mode,
    measure = misclassification_rate,
)

max_stars_repo_name: ChristianMichelsen/reference-free

max_issues_repo_name: ChristianMichelsen/reference-free

max_forks_repo_name: ChristianMichelsen/reference-free


text:
# Pendulum Swing-Up Direct Transcription
using LinearAlgebra
const MOI = MathOptInterface

# pendulum dynamics
function pendulum_dynamics(x,u)
    ẋ = zero(x)
    m = 1.
    l = 0.5
    b = 0.1
    lc = 0.5
    I = 0.25
    g = 9.81
    ẋ[1] = x[2]
    ẋ[2] = (u[1] - m*g*lc*sin(x[1]) - b*x[2])/I

    return ẋ
end

n = 2; m = 1 # states, controls

# Problem
x0 = [0.;0.]
xf = [pi;0.]

N = 51 # knot points
NN = (n+m+1)*(N-1) + n # number of decision variables

p = n*N + n + (N-2)# number of equality constraints: dynamics, xf, h_k = h_{k+1}

Q = Diagonal(zeros(n))
R = Diagonal(zeros(m))
Qf = Diagonal(zeros(n))

tf0 = 3.
dt = tf0/(N-1)

u_max = 3.
u_min = -3.

h_max = dt
h_min = dt

function line_trajectory(x0::Vector,xf::Vector,N::Int)
    t = range(0,stop=N,length=N)
    slope = (xf .- x0)./N
    x_traj = [slope*t[k] for k = 1:N]
    x_traj[1] = x0
    x_traj[end] = xf
    x_traj
end

X = line_trajectory(x0,xf,N)
U = [0.01*rand(m) for k = 1:N-1]
H = [dt*ones(1) for k = 1:N-1]

# dynamics constraints
fc(z) = pendulum_dynamics(z[1:n],z[n .+ (1:m)])
fc(x,u) = pendulum_dynamics(x,u)
∇fc(z) = ForwardDiff.jacobian(fc,z)
∇fc(x,u) = ∇fc([x;u])
dfcdx(x,u) = ∇fc(x,u)[:,1:n]
dfcdu(x,u) = ∇fc(x,u)[:,n .+ (1:m)]

F(y,x,u,h) = y - x - h*fc(0.5*(y+x),u)
dFdy(y,x,u,h) = I - 0.5*h*dfcdx(0.5*(y+x),u)
dFdx(y,x,u,h) = -I - 0.5*h*dfcdx(0.5*(y+x),u)
dFdu(y,x,u,h) = -h*dfcdu(0.5*(y+x),u)
dFdh(y,x,u,h) = -fc(0.5*(y+x),u)

function packZ(X,U,H)
    Z = [X[1];U[1];H[1]]

    for k = 2:N-1
        append!(Z,[X[k];U[k];H[k]])
    end
    append!(Z,X[N])

    return Z
end

unpackZ(Z) = let n=n, m=m, N=N, NN=NN
    X = [k != N ? Z[((k-1)*(n+m+1) + 1):k*(n+m+1)][1:n] : Z[((k-1)*(n+m+1) + 1):NN] for k = 1:N]
    U = [Z[((k-1)*(n+m+1) + 1):k*(n+m+1)][n .+ (1:m)] for k = 1:N-1]
    H = [Z[((k-1)*(n+m+1) + 1):k*(n+m+1)][(n+m+1)] for k = 1:N-1]
    return X,U,H
end

unpackZ_timestep(Z) = let n=n, m=m, N=N, NN=NN
    Z_traj = [k != N ? Z[((k-1)*(n+m+1) + 1):k*(n+m+1)] : Z[((k-1)*(n+m+1) + 1):NN] for k = 1:N]
end

g(x,u,h) = let Q=Q, R=R
    (x'*Q*x + u'*R*u) + h
end

dgdx(x,u,h) = let Q=Q
    2*Q*x
end

dgdu(x,u,h) = let R=R
    2*R*u
end

dgdh(x,u,h) = 1.

gf(x) = let Qf=Qf
    x'*Qf*x
end

dgfdx(x) = let Qf=Qf
    2*Qf*x
end

cost(Z) = let n=n, m=m, N=N, NN=NN, Q=Q, R=R, Qf=Qf
    J = 0.
    Z_traj = unpackZ_timestep(Z)

    for k = 1:N-1
        y = Z_traj[k+1][1:n]
        x = Z_traj[k][1:n]
        u = Z_traj[k][n .+ (1:m)]
        h = Z_traj[k][n+m+1]
        J += g(0.5*(y+x),u,h)
    end
    xN = Z_traj[N][1:n]
    J += gf(xN)

    return J
end

∇cost!(∇J,Z) = let n=n, m=m, N=N, NN=NN, Q=Q, R=R, Qf=Qf
    Z_traj = unpackZ_timestep(Z)
    k = 1
    y = Z_traj[k+1][1:n]
    x = Z_traj[k][1:n]
    u = Z_traj[k][n .+ (1:m)]
    h = Z_traj[k][n+m+1]
    xm = 0.5*(y+x)
    ∇J[1:(n+m+1+n)] = [0.5*dgdx(xm,u,h);dgdu(xm,u,h);dgdh(xm,u,h);0.5*dgdx(xm,u,h)]

    for k = 2:N-1
        y = Z_traj[k+1][1:n]
        x = Z_traj[k][1:n]
        u = Z_traj[k][n .+ (1:m)]
        h = Z_traj[k][n+m+1]
        xm = 0.5*(y+x)
        ∇J[((k-1)*(n+m+1) + 1):(k*(n+m+1)+n)] = [0.5*dgdx(xm,u,h);dgdu(xm,u,h);dgdh(xm,u,h);0.5*dgdx(xm,u,h)]
    end
    xN = Z_traj[N][1:n]
    ∇J[(NN - (n-1)):NN] = dgfdx(xN)

    return nothing
end

constraints!(g,Z) = let n=n, m=m, N=N, NN=NN, x0=x0, xf=xf,p=p
    Z_traj = unpackZ_timestep(Z)
    k = 1
    y = Z_traj[k+1][1:n]
    x = Z_traj[k][1:n]
    u = Z_traj[k][n .+ (1:m)]
    h = Z_traj[k][n+m+1]
    con = x-x0

    for k = 1:N-1
        y = Z_traj[k+1][1:n]
        x = Z_traj[k][1:n]
        u = Z_traj[k][n .+ (1:m)]
        h = Z_traj[k][n+m+1]
        append!(con,F(y,x,u,h))

        if k != N-1
            h₊ = Z_traj[k+1][n+m+1]
            append!(con,h - h₊)
        end
    end
    xN = Z_traj[N][1:n]
    append!(con,xN-xf)

    # g = con
    copyto!(view(g,1:p),con)
    return nothing
end

∇constraints_vec!(∇con,Z) = let n=n, m=m, N=N, NN=NN, x0=x0, xf=xf, p=p
    Z_traj = unpackZ_timestep(Z)

    shift = 0
    # con = x-x0
    copyto!(view(∇con,1:n^2),vec(Diagonal(ones(n))))
    shift += n^2

    for k = 1:N-1
        y = Z_traj[k+1][1:n]
        x = Z_traj[k][1:n]
        u = Z_traj[k][n .+ (1:m)]
        h = Z_traj[k][n+m+1]

        # dynamics constraint jacobians
        copyto!(view(∇con,shift .+ (1:(n+m+1+n)*(n))),vec([dFdx(y,x,u,h) dFdu(y,x,u,h) dFdh(y,x,u,h) dFdy(y,x,u,h)]))
        shift += (n+m+1+n)*n

        if k != N-1
            copyto!(view(∇con,shift .+ (1:2)),[1.;-1.])
            shift += 2
        end
    end
    xN = Z_traj[N][1:n]
    copyto!(view(∇con,shift .+ (1:n^2)),vec(Diagonal(ones(n))))

    return nothing
end

constraint_sparsity() = let n=n, m=m, N=N, NN=NN, x0=x0, xf=xf, p=p

    row = []
    col = []

    r = 1:n
    c = 1:n

    for cc in c
        for rr in r
            push!(row,convert(Int,rr))
            push!(col,convert(Int,cc))
        end
    end

    for k = 1:N-1
        # dynamics constraint jacobians
        c_idx = ((k-1)*(n+m+1)) .+ (1:(n+m+1+n))
        r_idx = (n + (k-1)*(n + 1)) .+ (1:n)

        for cc in c_idx
            for rr in r_idx
                push!(row,convert(Int,rr))
                push!(col,convert(Int,cc))
            end
        end

        if k != N-1
            c_idx = ((k-1)*(n+m+1) + n + m) .+ 1
            c2_idx = ((k)*(n+m+1) + n + m) .+ 1
            r_idx = (n + n + (k-1)*(n + 1)) .+ 1

            for cc in c_idx
                for rr in r_idx
                    push!(row,convert(Int,rr))
                    push!(col,convert(Int,cc))
                end
            end

            for cc in c2_idx
                for rr in r_idx
                    push!(row,convert(Int,rr))
                    push!(col,convert(Int,cc))
                end
            end
        end

    end
    k = N
    r = (p - (n-1)):p
    # println(r)
    r = ((n + (k-1)*(n+1)-1) .+ (1:n))
    # println(r)
    c = (NN - (n-1)):NN

    for cc in c
        for rr in r
            push!(row,convert(Int,rr))
            push!(col,convert(Int,cc))
        end
    end

    return collect(zip(row,col))
end

primal_bounds() = let n=n, m=m, N=N, NN=NN, u_max=u_max, u_min=u_min, h_max=h_max, h_min=h_min
    Z_low = [-Inf*ones(n);u_min*ones(m);h_min]
    Z_up = [Inf*ones(n);u_max*ones(m);h_max]

    for k = 2:N-1
        append!(Z_low,[-Inf*ones(n);u_min*ones(m);h_min])
        append!(Z_up,[Inf*ones(n);u_max*ones(m);h_max])
    end
    append!(Z_low,-Inf*ones(n))
    append!(Z_up,Inf*ones(n))

    return Z_low, Z_up
end

constraint_bounds() = let p=p
    c_low = zeros(p); c_up = zeros(p)

    return c_low, c_up
end

## set up optimization (MathOptInterface)
ZZ = packZ(X,U,H)

struct Problem <: MOI.AbstractNLPEvaluator
    enable_hessian::Bool
end

MOI.features_available(prob::Problem) = [:Grad, :Jac]
MOI.initialize(prob::Problem, features) = nothing
MOI.jacobian_structure(prob::Problem) = constraint_sparsity()
MOI.hessian_lagrangian_structure(prob::Problem) = []

function MOI.eval_objective(prob::Problem, Z)
    return cost(Z)
end

function MOI.eval_objective_gradient(prob::Problem, grad_f, Z)
    ∇cost!(grad_f, Z)
end

function MOI.eval_constraint(prob::Problem, g, Z)
    constraints!(g,Z)
end

function MOI.eval_constraint_jacobian(prob::Problem, jac, Z)
    ∇constraints_vec!(jac,Z)
end

MOI.eval_hessian_lagrangian(prob::Problem, H, x, σ, μ) = nothing

prob = Problem(false)

Z_low, Z_up = primal_bounds()
c_low, c_up = constraint_bounds()

nlp_bounds = MOI.NLPBoundsPair.(c_low,c_up)
block_data = MOI.NLPBlockData(nlp_bounds,prob,true)

solver = SNOPT7.Optimizer()

Z = MOI.add_variables(solver,NN)

for i = 1:NN
    zi = MOI.SingleVariable(Z[i])
    MOI.add_constraint(solver, zi, MOI.LessThan(Z_up[i]))
    MOI.add_constraint(solver, zi, MOI.GreaterThan(Z_low[i]))
    MOI.set(solver, MOI.VariablePrimalStart(), Z[i], ZZ[i])
end

# Solve the problem
MOI.set(solver, MOI.NLPBlock(), block_data)
MOI.set(solver, MOI.ObjectiveSense(), MOI.MIN_SENSE)
MOI.optimize!(solver)

# Get the solution
res = MOI.get(solver, MOI.VariablePrimal(), Z)

X,U,H = unpackZ(res)
Xblk = zeros(n,N)
for k = 1:N
    Xblk[:,k] = X[k]
end

Ublk = zeros(m,N-1)
for k = 1:N-1
    Ublk[:,k] = U[k]
end

plot(Xblk')
plot(Ublk')


text:
#####################
function read_results(
    results_path::String;
    control_vars::Array = [],
    control_val_equal::Array = [],
    vars_to_check::Array = [])

    """
	Description:
		read results from path & save to dict. returns dict and filtered list of parametrizations.
	"""

# reading files
    params_sets = CSV.read("$results_path/sim_params_comb_df.csv", DataFrame)

    if length(control_vars) > 0
        # filtering params:
        for (cond, val) in zip(control_vars, control_val_equal)
            params_sets = params_sets[params_sets[:,cond] .== val, :]
        end
        if length(vars_to_check) > 0
                @assert size(unique(params_sets[!, Symbol.(vars_to_check)]))[1] == size(params_sets)[1] "Incorrect filtering conditions specified. Multiple or none scenarios found."
        end
    end

    full_output = Dict()
    file_list = glob("simul_res_parset_id_*.json", results_path)
    for id in params_sets[:,:parset_id]
        if string(results_path,"/simul_res_parset_id_",id,".json") ∈ file_list
            full_output[id] = JSON.parsefile(string(results_path,"/simul_res_parset_id_",id,".json"))
            params = filter(r -> in(r.parset_id, [id]), params_sets)
            for col in names(params)
                full_output[id][col] = params[1,col]
            end
        else
            filter!(r -> !in(r.parset_id, [id]), params_sets)
        end
    end
    return full_output, params_sets
end

function develplot_by_var(response_var::String, variable::String, results::Dict;
    total_indic::Union{String, Nothing} = nothing,
    div::Union{Float64, Int} = 1,
    var_unit::String = "",
    max_it::Union{Int, Nothing} = nothing,
    min_it::Int = 1,
    x_lab::String = "iteration",
    y_lab::Union{String, Nothing} = nothing,
    title::Union{String, Nothing} = nothing,
    with_std::Bool = false,
    out_name::String="develplot",
    with_export::Bool = true)

    """
	Description:
		Plots line series with or without std
	"""
    full_output = deepcopy(results)

    #gathering a df
    df = DataFrame()
    stds = DataFrame()
    var_vals = []
    for (key, value) in full_output
        symb = value[variable]/div
        if !isnothing(total_indic) && length(value[total_indic])>1
            value[total_indic] = value[total_indic][1]
            if total_indic == "prct_of_agents_used_TTC"
                value[total_indic] = value[total_indic]/value["N_agents"]
            end
        end
        push!(var_vals, symb)
        if response_var in ["total_infected", "TTC_infected", "street_infected"]
            df[:,Symbol(symb)] = value[response_var][1]./(isnothing(total_indic) ? 1 : value[total_indic])
            stds[:,Symbol(symb)] = value[response_var][2]./(isnothing(total_indic) ? 1 : value[total_indic])
        else
            df[:,Symbol(symb)] = value[response_var]./(isnothing(total_indic) ? 1 : value[total_indic])
        end
    end

    var_vals = sort(collect(var_vals))

    pushfirst!(var_vals, rand())
    cs1 = ColorScheme(range(colorant"orange", colorant"purple", length = length(var_vals)))
    cs1 = Dict(zip(var_vals, cs1))
    var_vals = var_vals[2:end]
    max_it =  isnothing(max_it) ? size(df)[1] : max_it

    y_lab =  isnothing(y_lab) ? response_var : y_lab
    title =  isnothing(title) ? "$response_var with respect to $variable" : title

    gr(size = (800, 550))
    my_plot = Plots.plot()
    for i in var_vals
        rounding = round(i, digits = 1)
        plot!(df[min_it:max_it, Symbol(i)], linewidth = 2.5,
        legend = :topleft, lab = "$rounding $var_unit",
            linecolor = cs1[i], annotations = ([max_it + 240], [df[max_it, Symbol(i)]],
                    Plots.text("$rounding", "Courier", :purple, 8)))
        if with_std == true
            plot!((df[min_it:max_it, Symbol(i)] + stds[min_it:max_it, Symbol(i)]), linewidth = 2.5,
            legend = false,  linestyle = :dot, linecolor = cs1[i])

            plot!((df[min_it:max_it, Symbol(i)] - stds[min_it:max_it, Symbol(i)]), linewidth = 2.5,
            legend = false,  linestyle = :dot, linecolor = cs1[i])
        end
    end
    xlabel!(x_lab)
    ylabel!(y_lab)
    title!(title)
    if with_export == true
        #n = now()
        Plots.savefig(my_plot,"./viz/$out_name.pdf")
    end
    #my_plot
end

#####################
function point_in_time_by_var(response_vars::Array, variable::String, results::Dict;
    points_in_time::Array{} = [],
    total_indic::Union{String, Nothing} = nothing,
    div::Union{Float64, Int} = 1,
    time_div::Union{Float64, Int} = 1,
    time_unit::String = "",
    x_lab::Union{String, Nothing} = nothing,
    y_lab::Union{String, Nothing} = nothing,
    title::Union{String, Nothing} = nothing,
    out_name::String = "point_in_time",
    with_export::Bool = true)

    """
	Description:
		Line plots for a numerical parameter shot at a given point(or points) of time
	"""
    @assert min(length(response_vars),length(points_in_time)) ==  1 "Only one array longer than 1 allowed"

    pyplot()

    full_output = deepcopy(results)

    #gathering a df
    df = Dict()
    for (key, value) in full_output
        symb = value[variable]/div
        if !isnothing(total_indic) && length(value[total_indic])>1
            value[total_indic] = value[total_indic][1]
            if total_indic == "prct_of_agents_used_TTC"
                value[total_indic] = value[total_indic]*value["N_agents"]/100
            end
        end
        for response_var in response_vars
            if response_var ∉ keys(df)
                df[response_var] = DataFrame()
            end
            if response_var in ["total_infected", "TTC_infected", "street_infected"]
                df[response_var][:,Symbol(symb)] = 100value[response_var][1]./(isnothing(total_indic) ? 1 : value[total_indic])
            else
                df[response_var][:,Symbol(symb)] = value[response_var]./(isnothing(total_indic) ? 1 : value[total_indic])
            end
        end
    end

    points_in_time = length(points_in_time) > 0 ? points_in_time : [size(df[response_vars][1])[1]]
    points_in_time = sort(points_in_time)
    points_in_time = trunc.(Int, points_in_time)

    results = Dict()
    for response_var in response_vars
        results[response_var] = df[response_var][points_in_time, :]
        results[response_var] = DataFrame([[names(results[response_var])]; collect.(eachrow(results[response_var]))], [Symbol(variable); Symbol.("ITER_", points_in_time)])
        results[response_var][!,variable] = parse.(Float64, results[response_var][!,variable])
        sort!(results[response_var], Symbol(variable))
    end

    first_var = response_vars[1]
    x_lab =  isnothing(x_lab) ? variable : x_lab
    y_lab =  isnothing(y_lab) ? "infected" : y_lab
    title =  isnothing(title) ? "$first_var with respect to $variable" : title

    if length(points_in_time) >= length(response_vars)
        pushfirst!(points_in_time, -1)
        cs1 = ColorScheme(range(colorant"orange", colorant"purple", length = length(points_in_time)))
        cs1 = Dict(zip(points_in_time, cs1))
        points_in_time = points_in_time[2:end]
    else
        pushfirst!(response_vars, "aaa")
        cs1 = ColorScheme(range(colorant"orange", colorant"purple", length = length(response_vars)))
        cs1 = Dict(zip(response_vars, cs1))
        response_vars = response_vars[2:end]
    end

    #gr(size = (800, 600))
    my_plot = Plots.plot()

    for var in response_vars
        for i in points_in_time
            key_name = length(points_in_time) >= length(response_vars) ? i : var
            rounding = trunc(Int, i/time_div)
            lab = length(points_in_time) >= length(response_vars) ? "$rounding $time_unit" : var
            plot!(results[var][:, Symbol(variable)], results[var][!, Symbol("ITER_", i)], linewidth = 2.5,
            #legend = :topleft, lab = lab,
            legend = false,
            grid=false,
            linecolor = cs1[key_name])
        end
    end

    xlabel!(x_lab)
    ylabel!(y_lab)
    title!(title)
    if with_export == true
        #n = now()
        Plots.savefig(my_plot,"./viz/$out_name.eps")
    end
    #my_plot
end

function make_heatmap(response_var::String, variable1::String,
    variable2::String, results::Dict, params_sets::DataFrame;
    total_indic::Union{String, Nothing} = nothing,
    div1::Union{Float64, Int} = 1,
    div2::Union{Float64, Int} = 1,
    out_name::String="heatmap",
    with_export::Bool = true)

    """
	Description:
		Plots a heatmap
	"""
    full_output = deepcopy(results)

    #gathering a matrix
    var1 = []
    var2 = []
    resp_var = Array{Union{Float64, Int, Missing}}(missing,
    length(unique(params_sets[!, Symbol(variable1)])), length(unique(params_sets[!, Symbol(variable2)])))

    for (key, value) in full_output
        if (length(findall(x -> x == value[variable1]/div1, var1 )) == 0)
            push!(var1, value[variable1]/div1)
        end
        if (length(findall(x -> x == value[variable2]/div2, var2)) == 0)
            push!(var2, value[variable2]/div2)
        end
    end
    sort!(var1)
    sort!(var2)

    for (key, value) in full_output

        row_pos = findall(x -> x == value[variable1]/div1, var1)[1]
        col_pos = findall(x -> x == value[variable2]/div2, var2)[1]

        val = value[response_var][1]/(isnothing(total_indic) ? 1 : value[total_indic])
        resp_var[row_pos, col_pos] = val
    end

    #backend(:plotly)
    plotly_js()
    my_plot = Plots.heatmap(var2,var1,resp_var)
    xlabel!(variable2)
    ylabel!(variable1)
    title!(response_var)

    if with_export == true
        Plots.savefig(my_plot,"./viz/$out_name.eps")
    end
end

function plot_XY(variable1::String, variable2::String, results::Dict;
    point_in_time::Union{Float64, Int, Nothing} = nothing,
    total_indic::Union{String, Nothing} = nothing,
    x_lab::Union{String, Nothing} = nothing,
    y_lab::Union{String, Nothing} = nothing,
    title::Union{String, Nothing} = nothing,
    out_name::String = "XY",
    with_export::Bool = true)

    full_output = deepcopy(results)

    pyplot()

    #gathering a df
    df = DataFrame([Float64, Float64], Symbol.([variable1,variable2]))
    for (key, value) in full_output
        val = deepcopy(value)
        for var in [variable1, variable2]
            if length(val[var]) > 1
                val[var] = val[var][1]
            end
            if var in ["total_infected", "TTC_infected", "street_infected"]
                val[var] = isnothing(point_in_time) ? val[var][end] : val[var][Int(point_in_time)]
                val[var] = val[var]/(isnothing(total_indic) ? 1 : value[total_indic])
            end
            if var == "prct_of_agents_used_TTC"
                val[var] = val[var]/100
            end
        end
        push!(df, Dict(variable1 => val[variable1], variable2 => val[variable2]))
    end
    sort!(df, variable1)

    #gr(size = (800, 600))
    x_lab =  isnothing(x_lab) ? variable1 : x_lab
    y_lab =  isnothing(y_lab) ? variable2 : y_lab
    title =  isnothing(title) ? "$variable2 with respect to $variable2" : title

    my_plot = Plots.plot(df[!, variable1], df[!, variable2], linewidth = 2.5, legend = false,
    linecolor="purple", grid=false)
    xlabel!(x_lab)
    ylabel!(y_lab)
    title!(title)

    if with_export == true
        Plots.savefig(my_plot,"./viz/$out_name.eps")
    end
    #my_plot
end

function plot_XYY(variable1::String, variable2::String, x_var::String,  results::Dict;
    point_in_time::Union{Float64, Int, Nothing} = nothing,
    divx::Union{Float64, Int} = 1,
    total_indic::Union{String, Nothing} = nothing,
    x_lab::Union{String, Nothing} = nothing,
    y1_lab::Union{String, Nothing} = nothing,
    y2_lab::Union{String, Nothing} = nothing,
    title::Union{String, Nothing} = nothing,
    out_name::String = "XYY",
    with_export::Bool = true)

    pyplot()

    full_output = deepcopy(results)
    #gathering a df
    df = DataFrame([Float64, Float64, Float64], Symbol.([variable1, variable2, x_var]))
    for (key, value) in full_output
        val = deepcopy(value)
        val[x_var] = val[x_var]/divx
        for var in [variable1, variable2]
            if length(val[var]) > 1
                val[var] = val[var][1]
            end
            if var in ["total_infected", "TTC_infected", "street_infected"]
                val[var] = isnothing(point_in_time) ? val[var][end] : val[var][Int(point_in_time)]
                val[var] = val[var]/(isnothing(total_indic) ? 1 : value[total_indic])
            end
            if var == "prct_of_agents_used_TTC"
                val[var] = val[var]/100
            end
        end
        push!(df, Dict(variable1 => val[variable1], variable2 => val[variable2], x_var => val[x_var]))
    end
    sort!(df, x_var)

    #gr(size = (800, 600))
    x_lab =  isnothing(x_lab) ? x_var : x_lab
    y1_lab =  isnothing(y_lab) ? variable1 : y1_lab
    y2_lab =  isnothing(y_lab) ? variable2 : y2_lab
    title =  isnothing(title) ? "$variable1 & $variable2 with respect to $x_var" : title

    my_plot = Plots.plot(df[!, x_var], df[!, variable1], linewidth = 2.5, linecolor="purple",
    label=y1_lab,legend=:topleft, xlabel=x_lab)
    plot!(twinx(), df[!, x_var], df[!, variable2], xticks=:none, linecolor="black", grid=false, linewidth = 2.5, label=y2_lab)

    title!(title)

    if with_export == true
        Plots.savefig(my_plot,"./viz/$out_name.eps")
    end
    #my_plot
end

max_stars_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_issues_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_forks_repo_name: NykPol/EpidemicInUrbanNetworkToronto


text:
function choleskydec(g::AbstractGraph; kwargs...)
    matrix, permutation, nb_added_edge = choleskyfact(g)
    data = Dict("matrix" => matrix,
                "permutation" => permutation,
                "nb_added_edge" => nb_added_edge)
    return SimpleGraph(matrix), data
end

function sparsity_pattern(g::AbstractGraph)
    sparsity_pattern = adjacency_matrix(g)
    for i in 1:nv(g)
        sparsity_pattern[i, i] = length(neighbors(g, i)) + 1
    end
    return sparsity_pattern
end

function choleskyfact(g::AbstractGraph; kwargs...)
    A = sparsity_pattern(g)
    nb_edges_A = (nnz(A) - size(A, 1)) / 2
    # NOTE: AMD ordering by default
    F = cholesky(A; kwargs...) #NOTE: order = F.p
    # computing L + LT
    L = sparse(F.L)
    nb_edges_L = nnz(L) - size(A,1)
    nb_added_edge = nb_edges_L - nb_edges_A
    SP = L + L'
    #inverting permutation to get chordal extension of sparsity_pattern
    H = SP[invperm(F.p), invperm(F.p)]
    return H, F.p, nb_added_edge
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
function get_gm(obj::FieldDSSTNewtonianAttractionContext)
    return jcall(obj, "getGM", RealFieldElement, ())
end


text:
module RocketSomeOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: some()" begin

    println("Testing: operator some()")

    run_proxyshowcheck("Some", some())

    run_testset([
        (
            source = from(1:42) |> max() |> some(),
            values = @ts([ 42, c ]),
            source_type = Int,
        ),
        (
            source = from([ 2, 3, 4, nothing ]) |> some(),
            values = @ts([ 2, 3, 4, c ]),
            source_type = Int,
        ),
        (
            source = completed() |> skip_next(),
            values = @ts(c)
        ),
        (
            source = completed(Int) |> max() |> some(),
            values = @ts(c),
            source_type = Int
        ),
        (
            source = faulted(1) |> skip_next(),
            values = @ts(e(1))
        ),
        (
            source = never() |> skip_next(),
            values = @ts()
        )
    ])

    @testset begin
        @test_throws Exception of(1) |> some()
    end

end

end


text:
abstract Container <: Compose3DNode

type Context <: Container
	box::BoundingBox #Parent box.
	geometry_children::List{Geometry}
	material_children::List{Material}
	container_children::List{Container}
	light_children::List{Light}
	camera::Union{Camera, Void}
end

Context(x0::Length,y0::Length,z0::Length,width::Length,height::Length,depth::Length) =
    Context(
			BoundingBox(x0,y0,z0,width,height,depth),
			ListNull{Geometry}(),
			ListNull{Material}(),
			ListNull{Container}(),
			ListNull{Light}(),
			nothing
		)

Context(ctx::Context) = Context(
		ctx.box,
		ctx.geometry_children,
		ctx.material_children,
		ctx.container_children,
		ctx.light_children,
		ctx.camera
)

function copy(ctx::Context)
    return Context(ctx)
end

#Compositions.

function compose!(a::Context, b::Container)
    a.container_children = cons(b, a.container_children)
    return a
end

function compose!(a::Context, b::Geometry)
    a.geometry_children = cons(b, a.geometry_children)
    return a
end

function compose!(a::Context, b::Material)
    a.material_children = cons(b, a.material_children)
    return a
end

function compose!(a::Context, b::Light)
    a.light_children = cons(b, a.light_children)
    return a
end

function compose!(a::Context, b::Camera)
    a.camera = b
    return a
end

function compose(a::Context, b::Compose3DNode)
    a = copy(a)
    compose!(a, b)
    return a
end

# higher-order compositions
function compose!(a::Context, b, c, ds...)
    return compose!(compose!(a, b), c, ds...)
end

function compose!(a::Context, bs::AbstractArray)
    compose!(a, compose!(bs...))
end

function compose!(a::Context, bs::Tuple)
    compose!(a, compose!(bs...))
end

function compose!(a::Context)
    return a
end

function compose(a::Context, b, c, ds...)
    return compose(compose(a, b), c, ds...)
end

function compose(a::Context, bs::AbstractArray)
    compose(a, compose(bs...))
end

function compose(a::Context, bs::Tuple)
    compose(a, compose(bs...))
end

function compose(a::Context)
    return a
end

function draw(backend::Backend, root_container::Context)
    drawpart(backend, root_container, root_box(backend))
    finish(backend)
end

# Draw without finishing the backend
#
# Drawing is basically a depth-first traversal of the tree, pushing and popping
# properties, expanding context promises, etc. as needed.
#
function drawpart(backend::Backend, container::Container,
                  parent_box::Absolute3DBox)

    # used to collect property children
    properties = Array(Material, 0)
    ctx = container
    box = resolve(parent_box, ctx.box)

    child = ctx.material_children

    while !isa(child, ListNull)
        push!(properties, resolve(parent_box, child.head))
        child = child.tail
    end

    if !isempty(properties)
        push_property_frame(backend, properties)
    end

    child = ctx.geometry_children
    while !isa(child, ListNull)
        draw(backend, box, child.head)
        child = child.tail
    end

		child = ctx.light_children
    while !isa(child, ListNull)
        draw(backend, box, child.head)
        child = child.tail
    end

		if ctx.camera!=nothing
				draw(backend, box, ctx.camera)
		end

    child = ctx.container_children
    while !isa(child, ListNull)
        drawpart(backend, child.head, box)
        child = child.tail
    end

    if !isempty(properties)
        pop_property_frame(backend)
    end
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

let
    posterior = BAT.example_posterior()

    dummy_samples = bat_sample(BAT.bat_determ_rng(), posterior.prior, IIDSampling(nsamples=10)).result
    precompile(bat_eff_sample_size, map(typeof, (dummy_samples,)))
    precompile(show, map(typeof, (dummy_samples,)))

    precompile(SampledDensity, map(typeof, (posterior, dummy_samples)))

    for mcalg in (MetropolisHastings(), HamiltonianMC())
        precompile(bat_sample, map(typeof, (posterior, MCMCSampling(mcalg = mcalg))))
    end
end


text:
abstract Individual

type Individual{F}
  fitness::F
  stage::
end

type SimpleIndividual <: Individual
    fitness::F
    genome::
end

type Deme
  stages::Vector{Any}
end

# population:
# stages["genome"] = BitVector[]


text:
using GeneralizedDynamicsFromData
using OrderedCollections
using JLD2
using FileIO

experiment_name = "genetic_toggle_statistics_MSE_new"
repetitions = 100

# weight_decay = 1e-4
η = 1e-2
# η_decay_rate = 0.5
# η_decay_step = 100
# η_limit = 1e-4
loss_weights = [1/2, 1/2]

# construct_optimiser() = Flux.Optimiser(WeightDecay(weight_decay),
#                                        ExpDecay(η,η_decay_rate,η_decay_step,η_limit),
#                                        ADAM())
construct_optimiser() = ADAM(η)
construct_loss(θ, y, predict) = polar_loss(θ, y, loss_weights, predict)

net_config = OrderedDict([
    :inputs => 2,
    :outputs => 1,
    :neurons => 16,
    :layers => 1,
    :non_lin => tanh,
    :initialization => Flux.glorot_normal
])

problem_one_stable = Dict([:equation => genetic_toggle_switch,
                           :parameters => Float64[1.5, 1.5, 2.0, 2.0],
                           :u0 => Float64[2.0, 4.0],
                           :tspan => (0.0f0, 15.0f0),
                           :ts => 0.1,
                           :solver => Tsit5,
                           :optimizer => construct_optimiser,
                           :max_iter => 3000,
                           :loss => construct_loss]
                      )

for noise in [1e-2, 5e-3, 1e-3, 1e-4]
  summary, callbacks = repeat_experiment(problem_one_stable,
                                        net_config,
                                        repetitions;
                                        ε = noise,
                                        progress=false)
  filename_cb = joinpath("./data/", experiment_name*"_1s_"*string(noise)*"_all.h5")
  filename_sum = joinpath("./data/", experiment_name*"_1s_"*string(noise)*".jld2")
  callbacks_to_hdf5(callbacks, filename_cb)
  save(filename_sum, summary)
end

problem_two_stable = Dict([:equation => genetic_toggle_switch,
                           :parameters => Float64[3.5, 3.5, 2.0, 2.0],
                           :u0 => Float64[2.0, 4.0],
                           :tspan => (0.0f0, 15.0f0),
                           :ts => 0.1,
                           :solver => Tsit5,
                           :optimizer => construct_optimiser,
                           :max_iter => 3000,
                           :loss => construct_loss]
                         )

for noise in [1e-2, 5e-3, 1e-3, 1e-4]
  summary, callbacks = repeat_experiment(problem_two_stable,
                                        net_config,
                                        repetitions;
                                        ε = noise,
                                        progress=false)
  filename_cb = joinpath("./data/", experiment_name*"_2s_"*string(noise)*"_all.h5")
  filename_sum = joinpath("./data/", experiment_name*"_2s_"*string(noise)*".jld2")
  callbacks_to_hdf5(callbacks, filename_cb)
  save(filename_sum, summary)
end

max_stars_repo_name: pnieters/GeneralizedDynamicsFromData

max_issues_repo_name: pnieters/GeneralizedDynamicsFromData

max_forks_repo_name: pnieters/GeneralizedDynamicsFromData


text:
"""
# solvers for multi-robot motion planning

## included solvers
- PRM: probabilistic roadmap
- CBS: PRM + conflict-based search
- PP: PRM + prioritized planning
- RRT: rapidly-exploring random tree
- RRT-connect
- SSSP: simultaneous sampling-and-search planning
"""
module Solvers

using MRMP
import ..MRMP: AbsState, Node, dist, now, elapsed_sec, gen_uniform_sampling, get_mid_status
import Printf: @sprintf, @printf
import Base: @kwdef
import DataStructures: PriorityQueue, enqueue!, dequeue!
import Random: randperm

include("utils.jl")
include("prm.jl")
include("rrt.jl")
include("rrt_connect.jl")
include("pp.jl")
include("cbs.jl")
include("sssp.jl")

using .LibPP, .LibCBS, .LibRRT, .LibRRT_connect, .LibSSSP
export PRM_direct, PRMs, PRMs!, RRT, RRT_connect, PP, CBS, SSSP
end


text:
using Raylib_jll

struct RayColor
    r::Cuchar
    g::Cuchar
    b::Cuchar
    a::Cuchar
end

const RAYWHITE = RayColor(245, 245, 245, 255)
const LIGHTGRAY = RayColor(200, 200, 200, 255)
const DARKGRAY = RayColor(80, 80, 80, 255)
const MAROON = RayColor(190, 33, 55, 255)
const RED = RayColor(230, 41, 55, 255)

struct RayVector3
    x::Cfloat
    y::Cfloat
    z::Cfloat
end

struct RayCamera3D
    position::RayVector3
    target::RayVector3
    up::RayVector3
    fovy::Cfloat
    projection::Cint
end

function main()
    screenWidth = 800
    screenHeight = 450

    @ccall libraylib.InitWindow(screenWidth::Cint, screenHeight::Cint,
                             "raylib [core] example - 3d camera mode"::Cstring)::Cvoid

    camera = RayCamera3D(
        RayVector3(0, 10, 10),
        RayVector3(0, 0, 0),
        RayVector3(0, 1, 0),
        45,
        0,
    )

    cubePosition = RayVector3(0, 0, 0)

    @ccall libraylib.SetTargetFPS(60::Cint)::Cvoid

    while iszero(@ccall libraylib.WindowShouldClose()::Cint)
        @ccall libraylib.BeginDrawing()::Cvoid
        @ccall libraylib.ClearBackground(RAYWHITE::RayColor)::Cvoid

        @ccall libraylib.BeginMode3D(camera::RayCamera3D)::Cvoid
        @ccall libraylib.DrawCube(cubePosition::RayVector3, 2f0::Cfloat, 2f0::Cfloat, 2f0::Cfloat,
                               RED::RayColor)::RayColor

        @ccall libraylib.DrawCubeWires(cubePosition::RayVector3, 2f0::Cfloat, 2f0::Cfloat, 2f0::Cfloat,
                                    MAROON::RayColor)::Cvoid

        @ccall libraylib.DrawGrid(10::Cint, 1f0::Cfloat)::Cvoid

        @ccall libraylib.EndMode3D()::Cvoid
        @ccall libraylib.DrawText("Welcome to the third dimension!"::Cstring,
                               10::Cint, 40::Cint, 20::Cint, DARKGRAY::RayColor)::Cvoid
        @ccall libraylib.EndDrawing()::Cvoid
    end

    @ccall libraylib.CloseWindow()::Cvoid

end


text:
#function spectral_line(x::AbstractArray; center::Real=0.0, width::Real=0.0, depth::Real=1.0, max_dist::Real=5.0 )
#   return spectral_line.(x,center=center,width=width,depth=depth, max_dist=max_dist)
#end

function spectral_line(x::Real; center::Real=0.0, width::Real=0.0, depth::Real=1.0, max_dist::Real=5.0 )
   return (x-center)/width > max_dist ? 1.0 : 1.0-depth*exp(-0.5*((x-center)/width)^2)
end

function spectral_line(x::Real; center::AbstractArray, width::AbstractArray, depth::AbstractArray, max_dist::Real=5.0 )
   @assert length(center) == length(width) == length(depth)
   output = one(typeof(x))
   for i in 1:length(center)
      output *= spectral_line(x,center=center[i],width=width[i],depth=depth[i], max_dist=max_dist)
   end
   return output
end

function make_toy_data(obs_size::Tuple{Int64,Int64},num_obs::Int = 0)
   instrument = :HPF
   data = EchelleObservationSet(obs_size, num_obs, instrument=instrument)
   lambda_min = 4000.0
   lambda_max = 8000.0
   num_lines = 1000
   line_widths = fill(5.0,num_lines)
   line_depths = fill(0.5,num_lines)
   line_centers = lambda_min.+(lambda_max-lambda_min).*rand(num_lines)

   snr = 100
   (num_obs, num_orders) = obs_size
   for i in 1:num_orders
      data.lambda[:,i,:] .= range(lambda_min+(lambda_max-lambda_min)*(i-1)/num_orders,
         stop=lambda_min+(lambda_max-lambda_min)*(i)/num_orders,length=obs_size[1])
   end
   data.flux .= spectral_line.(data.lambda, center=line_centers, width=line_widths, depth=line_depths, max_dist=5.0)
   data.var .= 1.0/snr
   data.flux .*= 1.0 .+ randn(size(data.flux))./snr
   return data
end


text:
module FEASTsolver

using LinearAlgebra, FastGaussQuadrature

export feast, cauchynodes

"""
    feast(A, M0, λmin, λmax; cpts=8, tol=1.0e-12, maxiter=20)
    feast(A, Y, zs, ws, λmin, λmax; tol=1.0e-12, maxiter=20)

Determines all the eigenvalues and vectors of an Hermitian matrix `A` inside the interval `]λmin, λmax]`.

`tol` is the relative error tolerance. The error is estimated using the variation of the trace of the eigenvalues in `]λmin, λmax]`.
`maxiter` is the number of FEAST outer iterations.

`M0` is the size of the subspace search. It must be a bit larger than the number of eigenvalues in the interval.
It is also possible to provide an initial matrix `Y` as the initial search basis. In that case `M0=size(Y, 2)`.

`cpts` is the number of Gauss-Legendre points using in the approximation of the Cauchy contour integral.
The nodes `zs` and weights `ws` for the quadrature rule can also be passed explicitly.
"""
function feast(A, M0::Integer, λmin::Real, λmax::Real; cpts::Integer=8, tol::Real=1.0e-12, maxiter::Integer=20)
    # initial search space of size M0
    V0 = rand(complex(eltype(A)), size(A, 1), M0)
    # contour integral nodes and and weights
    zs, ws = cauchynodes(λmin, λmax, cpts)

    return feast(A, V0, zs, ws, λmin, λmax; tol=tol, maxiter=maxiter)
end

function feast(A, Y, zs, ws, λmin::Real, λmax::Real; tol::Real=1.0e-12, maxiter::Integer=20)

    # dimension of problem and of subspace search
    N, M0 = size(Y)

    # allocate output
    ritzvals = Vector{Complex{eltype(A)}}(undef, M0)
    ritzvecs = Matrix{Complex{eltype(A)}}(undef, M0, M0)

    # allocate computation arrays
    Q = similar(Y)
    q = Vector{eltype(Y)}(undef, N)
    AQ = Matrix{eltype(Y)}(undef, M0, M0)
    temp = similar(Y)
    BQ = Matrix{eltype(Y)}(undef, M0, M0)

    # convercenve history
    neighistory = zeros(Int, maxiter)
    errorhistory = zeros(Float64, maxiter)

    # error estimation
    previous = 0.0
    rerror = 0.0
    iter = 1

    # evaluate and store necessary factorizations
    F = [lu(z*I - A) for z in zs]

    while true

        # apply rational filter to Y basis
        fill!(Q, zero(eltype(Q)))

        for i in 1:length(zs)
            for n in axes(Y, 2)

                ldiv!(q, F[i], view(Y, :, n))

                Q[:, n] .+= ws[i]*q

            end
        end

        # form effective matrices in reduced subspace
        mul3mat!(AQ, Q, A, Q, temp)
        mul!(BQ, Q', Q)

        # diagonalize subspace problem
        ritzvals, ritzvecs = eigen!(Hermitian(AQ), Hermitian(BQ), sortby = λ -> (!(λmin*(1+tol) < real(λ) <= λmax*(1+tol)), real(λ), imag(λ)))

        # update eigenvectors
        mul!(Y, Q, ritzvecs)

        # determine number of eigenvalues is inside
        neig = 0
        for λ in ritzvals
            if λmin < real(λ) <= λmax
                neig += 1
            end
        end
        neighistory[iter] = neig

        # estimate relative error of the trace of the eigenvalues
        tr = sum(ritzvals[1:neig])
        rerror = abs(tr - previous)/abs(tr)
        previous = sum(ritzvals[1:neig])
        errorhistory[iter] = rerror

        if rerror >= tol
            if iter < maxiter # if tolerance is not satisfied and number of iterations small than max, go to next iteration
                iter += 1
            else # if tolerance is not satisfied but max number of iterations is exceded, throw warning and return result
                @warn "maximum number of iterations, $maxiter, reached"
                return ritzvals, Y, neighistory, errorhistory
            end
        else # if tolerance is satisfied, return result
            return ritzvals, Y, neighistory, errorhistory
        end

        end # end while

end

function mul3mat!(Out, A, B, C, temp)
    mul!(temp, B, C)
    mul!(Out, A', temp)
    return Out
end

"""
    cauchynodes(λmin, λmax, cpts)

Builds the nodes and weights for a `cpts`-points Gauss-Legendre quadrature rule along a circle in the complex plane which includes `]λmin, λmax]`.
"""
function cauchynodes(λmin, λmax, cpts)

    xs, ws = gausslegendre(cpts)

    c = (λmax + λmin)/2
    R = (λmax - λmin)/2

    weights = zeros(Complex{Float64}, 2*cpts)
    points = zeros(Complex{Float64}, 2*cpts)

    for k in 1:cpts
        weights[k] = ws[k]*R*cis(pi*(xs[k]+1)/2)/4
        points[k] = c + R*cis(pi*(xs[k]+1)/2)

        weights[k + cpts] = ws[k]*R*cis(-pi*(xs[k]+1)/2)/4
        points[k + cpts] = c + R*cis(-pi*(xs[k]+1)/2)
    end

    return points, weights
end

end


text:
include("imports.jl")

# Path to scratch
RESULTS_DIR = "results"

# CB data path
CB_PATH = "../data/cb_transformed_reduced.csv"
SIM_SMALL_PATH = "../data/simdat-nfac500.bson"
SIM_LARGE_PATH = "../data/simdat-nfac5000.bson"

# Read CB data
function load_cb_data(path_to_data)
  function separatesamples(Y)
    sample_ids = sort(unique(Y.sample_id))
    return Matrix.(Y[Y.sample_id .== i, Not(:sample_id)]
                   for i in sample_ids)
  end

  Y = coalesce.(CSV.read(path_to_data), NaN)
  separatesamples(Y), names(Y)
end

load_simdata(path_to_data) = BSON.load(path_to_data)[:y]

function flowsanitize(Y)
  @rput Y
  return R"""
  colnames(Y) = 1:NCOL(Y)
  flowCore::flowFrame(Y)
  """
end

"""
replace missing values (NaN) in yi with x
"""
function replaceMissing(yi, x)
  out = deepcopy(yi)
  out[isnan.(out)] .= x
  return out
end

function main(y, markernames, outpath, K=30)
  I = length(y)
  N = [size(yi, 1) for yi in y]
  Y = replaceMissing(vcat(y...), -6)
  J = size(Y, 2)
  ff_Y = flowsanitize(Y)

  println("Running FlowSOM ...")
  @time fsom = FlowSOM.FlowSOM(ff_Y,
                               colsToUse=1:J,  # columns to use
                               maxMeta=K,  # Meta clustering option
                               seed=42)  # Seed for reproducible results:

  idx_upper = cumsum(N)
  idx_lower = [1; idx_upper[1:end-1] .+ 1]
  idx = [idx_lower idx_upper]

  fsmeta = fsom[:metaclustering]
  fsclus = fsmeta[Int.(convert(Matrix, fsom[:FlowSOM][:map][:mapping])[:, 1])]

  mkpath(outpath)
  for i in 1:I
    # print output
    lami = fsclus[idx[i, 1]:idx[i, 2]]

    # Plot y
    Wi = zeros(K)
    Ni = length(lami)
    for lam_in in lami
      Wi[lam_in] += 1 / Ni
    end

    println("printing figure $i")
    plt.figure(figsize=(6, 6))
    PlotUtils.plot_yz.plot_y(y[i], Wi, lami, lw=3, vlim=(-4, 4),
                             cm=PlotUtils.blue2red.cm(9),
                             fs_xlab=PlotUtils.rcParams["font.size"],
                             fs_ylab=PlotUtils.rcParams["font.size"],
                             fs_lab=PlotUtils.rcParams["font.size"],
                             fs_cbar=PlotUtils.rcParams["font.size"],
                             markernames=markernames)
    plt.savefig("$(outpath)/y$(i).pdf", bbox_inches="tight")
    plt.close()

    # Write clusterings to disk
    writedlm("$(outpath)/lam$(i).txt", lami)

    # Write cluster sizes to disk
    writedlm("$(outpath)/W$(i).txt", Wi)
  end
end

### Main ###

# CB FlowSOM
println("Analyze CB ...")
y_cb, markernames_cb = load_cb_data(CB_PATH)
main(y_cb, markernames_cb, "results/cb")

# Sim Small FlowSOM
println("Analyze sim small ...")
main(load_simdata(SIM_SMALL_PATH), [], "results/sim-small")

# Sim Large FlowSOM
println("Analyze sim large ...")
main(load_simdata(SIM_LARGE_PATH), [], "results/sim-large")

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
struct SupervisedAlignment
    nn_stability::Int
    C_steps::Int
    balanced::Bool
    adjusted::Bool
    max_gamma::Float64
end

function SupervisedAlignment(;nn_stability, C_steps, balanced, adjusted, max_gamma)
    SupervisedAlignment(nn_stability, C_steps, balanced, adjusted, max_gamma)
end

function Base.string(m::SupervisedAlignment)
    return "SupervisedAlignment_adj=$(m.adjusted)"
end

function find_parameters(method::SupervisedAlignment, data, labels, oracle, solver, qm, seed)
    nnpool = NNPool(data, k=1)

    for i in 1:size(data, 2)
        add_to_pools!(nnpool, i, ask_oracle(oracle, i), method.nn_stability)
    end
    gamma = find_gamma(nnpool, data, method.nn_stability; max_γ=method.max_gamma, adjusted=method.adjusted)

    @show "Found gamma $gamma"
    max_C = find_max_feasible_C(data, gamma, solver)
    min_C = find_min_feasible_C(data, gamma, solver)

    C_range, model_range = calculate_C_grid(data, labels, oracle, gamma, min_C, max_C, method.C_steps, solver)
    qm = calculate_qm(nnpool, data, labels, model_range, qm)
    qm_best_idx = findmax(qm)[2]
    C = C_range[qm_best_idx]
    return (gamma, C)
end

max_stars_repo_name: holtri/lama

max_issues_repo_name: holtri/lama

max_forks_repo_name: holtri/lama


text:
include("next_permutation.jl")

function compute()

	a=collect(0:9)

	const N=999999
	for i=1:N
		if !next_permutation(a)
			return
		end
	end
    ans=0
    for c in a
       ans =ans *10 +c
    end

    return string(ans)

end

println(compute())


text:
using ArgParse
using Formatting
using JSON,TOML
using DataStructures

include("runCmd.jl")

s = ArgParseSettings()
@add_arg_table s begin

    "--casename"
        help = "Prefix of the file."
        arg_type = String
        required = true

    "--input-dir"
        help = "The folder where case folders are contained"
        arg_type = String
        required = true

    "--output-dir"
        help = "Output folder"
        arg_type = String
        required = true

    "--year-rng"
        help = "Range of years being diagnosed"
        arg_type = Int64
        nargs = 2
        required = true

    "--domain-file"
        help = "Domain file of sea ice model"
        arg_type = String
        required = true

    "--overwrite"
        action = :store_true

    "--remap-file-nn"
        help = "Remap weighting file"
        arg_type = String
        required = true

end

parsed = DataStructures.OrderedDict(parse_args(ARGS, s))

JSON.print(parsed, 4)

for y = parsed["year-rng"][1]:parsed["year-rng"][2], m=1:12

    date_str = format("{:04d}-{:02d}", y, m)

    old_file  = joinpath(parsed["input-dir"],  "$(parsed["casename"]).cice.h.$(date_str).nc")
    new_file1 = joinpath(parsed["output-dir"], "$(parsed["casename"]).cice_extra1.$(date_str).nc")
    new_file2 = joinpath(parsed["output-dir"], "$(parsed["casename"]).cice_extra2_rg.$(date_str).nc")

    if !isfile(new_file1) || !isfile(new_file2) || parsed["overwrite"]

        println("Generating file: $(new_file1)")
        pleaseRun(`ncap2 -h -O -v
                    -s 'vice=vicen001+vicen002+vicen003+vicen004+vicen005;'
                    -s 'aice=(aicen001+aicen002+aicen003+aicen004+aicen005)/100.0;'
                    $old_file $new_file1
        `)

        println("Generating file: $(new_file2)")
        pleaseRun(`ncremap -v aice,vice -m $(parsed["remap-file-nn"]) $new_file1 $new_file2`)
    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
immutable SkierUnif
    time::Int
    xlim::NTuple{2, Float64}
    ylim::NTuple{2, Float64}
end

eltype(::Type{SkierUnif}) = SkierState

function rand(rng::AbstractRNG, d::SkierUnif)
    x = rand(rng)*(last(d.xlim)-first(d.xlim))+first(d.xlim)
    y = rand(rng)*(last(d.ylim)-first(d.ylim))+first(d.ylim)
    psi = rand(rng)*2*pi
    return SkierState(d.time, [x,y], psi)
end

initial_state_distribution(p::PowseekerProblem) = SkierUnif(1, mdp(p).xlim, mdp(p).ylim)


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

mutable struct VirtualMachineScaleSetSku <: SwaggerModel
    resourceType::Any # spec type: Union{ Nothing, String } # spec name: resourceType
    sku::Any # spec type: Union{ Nothing, Sku } # spec name: sku
    capacity::Any # spec type: Union{ Nothing, VirtualMachineScaleSetSkuCapacity } # spec name: capacity

    function VirtualMachineScaleSetSku(;resourceType=nothing, sku=nothing, capacity=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetSku, Symbol("resourceType"), resourceType)
        setfield!(o, Symbol("resourceType"), resourceType)
        validate_property(VirtualMachineScaleSetSku, Symbol("sku"), sku)
        setfield!(o, Symbol("sku"), sku)
        validate_property(VirtualMachineScaleSetSku, Symbol("capacity"), capacity)
        setfield!(o, Symbol("capacity"), capacity)
        o
    end
end # type VirtualMachineScaleSetSku

const _property_map_VirtualMachineScaleSetSku = Dict{Symbol,Symbol}(Symbol("resourceType")=>Symbol("resourceType"), Symbol("sku")=>Symbol("sku"), Symbol("capacity")=>Symbol("capacity"))
const _property_types_VirtualMachineScaleSetSku = Dict{Symbol,String}(Symbol("resourceType")=>"String", Symbol("sku")=>"Sku", Symbol("capacity")=>"VirtualMachineScaleSetSkuCapacity")
Base.propertynames(::Type{ VirtualMachineScaleSetSku }) = collect(keys(_property_map_VirtualMachineScaleSetSku))
Swagger.property_type(::Type{ VirtualMachineScaleSetSku }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineScaleSetSku[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetSku }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetSku[property_name]

function check_required(o::VirtualMachineScaleSetSku)
    true
end

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


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	emissions(EP::Model, inputs::Dict, UCommit::Int)

This function creates expression to add the CO2 emissions by plants in each zone, which is subsequently added to the total emissions
"""
function emissions(EP::Model, inputs::Dict)

	println("Emissions Module (for CO2 Policy modularization")

	dfGen = inputs["dfGen"]

	G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
	T = inputs["T"]     # Number of time steps (hours)
	Z = inputs["Z"]     # Number of zones
	COMMIT = inputs["COMMIT"] # For not, thermal resources are the only ones eligible for Unit Committment

	@expression(EP, eEmissionsByPlant[y=1:G,t=1:T],
	 	if y in inputs["COMMIT"]
		 	dfGen[!,:CO2_per_MWh][y]*EP[:vP][y,t]+dfGen[!,:CO2_per_Start][y]*EP[:vSTART][y,t]
	 	else
		 	dfGen[!,:CO2_per_MWh][y]*EP[:vP][y,t]
	 	end
 	)

 	@expression(EP, eEmissionsByZone[z=1:Z, t=1:T], sum(eEmissionsByPlant[y,t] for y in dfGen[(dfGen[!,:Zone].==z),:R_ID]))

	return EP
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
using KernelAbstractions, CUDAKernels, CUDA, Test
using KernelAbstractions.Extras: @unroll

has_cuda_gpu() || exit()
CUDA.allowscalar(false)

const nreps = 3
const N = 2048
const T = Float32

const TILE_DIM = 32
const BLOCK_ROWS = 8

# Simple variants

@kernel function simple_copy_kernel!(output, @Const(input))
  I, J = @index(Global, NTuple)
  @inbounds output[I, J] = input[I, J]
end

@kernel function simple_transpose_kernel!(output, @Const(input))
    I, J = @index(Global, NTuple)
    @inbounds output[I, J] = input[I, J]
end

# Local memory variants

@kernel function lmem_copy_kernel!(output, @Const(input),
                                   ::Val{BANK}=Val(1)) where BANK
    I, J = @index(Global, NTuple)
    i, j = @index(Local,  NTuple)

    N = @uniform groupsize()[1]
    M = @uniform groupsize()[2]

    # +1 to avoid bank conflicts on shared memory
    tile = @localmem eltype(output) (N+BANK, M)

    @inbounds tile[i, j] = input[I, J]

    @synchronize

    @inbounds output[I, J] = tile[i, j]
end

@kernel function lmem_transpose_kernel!(output, @Const(input),
                                       ::Val{BANK}=Val(1)) where BANK
    gi, gj = @index(Group, NTuple)
    i, j = @index(Local,  NTuple)

    N = @uniform groupsize()[1]
    M = @uniform groupsize()[2]

    # +1 to avoid bank conflicts on shared memory
    tile = @localmem eltype(output) (N+BANK, M)

    # Manually calculate global indexes
    # Later on we need to pivot the group index
    I = (gi-1) * N + i
    J = (gj-1) * M + j

    @inbounds tile[i, j] = input[I, J]

    @synchronize

    # Pivot the group index
    I = (gj-1) * M + i
    J = (gi-1) * N + j

    @inbounds output[I, J] = tile[j, i]
end

# Local Memory + process multiple elements per lane

@kernel function coalesced_copy_kernel!(output, @Const(input),
                                        ::Val{BANK}=Val(1)) where BANK
    gi, gj = @index(Group, NTuple)
    i, j   = @index(Local, NTuple)

    TILE_DIM   = @uniform groupsize()[1]
    BLOCK_ROWS = @uniform groupsize()[2]

    # +1 to avoid bank conflicts on shared memory
    tile = @localmem eltype(output) (TILE_DIM+BANK, TILE_DIM)

    # Can't use @index(Global), because we use a smaller ndrange
    I = (gi-1) * TILE_DIM + i
    J = (gj-1) * TILE_DIM + j

    @unroll for k in 0:BLOCK_ROWS:(TILE_DIM-1)
        @inbounds tile[i, j+k] = input[I, J+k]
    end

    @synchronize

    @unroll for k in 0:BLOCK_ROWS:(TILE_DIM-1)
        @inbounds output[I, J+k] = tile[i, j+k]
    end
end

@kernel function coalesced_transpose_kernel!(output, @Const(input),
                                             ::Val{BANK}=Val(1)) where BANK
    gi, gj = @index(Group, NTuple)
    i, j   = @index(Local, NTuple)

    TILE_DIM   = @uniform groupsize()[1]
    BLOCK_ROWS = @uniform groupsize()[2]

    # +1 to avoid bank conflicts on shared memory
    tile = @localmem eltype(output) (TILE_DIM+BANK, TILE_DIM)

    # Can't use @index(Global), because we use a smaller ndrange
    I = (gi-1) * TILE_DIM + i
    J = (gj-1) * TILE_DIM + j

    @unroll for k in 0:BLOCK_ROWS:(TILE_DIM-1)
        @inbounds tile[i, j+k] = input[I, J+k]
    end

    @synchronize

    # Transpose block offsets
    I = (gj-1) * TILE_DIM + i
    J = (gi-1) * TILE_DIM + j

    @unroll for k in 0:BLOCK_ROWS:(TILE_DIM-1)
        @inbounds output[I, J+k] = tile[j+k, i]
    end
end

# Benchmark simple

for block_dims in ((TILE_DIM, TILE_DIM), (TILE_DIM*TILE_DIM, 1), (1, TILE_DIM*TILE_DIM))
    for (name, kernel) in (
                            ("copy",      simple_copy_kernel!(CUDADevice(), block_dims)),
                            ("transpose", simple_transpose_kernel!(CUDADevice(), block_dims)),
                          )
        NVTX.@range "Simple $name $block_dims" let
            input = CUDA.rand(T, (N, N))
            output = similar(input)

            # compile kernel
            ev = kernel(input, output, ndrange=size(output))
            CUDA.@profile begin
                for rep in 1:nreps
                  ev = kernel(input, output, ndrange=size(output), dependencies=(ev,))
                end
                wait(ev)
            end
        end
    end
end

# Benchmark localmem
for (name, kernel) in (
                        ("copy",      lmem_copy_kernel!(CUDADevice(), (TILE_DIM, TILE_DIM))),
                        ("transpose", lmem_transpose_kernel!(CUDADevice(), (TILE_DIM, TILE_DIM))),
                      )
    for bank in (true, false)
        NVTX.@range "Localmem $name ($TILE_DIM, $TILE_DIM) bank=$bank" let
            input = CUDA.rand(T, (N, N))
            output = similar(input)

            # compile kernel
            ev = kernel(input, output, Val(Int(bank)), ndrange=size(output))
            CUDA.@profile begin
                for rep in 1:nreps
                    ev = kernel(input, output, Val(Int(bank)), ndrange=size(output), dependencies=(ev,))
                end
                wait(ev)
            end
        end
    end
end

# Benchmark localmem + multiple elements per lane
for (name, kernel) in (
                        ("copy",      coalesced_copy_kernel!(CUDADevice(), (TILE_DIM, BLOCK_ROWS))),
                        ("transpose", coalesced_transpose_kernel!(CUDADevice(), (TILE_DIM, BLOCK_ROWS))),
                      )
    for bank in (true, false)
        NVTX.@range "Localmem + multiple elements $name ($TILE_DIM, $BLOCK_ROWS) bank=$bank" let
            input = CUDA.rand(T, (N, N))
            output = similar(input)

            # We want a number of blocks equivalent to (TILE_DIM, TILE_DIM)
            # but our blocks are (TILE_DIM, BLOCK_ROWS) so we need to remove
            # a factor from the size of the array otherwise we get to many blocks
            block_factor = div(TILE_DIM, BLOCK_ROWS)
            ndrange = (N, div(N, block_factor))

            # compile kernel
            ev = kernel(input, output, Val(Int(bank)), ndrange=ndrange)
            CUDA.@profile begin
                for rep in 1:nreps
                    ev = kernel(input, output, Val(Int(bank)), ndrange=ndrange, dependencies=(ev,))
                end
                wait(ev)
            end
        end
    end
end


text:
struct FirstOrderSommerfeld <: AbstractBoundaryCondition end
struct Higdon <: AbstractBoundaryCondition end
struct BetzMittra <: AbstractBoundaryCondition end


text:
using ATA
# For load data to input in custom obj_fun:
using CSV
using FileIO
using LinearAlgebra
cd("where your input files are")

# 1. Start ATA and add file with custom settings (Needed)
ata_model = start_ata(;
    settings_file = "settings_ata_custom.jl",
    bank_file = "data/bank.csv",
    bank_delim = ";",
)

# 2. Add categorical constraints (Optional)
add_constraints!(ata_model; constraints_file = "constraints.csv", constraints_delim = ";")

# 3. Add overlap maxima (Optional)
add_overlap!(ata_model; overlap_file = "overlap_matrix.csv", overlap_delim = ";")

# custom objective type, function and arguments
ata_model.obj.name = "custom"

ata_model.obj.fun = function (x::Matrix{Float64}, obj_args::NamedTuple)
    IIF = obj_args.IIF
    T = obj_args.T
    TIF = zeros(Float64, T)
    for t = 1:T
        K, I = size(IIF[t])
        # ungroup items
        xₜ = fs_to_items(x[:, t], I, obj_args.fs_items)
        if K > 1
            TIF[t] = Inf
            for k = 1:K
                TIF[t] = min(TIF[t], LinearAlgebra.dot(IIF[t][k, :], xₜ))
            end
        else
            TIF[t] = LinearAlgebra.dot(IIF[t][1, :], xₜ)[1]
        end
    end
    min_TIF = minimum(TIF)
    TIF = [min_TIF for t = 1:T]
    # Must return a vector of length T.
    # The resulting objective function is the minimum of all values in this vector.
    return TIF::Vector{Float64}
end

ata_model.obj.args = (
    T = ata_model.settings.T,
    IIF = FileIO.load("data/IIF.jld2", "IIF"),
    fs_items = ata_model.settings.fs.items,
)

print_infos(ata_model)
# Assembly settings

# SIMAN (Suggested for Large scale ATA):
# Set the solver, "siman" for simulated annealing, "jump" for MILP solver.
solver = "siman"

# SIMAN (Suggested for Large scale ATA):

start_temp = 0.0001
# Default: `0.1`. Values:  `[0, Inf]`.
# Starting temperature, set to minimum for short journeys (if 0 worse solutions will never be accepted).

geom_temp = 0.1
# Default: `0.1`. Values:  `[0, Inf)`.
# Decreasing geometric factor.

n_item_sample = Inf
# Default: 1. Values: `[1, Inf]`.
# Number of items to alter. Set to minimum for a shallow analysis,
# set to maximum for a deep analysis of the neighbourhoods.

n_test_sample = ata_model.settings.T
# Default: 1. Values: `[1, Inf]`.
# Number of tests to alter. Set to minimum for a shallow analysis, set to maximum for a deep analysis of the neighbourhoods.

opt_feas = 0.9
# Default: 0.0. Values: `[0, Inf)`.
# Optimality/feasibility balancer, if 0 only feasibility of solution is analysed. Viceversa, if 1, only optimality is considered (uncontrained model). All the other values are accepted but produce uninterpretable results.

n_fill = 1
# Default: 1. Values: `[0, Inf)`.
# Number of fill-up phases, usually 1 is sufficient, if start_temp is high it can be higher.
# If a starting_design is supplied, it should be set to 0.

verbosity = 1
# Default: 2. Values: `1` (minimal), `2` (detailed).
# Verbosity level. In the console '+' stands for improvement, '_' for accepting worse solution.
# The dots are the fill-up improvement steps.

#! Termination criteria:

max_time = 100.0
# Default: `1000.0`. Values: `[0, Inf)`.
# Time limit in seconds.

max_conv = 5
# Default: `2`. Values: `[1, Inf)`.
# Maximum convergence, stop when, after max_conv rounds no improvements have been found.
# Set to minimum for shallow analysis, increase it for deep analysis of neighbourhoods.

feas_nh = 1
# Default: `0`. Values: `[1, Inf)`.
# Maximum number of Feasibility neighbourhoods to explore, set to the minimum if the model is small or not highly constrained.

opt_nh = Inf
# Default: `5`. Values: `[1, Inf)`.
# Maximum number of Optimality neighbourhoods to explore, set to the minimum if the model is highly constrained.

# 4. assemble
assemble!(
    ata_model;
    solver = solver,
    max_time = max_time,
    start_temp = start_temp,
    geom_temp = geom_temp,
    n_item_sample = n_item_sample,
    n_test_sample = n_test_sample,
    verbosity = verbosity,
    max_conv = max_conv,
    opt_feas = opt_feas,
    n_fill = n_fill,
    feas_nh = feas_nh,
    opt_nh = opt_nh,# ,
    # optimizer_attributes = optimizer_constructor,
    # optimizer_constructor =optimizer_attributes
)

# All the settings and outputs from optimization are in ata_model object.
# See the struct in ATA.jl to understand how to retrieve all the information.
# A summary of the resulting tests is available in results_folder/results.txt
# If siman is chosen, the optimality and feasibility of the best neighbourhood
# is reported in "results/results_ata.jl"

ata_model.obj.name = "maximin"
print_results(ata_model; results_folder = "results")

#]add https://github.com/giadasp/ATAPlot.jl
using ATAPlot
plot_results(ata_model; plots_folder = "plots")


text:
## Aircraft stability analysis example
using AeroMDAO
using LinearAlgebra # For norm()

# Lifting surfaces
wing  = WingSection(span       = 4.0,
                    dihedral   = 5.0,
                    sweep_LE   = 15.0,
                    taper      = 0.4,
                    root_chord = 2.0,
                    root_twist = 0.0,
                    tip_twist  = -2.0,
                    root_foil  = naca4((2,4,1,2)),
                    tip_foil   = naca4((2,4,1,2)))
wing_mac  = mean_aerodynamic_center(wing)
wing_pos  = [0., 0., 0.]
wing_plan = plot_wing(wing;
                      position = wing_pos)

print_info(wing, "Wing")

htail = WingSection(span       = 1.0,
                    dihedral   = 0.0,
                    sweep_LE   = 15.0,
                    taper      = 0.6,
                    root_chord = 0.8,
                    root_twist = 0.0,
                    tip_twist  = 0.0,
                    root_foil  = naca4((0,0,1,2)),
                    tip_foil   = naca4((0,0,0,9)));
htail_mac  = mean_aerodynamic_center(htail)
htail_pos  = [5., 0., 0.]
α_h_i      = 0.
htail_plan = plot_wing(htail;
                       position = htail_pos)

print_info(htail, "Horizontal Tail")

vtail = HalfWingSection(span       = 0.8,
                        dihedral   = 0.0,
                        sweep_LE   = 8.0,
                        taper      = 0.6,
                        root_chord = 0.8,
                        root_twist = 0.0,
                        tip_twist  = 0.,
                        root_foil  = naca4((0,0,0,9)),
                        tip_foil   = naca4((0,0,0,9)))
vtail_mac  = mean_aerodynamic_center(vtail) # NEEDS FIXING FOR ROTATION
vtail_pos  = [5., 0., 0.]
vtail_plan = plot_wing(vtail;
                       position = vtail_pos,
                       angle    = π/2)

print_info(vtail, "Vertical Tail")

## Static stability
S, b, c = projected_area(wing), span(wing), mean_aerodynamic_chord(wing)
x_w = wing_pos + [ wing_mac[1], 0, 0 ]

# Horizontal tail volume coefficient
S_h = projected_area(htail)
x_h = htail_pos + [ htail_mac[1], 0, 0 ]
l_h = norm(x_h - x_w)
V_h = S_h * l_h / (S * c)
println("Horizontal TVC      V_h: $V_h")

# Vertical tail volume coefficient
S_v = projected_area(vtail)
x_v = vtail_pos + [ vtail_mac[1], 0, 0 ]
l_v = norm(x_v - x_w)
V_v = S_v * l_v / (S * b)
println("Vertical TVC        V_v: $V_v")

## Panelling and assembly
wing_panels  = 	panel_wing(wing, [20], 10;
                           position = wing_pos
                          )
htail_panels =	panel_wing(htail, [10], 5;
                           position = htail_pos,
                           angle    = deg2rad(α_h_i),
                           axis     = [0., 1., 0.]
                          )
vtail_panels = 	panel_wing(vtail, [10], 5;
                           position = vtail_pos,
                           angle    = π/2
                          )

aircraft = Dict("Wing"            => wing_panels,
                "Horizontal Tail" => htail_panels,
                "Vertical Tail"   => vtail_panels);

## Evaluate case
ac_name = "My Aircraft"
ρ 		= 1.225
ref     = x_w
V, α, β = 1.0, 3.0, 0.0
Ω 		= [0.0, 0.0, 0.0]
fs 	    = Freestream(V, α, β, Ω)

@time dv_data =
solve_stability_case(aircraft, fs;
                     rho_ref     = ρ,           # Reference density
                     r_ref       = ref,         # Reference point for moments
                     area_ref    = S,           # Reference area
                     span_ref    = b,           # Reference span
                     chord_ref   = c,           # Reference chord
                     name        = ac_name,     # Aircraft name
                     print       = true,        # Prints the results for the entire aircraft
                     print_components = true,   # Prints the results for each component
                    );

## Process data
labels = (collect ∘ keys)(dv_data) # Get aircraft component names from analysis
comp = labels[1]                   # Pick your component
nf, ff, dvs = dv_data[comp];       # Get the nearfield, farfield, and stablity derivative coefficients
print_case(dv_data, comp)          # Pretty-print the results

## Aerodynamic quantities of aircraft
nf_plane, ff_plane, dvs_plane = dv_data[labels[1]]

# Center of pressure
Cm   = nf_plane[5]	# Moment coefficient
CL   = nf_plane[3]	# Lift coefficient
cp   = -c * Cm / CL

# Neutral point
Cm_α = dvs_plane[5,1]
CL_α = dvs_plane[3,1]
np   = -c * Cm_α / CL_α

# Spiral stability
Cl_β = dvs_plane[4,2]
Cl_r = dvs_plane[4,5]
Cn_r = dvs_plane[6,5]
Cn_β = dvs_plane[6,2]

γ = Cl_β * Cn_r / (Cl_r * Cn_β) # Check degree-radian issue

# Locations
x_np = [ ref[1] .+ np; zeros(2) ]	# Neutral point
x_cp = [ ref[1] .+ cp; zeros(2) ]	# Center of pressure

println("Aerodynamic Center x_ac: $(x_w[1]) m")
println("Neutral Point      x_np: $(x_np[1]) m")
println("Center of Pressure x_cp: $(x_cp[1]) m")
println("Spiral Stability      γ: $γ")

## Plotting everything
using Plots
gr(dpi = 300)

#
z_limit = b
aircraft_plot =
plot(xaxis = "x", yaxis = "y", zaxis = "z",
     aspect_ratio = 1,
     camera = (30, 60),
     xlim = (0, z_limit),
    #  ylim = (-z_limit/2, z_limit/2),
     zlim = (-z_limit/2, z_limit/2),
     size = (1280, 720)
    )

plot!(wing_plan, label = "Wing")
plot!(htail_plan, label = "Horizontal Tail")
plot!(vtail_plan, label = "Vertical Tail")

scatter!(tuple(x_w...), color = :blue, label = "Wing MAC")
scatter!(tuple(x_h...), color = :red, label = "Horizontal Tail MAC")
scatter!(tuple(x_v...), color = :green, label = "Vertical Tail MAC")

scatter!(tuple(x_np...), color = :orange, label = "Neutral Point")
scatter!(tuple(x_cp...), color = :brown, label = "Center of Pressure")
# savefig(aircraft_plot, "Aircraft.png")
plot!()

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
using ExtensibleUnions
using Test


text:
using Pkg
using HTTP
using URIs

export webactivate

"""
    webactivate(workspace::AbstractSTring)
    webactivate(workspace::URI)
Activate the workspace whose Project and Manifest files are
located at the specified URI.
This is done by creating a local temporary directory, copying
the files there and activating it.
The path to the temporary directory is returned.
"""
function webactivate end

function webactivate(workspace::AbstractString)
    webactivate(URI(workspace))
end

function webactivate(workspace::URI)
    local_workspace = mktempdir()
    for f in ["Project.toml", "Manifest.toml"]
	from = resolvereference(workspace, f)
	to = joinpath(local_workspace, f)
	response = HTTP.request("GET", from)
	@assert response.status == 200
	open(to, "w") do f
	    write(f, String(response.body))
	end
    end
    Pkg.activate(local_workspace)
    local_workspace
end


text:
# These functions convert 5-tuples to/from Int16s

Byte = UInt8

"""
    code_2_byte(code::NTuple{5,Int})
Compress a Worlde result into a single byte.
"""
code_2_byte(code::NTuple{5,Int})::Byte = evalpoly(3, code)

"""
    byte_2_code(x)
Convert a (byte) into a 5-tuple, inverting the action of `code_2_byte`.
"""
byte_2_code(x::Integer) = Tuple(digits(x, base = 3, pad = five))

export all_pairs_compute, code_2_byte, byte_2_code

"""
    all_pairs_compute
Precompute a (compressed) Wordle score for all pairs of words.
"""
function all_pairs_compute()::Matrix{Byte}

    n_ans = length(ANS_LIST)
    n_guess = length(GUESS_LIST)
    M = zeros(Byte, n_guess, n_ans)

    PM = Progress(n_guess)
    for i = 1:n_guess
        a = GUESS_LIST[i]
        for j = 1:n_ans
            b = ANS_LIST[j]
            M[i, j] = code_2_byte(wordle_score(a, b))
        end
        next!(PM)
    end
    return M
end


text:
# Seedlink with command-line stations
config_file = path*"/SampleFiles/seedlink.conf"
st = 0.0
en = 60.0
seq = "500000"
sta = ["CC.SEP", "UW.HDW"]
pat = ["?????.D", "?????.D"]
trl = ".??.???.D"
sta_matrix = String.(reshape(split(sta[2],'.'), 1,2))

printstyled("  SeedLink\n", color=:light_green)
printstyled("  (SeedLink tests require up to 6 minutes)\n", color=:green)

# has_stream
printstyled("    has_stream\n", color=:light_green)
ta = Array{Union{Task,Nothing}, 1}(undef, 4)
ta[1] = @async has_stream(sta)[2]
ta[2] = @async has_stream(sta, pat, d='.')[2]
ta[3] = @async has_stream(join(sta, ','))[2]
ta[4] = @async has_stream(sta_matrix)[1]
SL_wait(ta, 1)

# has_stream
printstyled("    has_sta\n", color=:light_green)
ta[1] = @async has_sta(sta[1])[1]
ta[2] = @async has_sta(sta[1]*trl)[1]
ta[3] = @async has_sta(sta)[1]
ta[4] = @async has_sta(parse_charr(sta, '.', false))[1]
SL_wait(ta, 1)

# Attempting to produce errors
printstyled("    produce expected errors and warnings\n", color=:light_green)
redirect_stdout(out) do
  S1 = SeisData()

  @test_throws ErrorException seedlink!(S1, "DATA", [sta[1]], ["*****.X"])

  S2 = seedlink("DATA", [sta[1]], pat, x_on_err=false)
  write(S2.c[1], "BYE\r")
  close(S2.c[1])
  @test_throws ErrorException seedlink!(S2, "DATA", [replace(sta[1], "SEP" => "XOX")], ["?????.D"])

  S3 = seedlink("DATA", [replace(sta[1], "SEP" => "XOX")], ["*****.X"], x_on_err=false)
  write(S3.c[1], "BYE\r")
  close(S3.c[1])

  S4 = seedlink("DATA", hcat(sta_matrix, "***", "***", "X"), x_on_err=false)
  write(S4.c[1], "BYE\r")
  close(S4.c[1])
end

# DATA mode
printstyled("    DATA mode\n", color=:light_green)
printstyled("      link 1: command-line station list\n", color=:light_green)
T = SeisData()
redirect_stdout(out) do
  seedlink!(T, "DATA", sta, refresh=9.9, kai=7.0, v=1)
end

printstyled("      link 2: station file\n", color=:light_green)
redirect_stdout(out) do
  seedlink!(T, "DATA", config_file, refresh=13.3, v=3)
end
wait_on_data!(T, 60.0)

# To ensure precise timing, we'll pass d0 and d1 as strings
dt = en-st
(d0,d1) = parsetimewin(st,en)

# FETCH mode (indistinguishable from DATA mode for most users)
printstyled("    FETCH mode with seq (should fail)\n", color=:light_green)
redirect_stdout(out) do
  V = seedlink("FETCH", refresh=10.0, config_file, seq=seq, s=now()-Hour(1), v=2)
  wait_on_data!(V, 30.0)
end

printstyled("    FETCH mode\n", color=:light_green)
V = seedlink("FETCH", refresh=10.0, config_file)
printstyled("      link initialized\n", color=:light_green)
wait_on_data!(V, 30.0)

# SeedLink time mode (more complicated)
printstyled("    TIME mode\n", color=:light_green)
U = SeisData()
seedlink!(U, "TIME", sta, refresh=10.0, s=d0, t=d1, w=true)
printstyled("      first link initialized\n", color=:light_green)

# Seedlink with a config file
seedlink!(U, "TIME", config_file, refresh=10.0, s=d0, t=d1)
printstyled("      second link initialized\n", color=:light_green)

# Seedlink with a config string
redirect_stdout(out) do
  seedlink!(U, "TIME", "CC.VALT..???, UW.ELK..EHZ", refresh=10.0, s=d0, t=d1, v=3)
end
printstyled("      third link initialized\n", color=:light_green)
wait_on_data!(U, 60.0)


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id: HairpinIncr.mso								$
#*------------------------------------------------------------------
type Summary_Hairpin
	Summary_Hairpin()=begin
		new(
			area(Dict{Symbol,Any}(
				:Brief=>"Total Exchange Surface Area"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Total Duty",
				:Default=>7000,
				:Lower=>1e-8,
				:Upper=>1e10
			)),
			Results_Hairpin(Dict{Symbol,Any}(
				:Brief=>"Inner Side Summary"
			)),
			Results_Hairpin(Dict{Symbol,Any}(
				:Brief=>"Outer Side Summary"
			)),
			[:A,:Qtotal,:Inner,:Outer,]
		)
	end
	A::area
	Qtotal::power
	Inner::Results_Hairpin
	Outer::Results_Hairpin
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Summary_Hairpin
function atributes(in::Summary_Hairpin,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="to be documented"
	fields[:Info]="to be documented"
	drive!(fields,_)
	return fields
end
Summary_Hairpin(_::Dict{Symbol,Any})=begin
	newModel=Summary_Hairpin()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
module TensorNetworkTensors
import TensorOperations, LinearAlgebra
using TupleTools
using Parameters

const TT = TupleTools
const TO = TensorOperations
const LA = LinearAlgebra

export gatherby, groupby
include("auxiliaryfunctions.jl")

export AbstractTensor
"""
    AbstractTensor{T,N}
Abstract supertype for all tensornetwork-tensors.
"""
abstract type AbstractTensor{T,N} end

export DTensor
include("DTensors.jl")

export DAS, DASCharges, DASCharge, NDAS, NDASCharges, NDASCharge,
        U1, U1Charges, U1Charge, ZN, ZNCharges, ZNCharge
export Z2, Z2Charges, Z2Charge
export chargeindex, ⊗, ⊕, DASSector, InOut, charge, deleteat, permute
export allsectors, covariantsectors, invariantsectors
export DASTensor
export charges, sizes, in_out, tensor, symmetry, setcharges!, setsizes!, setin_out!,
    settensor!, isinvariant, initwithrand!, initwithzero!, setcharge!
export chargestype, chargetype, dimstype, symmetry

include("DASTensors.jl")

export toarray, todense

#Multiplication and Division with Numbers
export apply, apply!
include("tensorops.jl")

#Overload functions from TensorOperations
include("tensoroperations.jl")

#Reshaping: Splitting and fusing legs of tensors
export  fuselegs, fuselegs!,
        splitlegs, splitlegs!
include("splitfuse.jl")

#Factorizations of Tensors
export tensorsvd, tensorsvd!
export tensorqr, tensorrq
export tensoreig
export  svdtrunc_default,
        svdtrunc_discardzero,
        svdtrunc_maxχ,
        svdtrunc_maxcumerror,
        svdtrunc_maxerror
include("tensorfactorization.jl")

#Overload the functions required for use with KrylovKit
include("krylovkit.jl")
end # module


text:
using Printf
using Random
using Distances
using Distributions
using LinearAlgebra

using TripletEmbeddings
using Test

@testset "embeddings.jl" begin
    d = 2
    n = 10
    X = Embedding(d,n)

    @test size(X) == (d,n)
    @test size(Embedding([1 0; 0 1])) == (2,2)
    @test size(Embedding((1,10))) == (1,10)
    @test eltype(Embedding([1 1; 1 1])) <: AbstractFloat # This is needed for kernel(loss, X)

    @test nitems(X) == n
    @test ndims(X) == d

    X₀ = rand(n)
    X = Embedding(X₀)
    @test nitems(X) == n
    @test ndims(X) == 1

    @test V(d) == [0.5 -0.5; -0.5 0.5]

    @test  L(3, (1,2,3)) == [0.0 -1.0 1.0; -1.0 1.0 0.0; 1.0 0.0 -1.0]
end

@testset "triplets.jl" begin
    data = [1 2 3]
    test_triplets = [Triplet(Int16, 3, 2, 1); Triplet(Int16, 1, 2, 3)]
    triplets = Triplets(data)
    @test triplets == test_triplets

    data = [1 2 3; 1 2 3]
    triplets = Triplets(data)
    @test triplets == test_triplets

    @test size(triplets) == (2,)
    @test length(triplets) == 2

    triplets = [Triplet(Int32, 1,2,3), Triplet(1,2,5)]
    @test checktriplets(triplets) == false

    @test_throws ArgumentError Triplets(ones(1,10))

    c = 0
    triplets = sampletriplet.(fill(1000, 1000000))
    for t in triplets
        c += t[1] == t[2] || t[1] == t[3]
    end
    @test c == 0
end

@testset "STE.jl" begin
    σ = 1/sqrt(2)
    loss = STE(σ = σ)

    @test loss.σ == σ
    @test loss.constant == 1/σ^2

    data = [1 2 3; 1 2 3]
    triplets = Triplets(data)
end

function test_insample(d::Int)
    @assert d == 1 || d == 2

    n = 100
    X = rand(d, n)
    if d == 1
        Y = randn() * X .+ 5 * rand()
    else
        θ = rand(Uniform(0, 2π))
        R = [cos(θ) -sin(θ); sin(θ) cos(θ)]
        Y = broadcast(+, randn() * R * X , 5 * rand(d))
    end
    X = Embedding(X)
    Z, transform = procrustes(X, Y)

    return mse(Z, Y)
end

function test_outofsample()
    Y = broadcast(+, [1 1 -1 -1; 1 -1 1 -1], [2, 2])
    s = rand()
    θ = rand()
    R = [cos(θ) -sin(θ); sin(θ) cos(θ)]
    b = rand(2)
    X = [Y mean(Y, dims=2)]
    X = broadcast(+, s * R * X, b)

    _, transform = procrustes(X[:,1:end-1], Y)
    return mean(Y, dims=2) ≈ apply(transform, X[:,end])
end

@testset "procrustes.jl" begin
    @test isapprox(test_insample(1), 0.0; atol=eps())
    @test isapprox(test_insample(2), 0.0; atol=eps())
    @test test_outofsample()
end


text:
using SpinGlassNetworks
using SpinGlassTensors
using SpinGlassEngine

function bench()
    m = 16
    n = 16
    t = 8

    β = 1.

    L = n * m * t
    num_states = 100

    instance = "$(@__DIR__)/instances/chimera_droplets/2048power/001.txt"

    ig = ising_graph(instance)

    fg = factor_graph(
        ig,
        #spectrum=full_spectrum,
        cluster_assignment_rule=super_square_lattice((m, n, t))
    )

    #for transform ∈ all_lattice_transformations
    peps = PEPSNetwork{false}(m, n, fg, rotation(0), β=β, bond_dim=32)
    update_gauges!(peps, :rand)
    @time sol = low_energy_spectrum(peps, num_states)#, merge_branches(peps, 1.0))
    println(sol.energies[1:1])
end

bench()


text:
export import_csv
export read_nnimage
export make_centroids
export centroid2xy
export import_stack
export stack2xy

function import_csv(centroid_data_file::String)
    nnResult = read(centroid_data_file; header=false)
end

function read_nnimage(image_path::String)
    gray_im = Gray.(load(image_path))
end

function make_centroids(gray_im; threshold=0.5)
    bw_im = gray_im .> threshold
    markers = label_components(erode(bw_im))
    centroids = component_centroids(markers)
    popfirst!(centroids)
    return centroids
end

function centroid2xy(centroids)
    xy = [first.(centroids) last.(centroids)]
    atom_xy = permutedims(xy, (2,1))
end

function import_stack(nnResult_stack_path::String; range=[])
    if range == []
        nnResult_stack = h5read(nnResult_stack_path, "Dataset1")
    else
        nnResult_stack = h5read(nnResult_stack_path, "Dataset1", (:,:,range))
    end
    return nnResult_stack
end

function stack2xy(nnResult_stack)
    centroids_collection = map(x->make_centroids(nnResult_stack[:,:,x]), 1:size(nnResult_stack,3))
    xy_stack = map(centroid2xy, centroids_collection)
    return xy_stack
end