text:
# Simple data provider that load text
using Iterators
using MXNet

function build_vocabulary(corpus_fn::AbstractString, vocab_fn::AbstractString; max_vocab=10000)
  if isfile(vocab_fn)
    info("Vocabulary already exists, reusing $vocab_fn...")
    vocab = Dict{Char,Int}([w => i for (i,w) in enumerate(readall(vocab_fn))])
  else
    # count symbol frequency
    dict = Dict{Char,Int}()
    open(corpus_fn) do io
      for line in eachline(io)
        for c in line
          dict[c] = get(dict, c, 0) + 1
        end
      end
    end

    vocab = sort(collect(dict), by=x->-x.second)
    vocab = vocab[1:min(max_vocab,length(vocab))]
    open(vocab_fn, "w") do io
      for x in vocab
        print(io, x.first)
      end
    end

    vocab = Dict([x.first => i for (i,x) in enumerate(vocab)])
  end
  vocab[UNKNOWN_CHAR] = length(vocab)
  return vocab
end

#--CharSeqProvider
type CharSeqProvider <: mx.AbstractDataProvider
  text       :: AbstractString
  batch_size :: Int
  seq_len    :: Int
  vocab      :: Dict{Char,Int}

  prefix     :: Symbol
  n_layer    :: Int
  dim_hidden :: Int
end
#--/CharSeqProvider

function mx.get_batch_size(p :: CharSeqProvider)
  p.batch_size
end

#--provide
function mx.provide_data(p :: CharSeqProvider)
  [(symbol(p.prefix, "_data_$t"), (length(p.vocab), p.batch_size)) for t = 1:p.seq_len] ∪
  [(symbol(p.prefix, "_l$(l)_init_c"), (p.dim_hidden, p.batch_size)) for l=1:p.n_layer] ∪
  [(symbol(p.prefix, "_l$(l)_init_h"), (p.dim_hidden, p.batch_size)) for l=1:p.n_layer]
end
function mx.provide_label(p :: CharSeqProvider)
  [(symbol(p.prefix, "_label_$t"), (p.batch_size,)) for t = 1:p.seq_len]
end
#--/provide

#--eachbatch-part1
function mx.eachbatch(p :: CharSeqProvider)
  data_all  = [mx.zeros(shape) for (name, shape) in mx.provide_data(p)]
  label_all = [mx.zeros(shape) for (name, shape) in mx.provide_label(p)]

  data_jl = [copy(x) for x in data_all]
  label_jl= [copy(x) for x in label_all]

  batch = mx.DataBatch(data_all, label_all, p.batch_size)
  #...
  #--/eachbatch-part1

  #--eachbatch-part2
  #...
  function _text_iter()
    text = p.text

    n_batch = floor(Int, length(text) / p.batch_size / p.seq_len)
    text = text[1:n_batch*p.batch_size*p.seq_len] # discard tailing
    idx_all = 1:length(text)

    for idx_batch in partition(idx_all, p.batch_size*p.seq_len)
      for i = 1:p.seq_len
        data_jl[i][:] = 0
        label_jl[i][:] = 0
      end

      for (i, idx_seq) in enumerate(partition(idx_batch, p.seq_len))
        for (j, idx) in enumerate(idx_seq)
          c_this = text[idx]
          c_next = idx == length(text) ? UNKNOWN_CHAR : text[idx+1]
          data_jl[j][char_idx(vocab,c_this),i] = 1
          label_jl[j][i] = char_idx(vocab,c_next)-1
        end
      end

      for i = 1:p.seq_len
        copy!(data_all[i], data_jl[i])
        copy!(label_all[i], label_jl[i])
      end

      produce(batch)
    end
  end

  return Task(_text_iter)
end
#--/eachbatch-part2

# helper function to convert a char into index in vocabulary
function char_idx(vocab :: Dict{Char,Int}, c :: Char)
  if haskey(vocab, c)
    vocab[c]
  else
    vocab[UNKNOWN_CHAR]
  end
end


text:
using IterativeSolvers
using SparseArrays
using Random
using Printf

#########################
# Method Implementation #
#########################

@inline function omega(t, s)
    angle = sqrt(2.)/2
    ns = norm(s)
    nt = norm(t)
    ts = dot(t,s)
    rho = abs(ts/(nt*ns))
    om = ts/(nt*nt)
    if rho < angle
        om = om*convert(typeof(om),angle)/rho
    end
    om
end

#s = 8, , ,
function sidrs!(X0, A, C::T; s::Number = 4, tol=missing, droptol=(eps(real(eltype(C)))), maxiter=size(A, 2), verbose::Bool=false) where {T}
	if ismissing(tol)
		tol = 10sqrt(eps(real(eltype(C))))
	end
	verbose && @show tol
	X = X0

    verbose && @printf("=== idrs ===\n%4s\t%7s\n","iter","resnorm")
    R = C - A(X)
    normR = norm(R)
	iter = 1

    if normR <= tol # Initial guess is a good enough solution
        return X
    end

    Z = zero(C)
	d = size(C,1)
    P = T[Diagonal(rand(d) .- 0.5) + sprandn(d,d, 1/d) for k in 1:s]
    U = T[copy(C) for k in 1:s]
    G = T[copy(C) for k in 1:s] #
    Q = copy(Z)
    V = copy(Z)

    M = Matrix{eltype(C)}(I,s,s)
    f = zeros(eltype(C),s)
    c = zeros(eltype(C),s)

    om::eltype(C) = 1
    while normR > tol && iter ≤ maxiter
        for i in 1:s
            f[i] = dot(P[i], R)
        end
        for k in 1:s

            # Solve small system and make v orthogonal to P

            c = LowerTriangular(M[k:s,k:s])\f[k:s]
            V = c[1] * G[k]
            Q = c[1] * U[k]

            for i = k+1:s
                V += c[i-k+1] * G[i]
                Q += c[i-k+1] * U[i]
            end

            # Compute new U[:,k] and G[:,k], G[:,k] is in space G_j
            V = R - V

            U[k] = Q + om * V
            G[k] = A(U[k])

            # Bi-orthogonalise the new basis vectors

            for i in 1:k-1
                alpha = dot(P[i],G[k])/M[i,i]
                G[k] -= alpha * G[i]
                U[k] -= alpha * U[i]
            end

            # New column of M = P'*G  (first k-1 entries are zero)

            for i in k:s
                M[i,k] = dot(P[i],G[k])
            end

            #  Make r orthogonal to q_i, i = 1..k

            beta = f[k]/M[k,k]
            R -= beta * G[k]
            X += beta * U[k]

			droptol!(X, droptol)

            normR = norm(R)
            verbose && @printf("%3d\t%1.2e\t%1.2e\n",iter,normR,sparsity(X))
            if normR < tol || iter == maxiter
                iter == maxiter && warn("not converged")
                break
            end
            if k < s
                f[k+1:s] .-=  beta*M[k+1:s,k]
            end
            iter += 1
        end # k in 1:s

		#@show sparsity(R), sparsity(X), sparsity(Q)

        # Now we have sufficient vectors in G_j to compute residual in G_j+1
        # Note: r is already perpendicular to P so v = r
        copyto!(V, R)
        Q = A(V)
        om = omega(Q, R)
        R -= om .* Q
        X += om .* V
		droptol!(X, droptol)

        normR = norm(R)

        iter += 1
    end
    verbose && print("\n")
    X
end

slyap_op(B) = X -> B*X + X*B'
linl_op(B) = X -> B*X
linr_op(B) = X -> X*B
slyapunov(B, C; tol=missing, droptol = 1e-7) = sidrs!(copy(C), slyap_op(B), -C; maxiter=prod(size(B)), verbose=true, tol=tol, droptol = droptol)
slsolve(B, Y; tol=missing, droptol = 1e-7) = sidrs!(copy(Y), linl_op(B), Y; maxiter=prod(size(B)), verbose=true, tol=tol, droptol = droptol)
srsolve(B, Y; tol=missing, droptol = 1e-7) = sidrs!(copy(Y), linr_op(B), Y; maxiter=prod(size(B)), verbose=true, tol=tol, droptol = droptol)

"""
Compute lyaponov equation using the Kronecker product for testing.
"""
kronlyap(B,C) = -reshape(inv(kron(B, Diagonal(I,size(B,1))) + kron(Diagonal(I,size(B,1)), B))*vec(C), size(C))

ext: jl


text:
include(joinpath(Pkg.dir("TrackRoots"), "src", "utils.jl"))
include(joinpath(Pkg.dir("TrackRoots"), "src", "stages.jl"))
include(joinpath(Pkg.dir("TrackRoots"), "src", "startPoints.jl"))
include(joinpath(Pkg.dir("TrackRoots"), "src", "calibrates.jl"))
include(joinpath(Pkg.dir("TrackRoots"), "src", "tracks.jl"))
include(joinpath(Pkg.dir("TrackRoots"), "src", "saves.jl"))
include(joinpath(Pkg.dir("TrackRoots"), "src", "nodes.jl"))
#=function main(stages::Vector{Stage}, startpoints::Vector{Vector{Mark}}, Δx::Float64, output::IO)
    calibstages = stages2calib(stages, Δx)
    output == STDERR && info("Calibrated the data. Tracking the roots…")
    tracks = trackroot.(calibstages, startpoints)
    output == STDERR && info("Tracked the roots. Saving the results…")
    saveit(calibstages, tracks, output)
    output == STDERR && info("Done!")
end=#
rm("/home/yakir/.julia/datadeps/all/9/pos1,3-dr5_pos2,4-cle44-luc-1-5/stage 2/", force=true, recursive=true)
mkdir("/home/yakir/.julia/datadeps/all/9/pos1,3-dr5_pos2,4-cle44-luc-1-5/stage 2")
ndfile = "/home/yakir/.julia/datadeps/all/9/pos1,3-dr5_pos2,4-cle44-luc-1-5.nd"
stages = nd2stages(ndfile)
Δx = 0.04993757802746567
startpoints = [StaticArrays.SArray{Tuple{2},Float64,1,2}[]
,StaticArrays.SArray{Tuple{2},Float64,1,2}[[511.26, 550.57]]
,StaticArrays.SArray{Tuple{2},Float64,1,2}[]
,StaticArrays.SArray{Tuple{2},Float64,1,2}[]]
# main(stages, startpoints, Float64(Δx), STDERR)
calibstages = stages2calib(stages, Δx)
tracks = trackroot.(calibstages, startpoints)
saveit(calibstages, tracks, STDERR)


text:
x = [-1, 1, 3, 3, -1]
y = [2, 0, -5, 2, -5]

@test_throws MethodError scatterplot()
@test_throws MethodError scatterplot(sin, x)
@test_throws MethodError scatterplot([sin], x)
@test_throws DimensionMismatch scatterplot([1, 2], [1, 2, 3])
@test_throws DimensionMismatch scatterplot([1, 2, 3], [1, 2])
@test_throws DimensionMismatch scatterplot([1, 2, 3], 1:2)
@test_throws DimensionMismatch scatterplot(1:3, [1, 2])
@test_throws DimensionMismatch scatterplot(1:3, 1:2)

@testset "positional types" begin
    for p in (
        @inferred(scatterplot(x, y)),
        @inferred(scatterplot(float.(x), y)),
        @inferred(scatterplot(x, float.(y))),
    )
        @test p isa Plot
        test_ref("references/scatterplot/default.txt", @show_col(p))
    end

    for p in (@inferred(scatterplot(y)), @inferred(scatterplot(float.(y))))
        @test p isa Plot
        test_ref("references/scatterplot/y_only.txt", @show_col(p))
    end

    p = @inferred scatterplot(6:10)
    @test p isa Plot
    test_ref("references/scatterplot/range1.txt", @show_col(p))

    p = @inferred scatterplot(11:15, 6:10)
    @test p isa Plot
    test_ref("references/scatterplot/range2.txt", @show_col(p))

    p = @inferred scatterplot(x .* 1e3 .+ 15, y .* 1e-3 .- 15)
    test_ref("references/scatterplot/scale1.txt", @show_col(p))
    p = @inferred scatterplot(x .* 1e-3 .+ 15, y .* 1e3 .- 15)
    test_ref("references/scatterplot/scale2.txt", @show_col(p))
    miny = -1.2796649117521434e218
    maxy = -miny
    p = @inferred scatterplot([1], [miny], xlim = (1, 1), ylim = (miny, maxy))
    test_ref("references/scatterplot/scale3.txt", @show_col(p))
    p = @inferred scatterplot([1], [miny], xlim = [1, 1], ylim = [miny, maxy])
    test_ref("references/scatterplot/scale3.txt", @show_col(p))
end

@testset "keyword arguments" begin
    p = @inferred scatterplot(x, y, xlim = (-1.5, 3.5), ylim = (-5.5, 2.5))
    test_ref("references/scatterplot/limits.txt", @show_col(p))
    p = @inferred scatterplot(x, y, xlim = [-1.5, 3.5], ylim = [-5.5, 2.5])
    test_ref("references/scatterplot/limits.txt", @show_col(p))

    p = @inferred scatterplot(x, y, grid = false)
    test_ref("references/scatterplot/nogrid.txt", @show_col(p))

    p = @inferred scatterplot(x, y, color = :blue, name = "points1")
    test_ref("references/scatterplot/blue.txt", @show_col(p))

    p = @inferred scatterplot(
        x,
        y,
        name = "points1",
        title = "Scatter",
        xlabel = "x",
        ylabel = "y",
    )
    @test p isa Plot
    test_ref("references/scatterplot/parameters1.txt", @show_col(p))

    @test @inferred(scatterplot!(p, [0.5, 1, 1.5], name = "points2")) === p
    test_ref("references/scatterplot/parameters2.txt", @show_col(p))

    @test @inferred(scatterplot!(p, [-0.5, 0.5, 1.5], [0.5, 1, 1.5], name = "points3")) ===
          p
    test_ref("references/scatterplot/parameters3.txt", @show_col(p))
    test_ref("references/scatterplot/nocolor.txt", @show_nocol(p))

    p = scatterplot(x, y, title = "Scatter", canvas = DotCanvas, width = 10, height = 5)
    @test p isa Plot
    test_ref("references/scatterplot/canvassize.txt", @show_col(p))
end

@testset "markers" begin
    p = scatterplot(
        x,
        y,
        title = "Vector of markers",
        marker = [:circle, "!", '.', :star5, :vline],
        color = [:red, :green, :yellow, :blue, :cyan],
    )
    test_ref("references/scatterplot/markers_vector.txt", @show_col(p))

    scatterplot(x, y, marker = :circle)

    n = collect(1:length(UnicodePlots.MARKERS))
    m = collect(keys(UnicodePlots.MARKERS))
    for canvas in (BrailleCanvas, DotCanvas, AsciiCanvas)
        p = scatterplot(n, n, title = "Supported markers", marker = m)
        test_ref("references/scatterplot/markers_$(canvas).txt", @show_col(p))
    end
end

@testset "densityplot" begin
    seed!(RNG, 1338)
    dx, dy = randn(RNG, 1000), randn(RNG, 1000)

    p = @inferred densityplot(dx, dy)
    @test @inferred(densityplot!(p, dx .+ 2, dy .+ 2)) === p
    @test p isa Plot
    test_ref("references/scatterplot/densityplot.txt", @show_col(p))

    p = @inferred densityplot(
        dx,
        dy,
        name = "foo",
        color = :red,
        title = "Title",
        xlabel = "x",
    )
    @test @inferred(densityplot!(p, dx .+ 2, dy .+ 2, name = "bar")) === p
    @test p isa Plot
    test_ref("references/scatterplot/densityplot_parameters.txt", @show_col(p))
end


text:
#####
##### Generate a flame graph for self-play using ProfileSVG
#####

using AlphaZero
using Profile
using ProfileSVG

experiment = Examples.experiments["tictactoe"] |> Scripts.dummy_run_experiment
network = experiment.mknet(experiment.gspec, experiment.netparams)
env = AlphaZero.Env(experiment.gspec, experiment.params, network)

AlphaZero.self_play_step!(env, nothing) # To compile every function
Profile.clear()
@profile AlphaZero.self_play_step!(env, nothing)
ProfileSVG.save("self-play-profile.svg")


text:
# Load necessary packages
using OrdinaryDiffEq
using Trixi
using Plots

###############################################################################
# semidiscretization of the linear advection equation

# Create equations
advectionvelocity = (1.0, 0.1)
equations = LinearScalarAdvectionEquation2D(advectionvelocity)

# Create DGSEM solver for polynomial degree = 3
solver = DGSEM(3, flux_lax_friedrichs)

# Create a uniformely refined mesh with periodic boundaries
coordinates_min = (-1, -1) # minimum coordinates
coordinates_max = ( 1,  1) # maximum coordinates
mesh_amr = TreeMesh(coordinates_min, coordinates_max,
                    initial_refinement_level=2,
                    n_cells_max=30_000)

# Create new initial condition for cosine pulse
function cosine_pulse(x, t, equations::LinearScalarAdvectionEquation2D)
  halfwidth = 0.3
  radius = sqrt(x[1]^2 + x[2]^2)

  if radius > halfwidth
    u = 0.0
  else
    u = 0.2 * cos(pi/2 * radius / halfwidth)
  end

  return Trixi.@SVector [u]
end

# Create semidiscretization with all spatial discretization-related components
semi = SemidiscretizationHyperbolic(mesh_amr, equations,
                                    cosine_pulse,
                                    solver)

###############################################################################
# ODE solvers, callbacks etc.

# Create ODE problem with time span from 0.0 to 1.7
ode = semidiscretize(semi, (0.0, 1.7));

# Use a simple AMR controller that to refine cells between level 4 and level 6, based on the maximum value of the solution in an element
amr_controller = ControllerThreeLevel(semi,
                                      IndicatorMax(semi, variable=first),
                                      base_level=4,
                                      med_level=5, med_threshold=0.05,
                                      max_level=6, max_threshold=0.15);

# Create the AMR callback that will be called every 5 time steps
amr_callback = AMRCallback(semi, amr_controller,
                           interval=5);

# Create visualization callback for in-situ visualization
visualization_callback = VisualizationCallback(interval=50, seriescolor=:heat);

# Create callback set to combine all created callbacks
callback_set = CallbackSet(visualization_callback, amr_callback)

###############################################################################
# run the simulation

# Create and solve ODE problem, using the previously constructed callback set
sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            callback=callback_set, dt=6.25e-3);

nothing


text:
push!(LOAD_PATH, abspath("./src"))
push!(LOAD_PATH, abspath("../MPSCircuit.jl/src"))
using Yao, MSQR, MPSCircuit

# MS Test of Cluster State.
nMeasure = 20000
nBitT = 5
vBit = 1
rBit = 1

regTar = rand_state(nBitT)
MPSGen = MPSC("CS", nBitT, vBit, rBit)
circuit = MScircuit(nBitT, vBit, rBit, MPSGen.mpsBlocks)
test1 = MSTtest(regTar, circuit, MPSGen.cExtend, nMeasure=nMeasure)
ExpectOverlap = test1.Eoverlap
ActualOverlap = test1.Aoverlap
println("ExpectOverlap: $(ExpectOverlap)\nActualOverlap: $(ActualOverlap)")
@show isapprox.(ActualOverlap, ExpectOverlap, atol=0.05*ExpectOverlap)
## MS Test of Extended MPS Cluster-state circuit.
ExtendOverlap = SWAPtest(zero_state(nBitT)|>MPSGen.cExtend, regTar, nMeasure=nMeasure).overlap
println("ExtendOverlap: $(ExtendOverlap)\n")
@show isapprox.(ExtendOverlap, ExpectOverlap, atol=0.05*ExpectOverlap)

# MS Test of differentiable quantum state.
nMeasure = 20000
nBitT = 4
vBit = 2
rBit = 1
depth = 5

regTar = rand_state(nBitT)
MPSGen = MPSC(("DC",depth),nBitT,vBit,rBit)
circuit = MScircuit(nBitT, vBit, rBit, MPSGen.mpsBlocks)
test2 = MSTtest(regTar, circuit, MPSGen.cExtend, nMeasure=nMeasure)
ExpectOverlap = test2.Eoverlap
ActualOverlap = test2.Aoverlap
println("ExpectOverlap: $(ExpectOverlap)\nActualOverlap: $(ActualOverlap)")
@show isapprox.(ActualOverlap, ExpectOverlap, atol=0.05*ExpectOverlap)
## MS Test of Extended MPS differentiable circuit.
ExtendOverlap = SWAPtest(zero_state(nBitT)|>MPSGen.cExtend, regTar, nMeasure=nMeasure).overlap
println("ExtendOverlap: $(ExtendOverlap)\n")
@show isapprox.(ExtendOverlap, ExpectOverlap, atol=0.05*ExpectOverlap)


text:
const TypesWithShow = Union{AbstractVariable, AbstractMonomial, AbstractTerm,
                            AbstractPolynomial, RationalPoly}
function Base.show(io::IO, mime::MIME"text/latex", p::TypesWithShow)
    print(io, "\$\$ ")
    _show(io, mime, p)
    print(io, " \$\$")
end

# If the MIME is not specified, IJulia thinks that it supports images, ...
# and then use the result of show and tries to interpret it as an svg, ...
# We need the two methods to avoid ambiguity
function Base.show(io::IO, mime::MIME"text/plain", p::TypesWithShow)
    _show(io, mime, p)
end
function Base.show(io::IO, mime::MIME"text/print", p::TypesWithShow)
    _show(io, mime, p)
end

Base.print(io::IO, p::TypesWithShow) = show(io, MIME"text/print"(), p)
Base.show(io::IO, p::TypesWithShow) = show(io, MIME"text/plain"(), p)

# VARIABLES
function _show(io::IO, mime::MIME, var::AbstractVariable)
    base, indices = name_base_indices(var)
    if isempty(indices)
        print(io, base)
    else
        print(io, base)
        print_subscript(io, mime, indices)
    end
end
_show(io::IO, mime::MIME"text/print", var::AbstractVariable) = print(io, name(var))

function print_subscript(io::IO, ::MIME"text/latex", index)
    print(io, "_{", join(index, ","), "}")
end
function print_subscript(io::IO, mime, indices)
    if length(indices) == 1
        print(io, unicode_subscript(indices[1]))
    else
        print(io, join(unicode_subscript.(indices), "\u208B"))
    end
end

const unicode_subscripts = ("₀","₁","₂","₃","₄","₅","₆","₇","₈","₉")
unicode_subscript(i) = join(unicode_subscripts[d+1] for d in reverse(digits(i)))

# MONOMIALS
function _show(io::IO, mime, m::AbstractMonomial)
    if isconstant(m)
        print(io, '1')
    else
        printed_var = false
        vars = variables(m)
        n = length(vars)
        for (i, var, exp) in zip(1:n, vars, exponents(m))
            if !iszero(exp)
                if mime isa MIME"text/print" && printed_var && i > 0 &&
                    print(io,"*")
                end
                _show(io, mime, var)
                printed_var = true
                if !isone(exp)
                    print_exponent(io, mime, exp)
                end
            end
        end
    end
end

print_exponent(io::IO, ::MIME"text/latex", exp) = print(io, "^{", exp, "}")
print_exponent(io::IO, ::MIME"text/print", exp) = print(io, "^", exp)
function print_exponent(io::IO, mime, exp)
    print(io, join(unicode_superscript.(reverse(digits(exp)))))
end

const unicode_superscripts = ("⁰","¹","²","³","⁴","⁵","⁶","⁷","⁸","⁹")
unicode_superscript(i) = unicode_superscripts[i+1]

# TERM
function _show(io::IO, mime, t::AbstractTerm)
    if isconstant(t)
        print_coefficient(io, coefficient(t))
    else
        if should_print_coefficient(coefficient(t))
            if !should_print_coefficient(-coefficient(t))
                print(io, '-')
            else
                print_coefficient(io, coefficient(t))
                if !iszero(t)
                    # Print a multiplication sign between coefficent and monmomial
                    # depending on the mime type
                    print_maybe_multiplication_sign(io, mime)
                end
            end
        end
        if !iszero(t)
            _show(io, mime, monomial(t))
        end
    end
end

"""
    print_maybe_multiplication_sign(io, mime)

Prints a multiplication sign depending on the `mime` type.
"""
print_maybe_multiplication_sign(io::IO, ::MIME"text/print") = print(io, "*")
print_maybe_multiplication_sign(io::IO, mime) = nothing

should_print_coefficient(x) = true  # By default, just print all coefficients
should_print_coefficient(x::Number) = !isone(x) # For numbers, we omit any "one" coefficients
print_coefficient(io::IO, coeff::Real) = print(io, coeff)
print_coefficient(io::IO, coeff) = print(io, "(", coeff, ")")

# POLYNOMIAL
function _show(io::IO, mime, p::AbstractPolynomial{T}) where T
    ts = terms(p)
    if isempty(ts)
        print(io, zero(T))
    else
        _show(io, mime, first(ts))
        for t in Iterators.drop(ts, 1)
            if isnegative(coefficient(t))
                print(io, " - ")
                _show(io, mime, abs(coefficient(t)) * monomial(t))
            else
                print(io, " + ")
                _show(io, mime, t)
            end
        end
    end
end

isnegative(x::Real) = x < 0
isnegative(x) = false

# RATIONAL POLY
function _show(io::IO, mime, p::RationalPoly)
    print(io, mime isa MIME"text/latex" ? "\\frac{" : "(")
    _show(io, mime, p.num)
    print(io, mime isa MIME"text/latex" ? "}{" : ") / (")
    _show(io, mime, p.den)
    print(io, mime isa MIME"text/latex" ? "}" : ")")
end


text:
using PushGP: PushInterpreter, execute!, stack, last, length, Add, Mul, Sub, Div, AfromB

@testset "Examples used in PushGP.jl README" begin

@testset "Example 1" begin
    simple_program = [2, 5, Add{Int}()]

    i = PushInterpreter()
    ret = execute!(i, simple_program)
    @test ret === nothing

    intstack = stack(i, Int)
    @test isa(intstack, PushGP.Stack{Int})
    @test last(intstack) == 7
    @test length(intstack) == 1
end

@testset "Example 2" begin
    # We assume we are continuing after example 1 above but to ensure tests
    # can be run out of order we rerun its example
    i = PushInterpreter()
    simple_program = [2, 5, Add{Int}()]
    execute!(i, simple_program)

    # Now we continue with example 2:
    execute!(i, [6.5, 5, Sub{Int}(), 2.0, Mul{Float64}(), 4.0, Div{Float64}(), AfromB{String,Float64}(), 1.7])
    @test last(stack(i, Int)) == 2 # 7 (from before) - 5
    @test length(stack(i, Int)) == 1
    @test last(stack(i, String)) == "3.25" # string((2.0 * 6.5) / 4.0)
    @test length(stack(i, String)) == 1
    @test last(stack(i, Float64)) == 1.7
    @test length(stack(i, Float64)) == 1
end
end # @testset "Examples from PushGP docs" begin


text:
using NURBS
using Base.Test

@testset "bsurface Building" begin

	arr = [-100.0 -100.0 -100.0 -100.0 -100.0 -50.0 -50.0 -50.0 -50.0 -50.0 0.0 0.0 0.0 0.0 0.0 50.0 50.0 50.0 50.0 50.0 100.0 100.0 100.0 100.0 100.0;-100.0 -50.0 0.0 50.0 100.0 -100.0 -50.0 0.0 50.0 100.0 -100.0 -50.0 0.0 50.0 100.0 -100.0 -50.0 0.0 50.0 100.0 -100.0 -50.0 0.0 50.0 100.0;0.0 0.0 0.0 0.0 0.0 0.0 25.0 50.0 25.0 0.0 0.0 25.0 50.0 25.0 0.0 0.0 25.0 150.0 25.0 0.0 0.0 0.0 0.0 0.0 0.0]

	res = [-100.0 -100.0 -100.0 -100.0 -25.9259 -25.9259 -25.9259 -25.9259 25.9259 25.9259 25.9259 25.9259 100.0 100.0 100.0 100.0; -100.0 -25.9259 25.9259 100.0 -100.0 -25.9259 25.9259 100.0 -100.0 -25.9259 25.9259 100.0 -100.0 -25.9259 25.9259 100.0; 0.0 0.0 0.0 0.0 0.0 34.8422 34.8422 0.0 0.0 51.3032 51.3032 0.0 0.0 0.0 0.0 0.0]

    resu = [
        -50.0 -50.0 -50.0 -50.0 -16.6667 -16.6667 -16.6667  -16.6667 16.6667 16.6667 16.6667 16.6667 50.0 50.0 50.0 50.0;
        -50.0 -16.6667 16.6667 50.0 -50.0 -16.6667 16.6667 50.0 -50.0 -16.6667 16.6667 50.0 -50.0 -16.6667 16.6667 50.0;
        20.8333 32.6646 32.6646 20.8333 25.6687 41.7905 41.7905 25.6687 31.0185 60.2176 60.2176 31.0185 31.9444 70.9362 70.9362 31.9444
    ]

    @testset "B-spline surface using open uniform knot vectors" begin
        @test typeof(bsplsurf(arr, 4, 4, 5, 5, 100, 100)) == Tuple{Array{Float64, 2}, Array{Array{Int64, 1}, 1}, Array{Array{Int64, 1}, 1}}
        @test isapprox(bsplsurf(arr, 4, 4, 5, 5, 4, 4)[1], res; atol = 1e-3)
    end

    @testset "B-spline surface using periodic uniform knot vectors" begin
        @test typeof(bsplsurfu(arr, 4, 4, 5, 5, 100, 100)) == Tuple{Array{Float64, 2}, Array{Array{Int64, 1}, 1}, Array{Array{Int64, 1}, 1}}
        @test isapprox(bsplsurfu(arr, 4, 4, 5, 5, 4, 4)[1], resu; atol = 1e-3)
    end

    @testset "B-spline surface and its derivatives using open uniform knot vectors" begin
        @test typeof(dbsurf(arr, 4, 4, 5, 5, 100, 100)) == Tuple{Array{Float64,2}, Array{Float64,2}, Array{Float64,2}, Array{Float64,2}, Array{Float64,2}, Array{Float64,2}, Array{Array{Int64, 1}, 1}, Array{Array{Int64, 1}, 1}}
        @test isapprox(dbsurf(arr, 4, 4, 5, 5, 4, 4)[1], res; atol = 1e-3)
    end

    @testset "B-Spline surface over open knot comparison" begin
        @test bsplsurf(arr, 4, 4, 5, 5, 100, 100)[1] == dbsurf(arr, 4, 4, 5, 5, 100, 100)[1]
        @test bsplsurf(arr, 4, 4, 5, 5, 100, 100)[2] == dbsurf(arr, 4, 4, 5, 5, 100, 100)[7]
        @test bsplsurf(arr, 4, 4, 5, 5, 100, 100)[3] == dbsurf(arr, 4, 4, 5, 5, 100, 100)[8]
    end
end


text:
using WebPlotDigitizer
using Test
using Dates

testdir(args...) = abspath(joinpath(dirname(pathof(WebPlotDigitizer)),"..","test",args...))
datadir(args...) = testdir("data",args...)

function testDate(dt::Date)
    res = WebPlotDigitizer.parse_import_value(Dates.datetime2unix(DateTime(dt))*1000.0,"2020/10/10")::Date
    res === res
end
function testDate(dt::DateTime)
    res = WebPlotDigitizer.parse_import_value(Dates.datetime2unix(DateTime(dt))*1000.0,"2020/10/10 10:00")::DateTime
    res === res
end
function testDate(dt::Time)
    res = WebPlotDigitizer.parse_import_value(Dates.datetime2unix(today()+dt)*1000.0,"10:00")::Time
    res === res
end

@testset "WebPlotDigitizer" begin
    @testset "Loading" begin
        wpd = WebPlotDigitizer.load_project(datadir("default.tar"))
        @test collect(keys(wpd.axes)) == ["XY"]
        @test collect(keys(wpd.axes["XY"].data)) == ["Dataset 1", "Dataset 2", "Dataset 3", "Dataset 4"]
        @test_throws ErrorException WebPlotDigitizer.load_from_tar(datadir("no_wpd.tar"))
        newname = WebPlotDigitizer.load_project(datadir("newname.tar"))
        @test collect(keys(newname.axes)) == ["XY", "XY 2"]
        @test collect(keys(newname.axes["XY"].data)) == ["Dataset 1", "Dataset 2", "Dataset 3"]
        @test collect(keys(newname.axes["XY 2"].data)) == ["Dataset 4"]
    end
    @testset "default.tar" begin
        wpd = WebPlotDigitizer.load_project(datadir("default.tar"))
        @test wpd["XY"] isa WebPlotDigitizer.XYAxes
        @test wpd["XY"].isLogX == false
        @test wpd["XY"].isLogY == false
        wpd
    end
    @testset "basictypes.tar" begin
        wpd = WebPlotDigitizer.load_project(datadir("basictypes.tar"))
        @test all(wpd["XY Date"]["Dataset 4"][:,1] .== [Date(2020,1,1), Date(2020,1,3), Date(2020,2,1), Date(2020,2,13)])
        @test all(isapprox.(wpd["XY Date"]["Dataset 4"][:,2], [7,2,3,1], atol=0.1))
        @test all(isapprox.(wpd["XY"]["Dataset 0"][:,1], [0,1,3,5], atol=0.1))
        @test all(isapprox.(wpd["XY"]["Dataset 0"][:,2], [7,2,3,1], atol=0.1))
        @test all(isapprox.(wpd["XY Log"]["Dataset 3"][:,1], [1, 10, 100, 1000], atol=0.01))
        @test all(isapprox.(wpd["XY Log"]["Dataset 3"][:,2], [7,2,3,1], atol=0.1))
        @test all(isapprox.(wpd["Bar"]["Dataset 1"][:,1], [10,2,5,11], atol=0.1))
        @test all(isapprox.(wpd["Polar"]["Dataset 2"][:,2], rad2deg.([0,0.3,0.6,0.7]*2π), atol=0.5))
        @test all(isapprox.(wpd["Polar"]["Dataset 2"][:,1], [1,2,3,2], atol=0.1))
    end
    @testset "dateandtime.tar" begin
        wpd = WebPlotDigitizer.load_project(datadir("dateandtime.tar"))
        x1 = [Time(10,1,3), Time(10,2,3), Time(11,3,4), Time(13,10,3)]
        y1 = [1, 10, -3, 3]
        x2 = Date(2020,9,1) .+ x1
        y2 = [-1, 6, 3, -10]

        @test all((wpd["Time"]["D1"][:,1] .- x1) .< Dates.Second(10))
        @test all(isapprox.(wpd["Time"]["D1"][:,2], y1, atol=0.1))
        @test all((wpd["Datetime"]["D2"][:,1] .- x2) .< Dates.Second(10))
        @test all(isapprox.(wpd["Datetime"]["D2"][:,2], y2, atol=0.1))
    end
    @testset "Sorting" begin
        wpd = WebPlotDigitizer.load_project(datadir("default.tar"))
        @test issorted(wpd["XY"]["Dataset 2"][:,2]) == false
        sortby!(wpd["XY"],2)
        @test issorted(wpd["XY"]["Dataset 2"][:,2])
    end
    @testset "Date parser" begin
        @test testDate(Date(2020,10,1))
        @test testDate(DateTime(2020,10,1,10,1,2))
        @test testDate(Time(10,1))
    end
end


text:
# mtDNA model by de Oliveira et al.
include("common.jl")

@with_kw struct MTDNAParams
    MT_DNA = 0.75  # Default Mitocondrial DNA content
    α = 5.625e-11  # Rate constant of mtDNA reapir
    β = 6.940e-11  # reaction rate constant of hydroxyl radical and mtDNA
    γ = 8.100e-9  # reaction rate constant of DOX and mtDNA
    κ = 0.045  # half saturation concentration of mtDNA repair
    K_B = 1.42   # rate constant of mt protein synthesis
    K_PROT = 0.32  # half saturation concentration of mt protein synthesis
    MT_PROT = K_B * _mm(MT_DNA, K_PROT)
end


text:
using CLIMA.PlanetParameters: Omega
export Gravity, RayleighSponge, Subsidence, GeostrophicForcing, Coriolis

# kept for compatibility
# can be removed if no functions are using this
function atmos_source!(f::Function, atmos::AtmosModel, source::Vars, state::Vars, aux::Vars, t::Real)
  f(source, state, aux, t)
end
function atmos_source!(::Nothing, atmos::AtmosModel, source::Vars, state::Vars, aux::Vars, t::Real)
end
# sources are applied additively
function atmos_source!(stuple::Tuple, atmos::AtmosModel, source::Vars, state::Vars, aux::Vars, t::Real)
  map(s -> atmos_source!(s, atmos, source, state, aux, t), stuple)
end

abstract type Source
end

struct Gravity <: Source
end
function atmos_source!(::Gravity, atmos::AtmosModel, source::Vars, state::Vars, aux::Vars, t::Real)
  if atmos.ref_state isa HydrostaticState
    source.ρu -= (state.ρ - aux.ref_state.ρ) * aux.orientation.∇Φ
  else
    source.ρu -= state.ρ * aux.orientation.∇Φ
  end
end

struct Coriolis <: Source
end
function atmos_source!(::Coriolis, atmos::AtmosModel, source::Vars, state::Vars, aux::Vars, t::Real)
  # note: this assumes a SphericalOrientation
  source.ρu -= SVector(0, 0, 2*Omega) × state.ρu
end

struct GeostrophicForcing{FT} <: Source
  f_coriolis::FT
  u_geostrophic::FT
  v_geostrophic::FT
end
function atmos_source!(s::GeostrophicForcing, atmos::AtmosModel, source::Vars, state::Vars, aux::Vars, t::Real)
  u_geo = SVector(s.u_geostrophic, s.v_geostrophic, 0)
  ẑ = vertical_unit_vector(atmos.orientation, aux)
  fkvector = s.f_coriolis * ẑ
  source.ρu -= fkvector × (state.ρu .- state.ρ*u_geo)
end

"""
    RayleighSponge{FT} <: Source

Rayleigh Damping (Linear Relaxation) for top wall momentum components
Assumes laterally periodic boundary conditions for LES flows. Momentum components
are relaxed to reference values (zero velocities) at the top boundary.
"""
struct RayleighSponge{FT} <: Source
  "Maximum domain altitude (m)"
  z_max::FT
  "Altitude at with sponge starts (m)"
  z_sponge::FT
  "Sponge Strength 0 ⩽ α_max ⩽ 1"
  α_max::FT
  "Relaxation velocity components"
  u_relaxation::SVector{3,FT}
  "Sponge exponent"
  γ::FT
end
function atmos_source!(s::RayleighSponge, atmos::AtmosModel, source::Vars, state::Vars, aux::Vars, t::Real)
  z = altitude(atmos.orientation, aux)
  if z >= s.z_sponge
    r = (z - s.z_sponge)/(s.z_max-s.z_sponge)
    β_sponge = s.α_max * sinpi(r/2)^s.γ
    source.ρu -= β_sponge * (state.ρu .- state.ρ*s.u_relaxation)
  end
end

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
const CenteredOrUpwind = Union{AbstractCenteredAdvectionScheme, AbstractUpwindBiasedAdvectionScheme}

#####
##### Flat Topologies
#####

momentum_flux_uu(i, j, k, grid::AbstractGrid{FT, Flat}, scheme::CenteredOrUpwind, u)    where FT = zero(FT)
momentum_flux_vu(i, j, k, grid::AbstractGrid{FT, Flat}, scheme::CenteredOrUpwind, u, v) where FT = zero(FT)
momentum_flux_wu(i, j, k, grid::AbstractGrid{FT, Flat}, scheme::CenteredOrUpwind, u, w) where FT = zero(FT)

momentum_flux_uv(i, j, k, grid::AbstractGrid{FT, TX, Flat}, scheme::CenteredOrUpwind, u, v) where {FT, TX} = zero(FT)
momentum_flux_vv(i, j, k, grid::AbstractGrid{FT, TX, Flat}, scheme::CenteredOrUpwind, v)    where {FT, TX} = zero(FT)
momentum_flux_wv(i, j, k, grid::AbstractGrid{FT, TX, Flat}, scheme::CenteredOrUpwind, v, w) where {FT, TX} = zero(FT)

momentum_flux_uw(i, j, k, grid::AbstractGrid{FT, TX, TY, Flat}, scheme::CenteredOrUpwind, u, v) where {FT, TX, TY} = zero(FT)
momentum_flux_vw(i, j, k, grid::AbstractGrid{FT, TX, TY, Flat}, scheme::CenteredOrUpwind, v, w) where {FT, TX, TY} = zero(FT)
momentum_flux_ww(i, j, k, grid::AbstractGrid{FT, TX, TY, Flat}, scheme::CenteredOrUpwind, w)    where {FT, TX, TY} = zero(FT)


text:
using SimpleGraphs: vlist, elist, out_neighbors

# Type piracy, this should really go in SimpleGraphs
# OTOH there's not anything else this could really mean
# function SimpleGraphs.components(g::SimpleDigraph)
#     SimpleGraphs.components(convert(SimpleGraph, g))
# end

function _before(g::SimpleDigraph, v)
    parents = in_neighbors(g, v)
    for i in parents
        append!(parents, _before(g, i))
    end
    return parents
end

before(g::SimpleDigraph, v; inclusive = true) = inclusive ? push!(_before(g, v), v) : _before(g, v)
function before(g::SimpleDigraph, vs...; inclusive = true)
    parents = inclusive ? collect(vs) : Symbol[]
    for v in vs
        append!(parents, before(g, v, inclusive = inclusive))
    end
    if !inclusive
        setdiff!(parents, vs)
    end
    return unique!(parents)
end

notbefore(g::SimpleDigraph, vs...; inclusive = false) = setdiff(vlist(g), before(g, vs...; inclusive = !inclusive))

function _after(g::SimpleDigraph, v)
    children = out_neighbors(g, v)
    for i in children
        append!(children, _after(g, i))
    end
    return children
end

after(g::SimpleDigraph, v; inclusive = true) = inclusive ? push!(_after(g, v), v) : _after(g, v)

function after(g::SimpleDigraph, vs...; inclusive = true)
    children = inclusive ? collect(vs) : Symbol[]
    for v in vs
        append!(children, after(g, v, inclusive = inclusive))
    end
    if !inclusive
        setdiff!(children, vs)
    end
    return unique!(children)
end

notafter(g::SimpleDigraph, vs...; inclusive = false) = setdiff(vlist(g), after(g, vs...; inclusive = !inclusive))

function assemblefrom(m::Model, params, args)
    theModule = getmodule(m)
    m_init = Model(theModule, args, NamedTuple(), NamedTuple(), nothing)
    m = foldl(params; init=m_init) do m0,v
        merge(m0, Model(theModule, findStatement(m, v)))
    end
    return m
end

getReturn(am::AbstractModel) = Model(am).retn

function setReturn(m::Model, x)
    theModule = getmodule(m)
    m0 = assemblefrom(m, parameters(m), arguments(m))
    isnothing(x) && return m0
    return merge(m0, Model(theModule, Return(x)))
end

function trim_args!(args, m, params)
    g = digraph(m)
    intersect!(args, setdiff(vcat((parents(g, p) for p in params)...), params))
end

sourcenodes(g::SimpleDigraph) = setdiff(vlist(g), [last(e) for e in elist(g)])
sinknodes(g::SimpleDigraph) = setdiff(vlist(g), [first(e) for e in elist(g)])

"""
    after(m::Model, xs...; strict=false)

Transforms `m` by moving `xs` to arguments. If `strict=true`, only descendants of `xs` are retained in the body. Otherwise, the remaining variables in the body are unmodified. Unused arguments are trimmed.

`predictive(m::Model, xs...) = after(m, xs..., strict = true)`

`Do(m::Model, xs...) = after(m, xs..., strict = false)`

# Example
```jldoctest
julia> m = @model (n, k) begin
           β ~ Gamma()
           α ~ Gamma()
           θ ~ Beta(α, β)
           x ~ Binomial(n, θ)
           z ~ Binomial(k, α / (α + β))
       end;

julia> Soss.after(m, :θ, strict=false) # same as Do(m, :θ)
@model (n, k, θ) begin
        β ~ Gamma()
        α ~ Gamma()
        x ~ Binomial(n, θ)
        z ~ Binomial(k, α / (α + β))
    end

julia> Soss.after(m, :θ, strict = true) # same as predictive(m, :θ)
@model (n, θ) begin
        x ~ Binomial(n, θ)
    end
```
"""
function after(m::Model, xs...; strict = false)
    g = digraph(m)
    for x in xs
        map(a->delete!(g, a, x), in_neighbors(g, x))
    end
    parms = Symbol[]
    for v in sinknodes(g)
        sinkparents = before(g, v)
        !strict || any(in(sinkparents).(xs)) ? append!(parms, sinkparents) : nothing
    end
    args = arguments(m) ∪ xs # Will trim later.
    parms = setdiff(parms, args)
    trim_args!(args, m, parms)
    m1 = setReturn(assemblefrom(m, parms, args), getReturn(m))
    newargs = arguments(m1) ∪ setdiff(variables(m) ∩ variables(getReturn(m1)), variables(m1))
    return merge(m1, Model(getmodule(m), newargs, quote end))
end

"""
    before(m::Model, xs...; inclusive=true, strict=true)

Transforms `m` by retaining all ancestors of any of `xs` if `strict=true`; if `strict=false`, retains all variables that are not descendants of any `xs`. Note that adding more variables to `xs` cannot result in a larger model. If `inclusive=true`, `xs` is considered to be an ancestor of itself and is always included in the returned `Model`. Unneeded arguments are trimmed.

`prune(m::Model, xs...) = before(m, xs..., inclusive = false, strict = false)`

`prior(m::Model, xs...) = before(m, xs..., inclusive = true, strict = true)`

# Examples
```jldoctest
julia> m = @model (n, k) begin
           β ~ Gamma()
           α ~ Gamma()
           θ ~ Beta(α, β)
           x ~ Binomial(n, θ)
           z ~ Binomial(k, α / (α + β))
       end;

julia> Soss.before(m, :θ, inclusive = true, strict = true)
@model begin
        β ~ Gamma()
        α ~ Gamma()
        θ ~ Beta(α, β)
    end

julia> Soss.before(m, :θ, inclusive = true, strict = false)
@model k begin
        β ~ Gamma()
        α ~ Gamma()
        θ ~ Beta(α, β)
        z ~ Binomial(k, α / (α + β))
    end

julia> Soss.before(m, :θ, inclusive = false, strict = true) # same as Soss.prior(m, :θ)
@model begin
        β ~ Gamma()
        α ~ Gamma()
    end

julia> Soss.before(m, :θ, inclusive=false, strict=false) # same as Soss.prune(m, :θ)
@model k begin
        β ~ Gamma()
        α ~ Gamma()
        z ~ Binomial(k, α / (α + β))
    end
```
"""
function before(m::Model, xs...; inclusive = true, strict = true)
    g = digraph(m)
    vars = strict ? before(g, xs..., inclusive = inclusive) : notafter(g, xs..., inclusive = inclusive)
    args = arguments(m) ∩ vars
    params = setdiff(vars, args)
    trim_args!(args, m, params)
    return assemblefrom(m, params, args)
end


text:
using Markdown
using InteractiveUtils

using ExtremeStats, PyCall, Conda, DataFrames

Conda.pip_interop(true, "/Users/david/.julia/conda/3")

Conda.pip("install","numerapi")

Conda.pip("install","pandas")

begin
	numerapi = pyimport("numerapi")
	pandas = pyimport("pandas")
end

napi = numerapi.NumerAPI()

begin
	raw_data = pandas.json_normalize(napi.daily_user_performances("david_plutus"), sep="-")
	raw_data.to_csv("data.csv")
end


text:
#using Distributed
using ClearStacktrace

using Plots

using SpikeNetOpt
SNO = SpikeNetOpt
#@show(varinfo(SNO))

using SpikingNeuralNetworks
using Evolutionary
SNN = SpikingNeuralNetworks
SNN.@load_units
unicodeplots()

##
##

##
# Ground truths
##
global Ne = 200;
global Ni = 50
global σee = 1.0
global pee = 0.5
global σei = 1.0
global pei = 0.5
MU = 10

global E
global spkd_ground
global GT = 26

g, Cg = SpikeNetOpt.make_net_from_graph_structure(GT)

P, C = SpikeNetOpt.make_net_SNN(Ne, Ni, σee = 0.5, pee = 0.8, σei = 0.5, pei = 0.8)
E, I = P
EE, EI, IE, II = C
SNN.monitor([E, I], [:fire])
sim_length = 1000
@inbounds for t = 1:sim_length*ms
    E.I = vec([11.5 for i = 1:sim_length])#vec(E_stim[t,:])#[i]#3randn(Ne)
    SNN.sim!(P, C, 1ms)

end
#spkd_ground = SpikeNetOpt.get(P[1])
spkd_ground = SpikeNetOpt.get_trains(P[1])
#@show(spkd_ground)
display(SNN.raster(P[1]))
SNN.raster(P[1]) |> display

#sgg = [convert(Array{Float32,1}, sg) for sg in spkd_ground]
P
#P, C = SpikeNetOpt.make_net(GT)
##
# Ground truth for optimization
##
#spkd_ground = get_trains(Pg[1])
##
# Not really necessary
##
#sgg = [convert(Array{Float32,1}, sg) for sg in spkd_ground]

META_HEUR_OPT = true

function loss(model)
    @show(Ne,Ni)
    @show(model)

    σee = model[1]
    pee = model[2]
    σei = model[3]
    pei = model[4]
    P1, C1 = SNO.make_net_SNN(Ne, Ni, σee = σee, pee = pee, σei = σei, pei = pei)#,a=a)
    E1, I1 = P1
    SNN.monitor([E1, I1], [:fire])
    sim_length = 1000
    @inbounds for t = 1:sim_length*ms
        E1.I = vec([11.5 for i = 1:sim_length])
        SNN.sim!(P1, C1, 1ms)
    end

    spkd_found = get_trains(P1[1])
    println("Ground Truth \n")
    SNN.raster([E]) |> display
    println("Best Candidate \n")
    SNN.raster([E1]) |> display
    error = SNO.spike_train_difference(spkd_ground, spkd_found)
    error = sum(error)
    @show(error)
    error
end

if META_HEUR_OPT
    D = 10
    bounds = [3ones(D) 40ones(D)]'
    a = view(bounds, 1, 1)
    b = view(bounds, 1, 2)
    information = Information(f_optimum = 0.0)
    options = Options( seed = 1, iterations=10, f_calls_limit =10)

    D = size(bounds, 2)
    nobjectives=1

    options = Options( seed = 1, iterations=10000, f_calls_limit = 25000)
    npartitions = nobjectives == 2 ? 100 : 12

    methods = [
            SMS_EMOA(N = 5, n_samples=5, options=options),
            NSGA2(options=options)
            ]

    for method in methods
        f_calls = 0
        result = optimize(SpikeNetOpt.loss, bounds, method)
        #result = optimize(SpikeNetOpt.loss, bounds)#, method)
        @show(result)

    end

end

# ## Visualize the result

#contour(0:.1:1, 0:.1:1, (x, y) -> model([x,y])[].data, fill=true)


text:
"""
    Const("anything")

DataType which behaves constant among `map`, `foreach` and the like. Just like an empty container,
however with additional information about which kind of empty.
"""
struct Const{T}
  value::T
  Const(value) = new{typeof(value)}(value)
end

# == controversy https://github.com/JuliaLang/julia/issues/4648
Base.:(==)(a::Const, b::Const) = a.value == b.value

function Base.show(io::IO, x::Const)
  print(io, "Const($(repr(x.value)))")
end

"""
    isconst(Const(3)) -> true
    isconst("anythingelse") -> false

returns true only if given an instance of  [`DataTypesBasic.Const`](@ref)
"""
Base.isconst(::Const) = true
Base.isconst(other) = false

# length 0 corresponds directly to Base.iterate
Base.length(::Const) = 0

# for convenience we support different access functions so that users do not have to remember the name of the single field "value"
Base.get(c::Const) = c.value
Base.getindex(c::Const) = c.value

# const just does nothing, i.e. leaves everything constant
Base.iterate(c::Const) = nothing
Base.foreach(f, c::Const) = nothing
Base.map(f, c::Const) = c
Iterators.flatten(c::Const) = c

Base.eltype(::Type{<:Const}) = Any

Base.convert(::Type{Const{T}}, x::Const{T}) where {T} = x
Base.convert(::Type{Const{T}}, x::Const) where {T} = Const(convert(T, x.value))

# Const is covariate
# promote_rule only works on concrete types, more general checks Type{<:Const} may overwrite
# unintentionally more specific promote_rule types
Base.promote_rule(::Type{Const{T}}, ::Type{Const{S}}) where {T, S<:T} = Const{T}
Base.promote_rule(::Type{Const{T}}, ::Type{Const}) where T = Const
Base.promote_rule(::Type{Const}, ::Type{Const}) = Const

# we need this for safety, if someone overwrites promote_typejoin for Unions with Const
Base.promote_typejoin(::Type{Const{T}}, ::Type{Const{T}}) where T = Const{T}
Base.promote_typejoin(::Type{<:Const}, ::Type{<:Const}) = Const


text:
# ************************************************************************
# FAUST Architecture File
# Copyright (C) 2021 GRAME, Centre National de Creation Musicale
# ---------------------------------------------------------------------

# This is sample code. This file is provided as an example of minimal
# FAUST architecture file. Redistribution and use in source and binary
# forms, with or without modification, in part or in full are permitted.
# In particular you can create a derived work of this FAUST architecture
# and distribute that work under terms of your choice.

# This sample code 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.
# ************************************************************************

# Architectures files
#include("../dsp/dsp.jl")

abstract type audio end

function init!(driver::audio, name::String, dsp::dsp)
end

# Blocking version
function run!(driver::audio)
end

function start!(driver::audio)
end

function stop!(driver::audio)
end


text:
using Winston

# optionally call figure prior to plotting to set the size
figure(width=600, height=400)
# plot some data
pl = plot(cumsum(rand(500) .- 0.5), "r", cumsum(rand(500) .- 0.5), "b")
# display the plot (not done automatically!)
display(pl)

# by default display will not wait and the plot will vanish as soon as it appears
# using readline is a blunt wait to allow the user to choose when to continue
println("Press enter to continue: ")
readline(STDIN)

# save the current figure
savefig("winston.svg")
# .eps, .pdf, & .png are also supported
# we used svg here because it respects the width and height specified above

# TODO: needs more expansive explanation

max_stars_repo_name: floswald/JuliaByExample

max_issues_repo_name: floswald/JuliaByExample

max_forks_repo_name: floswald/JuliaByExample


text:
# Haskell

function haskellsetup()
    quote
        const haskell_filename = Ref{AbstractString}("code.hs")
        function set_haskell_filename!(filename::AbstractString)
            return haskell_filename[] = filename
        end

        const haskell_chunks = Dict{AbstractString,Vector{String}}()
        function haskell_add_code!(haskell_code::AbstractString)
            return push!(get!(haskell_chunks, haskell_filename[], String[]), haskell_code)
        end
        function haskell_get_code()
            iobuffer = IOBuffer()
            allcode = Dict{AbstractString,String}()
            for (filename, chunks) in haskell_chunks
                for chunk in chunks
                    println(iobuffer, chunk)
                end
                allcode[filename] = String(take!(iobuffer))
            end
            return allcode
        end
        function haskell_write_code!()
            println("Generating Haskell code:")
            allcode = haskell_get_code()
            for (filename, content) in allcode
                println("Generating $filename...")
                write(filename, content)
            end
            println("Done.")
            return nothing
        end
    end
end
export haskellsetup

function haskellnewfile(filename::AbstractString, haskell_stmts::AbstractString)
    haskell_code = """
               {-# LANGUAGE ForeignFunctionInterface #-}
               module $(split(filename, ".")[1]) where
               import Foreign.C.Types

               $haskell_stmts
               """
    quote
        set_haskell_filename!($filename)
        haskell_add_code!($haskell_code)
    end
end
export haskellnewfile

function haskellcode(haskell_stmts::AbstractString)
    haskell_code = """
               $haskell_stmts
               """
    quote
        haskell_add_code!($haskell_code)
    end
end
export haskellcode

function haskellfunction(name::FnName, result::FnResult, arguments::AbstractVector{FnArg}, haskell_stmts::AbstractString)
    julia_code = quote
        function $(name.julia_name)($([:($(arg.julia_name)::$(arg.initial_julia_type)) for arg in arguments]...))
            res = ccall(($(name.cxx_name), $(name.cxx_library)), $(result.julia_type),
                        ($([arg.julia_type for arg in arguments if !arg.skip]...),),
                        $([arg.convert_from_initial(arg.julia_name) for arg in arguments if !arg.skip]...))
            return $(result.convert_to_final(:res))::$(result.final_julia_type)
        end
    end
    julia_code = clean_code(julia_code)
    simple_julia_code = simplify_code(julia_code)

    haskell_code = """
        {-
        $(string(simple_julia_code))
        -}
        $(name.cxx_name) :: $(join([arg.cxx_type for arg in arguments if !arg.skip], " -> ")) -> $(result.cxx_type)
        $(name.cxx_name) $(join([arg.cxx_name for arg in arguments if !arg.skip], " ")) =
            $(join(split(haskell_stmts, "\n"), "\n    "))
        foreign export ccall $(name.cxx_name) :: $(join([arg.cxx_type for arg in arguments if !arg.skip], " -> ")) -> $(result.cxx_type)
        """
    quote
        haskell_add_code!($haskell_code)
        $julia_code
    end
end
export haskellfunction


text:
using ArgDecomposition
using Test

@testset "ArgDecomposition.jl" begin

  for t ∈ (
    Val(4), Float32, 8.9, (a = 3, b = (3,2), d = Float64, e = (1f0, '2', 3.0, 4, 0x05)),
    (a = 3, b = (3,(2,21f0,3.0)), c = "hello, world!", d = Float64, e = (1f0, '2', 3.0, 4, 0x05))
  )

    ft = flatten_to_tuple(t)
    @test ft isa Tuple
    rt = reassemble_tuple(typeof(t), ft)
    @test rt === t
  end
end


text:
import ScikitLearnBase: transform

export transform

function transform(model::MatFacModel, A::AbstractMatrix;
                   inst_reg_weight::Real=1.0, feat_reg_weight::Real=1.0,
                   max_iter::Integer=100,
                   lr::Real=0.001, momentum::Real=0.5,
                   abs_tol::Real=1e-3, rel_tol::Real=1e-7,
                   loss_iter::Integer=10,
                   X_new::Union{Nothing,AbstractMatrix}=nothing,
                   new_inst_reg_mats::Union{Nothing,AbstractVector}=nothing,
                   K_opt_X::Union{Nothing,Integer}=nothing)

    # Setup
    iter = 0
    cur_loss = Inf
    new_loss = Inf

    K = size(model.X, 1)
    M_old = size(model.X, 2)
    M_new = size(A,1)
    M_comb = M_old + M_new
    N = size(A,2)

    if X_new == nothing
        X_new = 0.001*randn(K, M_new)
    end
    X_comb = hcat(model.X, X_new)

    println("X_COMB: ", size(X_comb))

    @assert size(X_comb,1) == size(model.Y,1)
    @assert size(X_new, 2) == size(A,1)

    # Distinguish between (1) K, the hidden dimension;
    #                     (2) the subsets 1:K_opt_X, 1:K_opt_Y of 1:K that is optimized;
    #                     (3) the subset of optimized dims that are regularized.
    # The regularized dim is determined by the lists of matrices
    # provided during MatFacModel construction.
    # K_opt_X, K_opt_Y have default value K.
    if K_opt_X == nothing
        K_opt_X = K
    end
    @assert K_opt_X <= K

    lr = Float32(lr)

    # Move data and model to the GPU.
    # Float16 is sufficiently precise for the data.
    X_comb_d = CuArray{Float32}(X_comb)
    Y_d = CuArray{Float32}(model.Y)
    A_d = CuArray{Float16}(A)
    XY_d = CUDA.zeros(Float16, M_new, N)

    println("X_COMB_D: ", size(X_comb_d))
    println("Y_D: ", size(Y_d))
    println("A_D: ", size(A_d))
    println("XY_D: ", size(XY_d))

    # Some frequently-used views of the X and Y arrays
    X_new_view = view(X_comb_d, :, (M_old+1):(M_old+M_new))
    X_opt_view = view(X_comb_d, 1:K_opt_X, (M_old+1):(M_old+M_new))
    Y_grad_view = view(Y_d, 1:K_opt_X, :)

    println("X_NEW_VIEW: ", size(X_new_view))
    println("X_OPT_VIEW: ", size(X_opt_view))
    println("Y_GRAD_VIEW: ", size(Y_grad_view))

    # Some bookkeeping for missing values.
    obs_mask = (!isnan).(A_d)
    missing_mask = (isnan.(A_d) .* Float16(0.5))
    # Convert NaNs in the data --> something benign so CuArray arithmetic works
    mask_func(x) = isnan(x) ? Float16(0.5) : Float16(x)
    map!(mask_func, A_d, A_d)

    # Scaling factors for the columns
    scales = [loss.scale for loss in model.losses]
    feature_scales = CuArray{Float32}(transpose(scales))

    # Bookkeeping for the loss functions
    ql_idx = CuVector{Int64}(findall(typeof.(model.losses) .== QuadLoss))
    ll_idx = CuVector{Int64}(findall(typeof.(model.losses) .== LogisticLoss))
    pl_idx = CuVector{Int64}(findall(typeof.(model.losses) .== PoissonLoss))

    XY_ql_view = view(XY_d, :, ql_idx)
    XY_ll_view = view(XY_d, :, ll_idx)
    XY_pl_view = view(XY_d, :, pl_idx)
    A_ql_view = view(A_d, :, ql_idx)
    A_ll_view = view(A_d, :, ll_idx)
    A_pl_view = view(A_d, :, pl_idx)

    # If regularizer matrices are not provided,
    # then let them be the identity
    if new_inst_reg_mats == nothing
        new_inst_reg_mats = [sparse(I, M_comb, M_comb) for i=1:K_opt_X]

    # Otherwise, make sure they're the right size
    # for the combined (old, new) X matrix.
    else
        @assert typeof(new_inst_reg_mats) <: AbstractVector
        for mat in new_inst_reg_mats
            @assert typeof(mat) <: AbstractMatrix
            @assert size(mat) == (M_comb, M_comb)
        end
    end
    # Convert regularizers to CuSparse matrices
    inst_reg_mats_d = [CUDA.CUSPARSE.CuSparseMatrixCSC{Float32}(mat .* inst_reg_weight) for mat in new_inst_reg_mats]

    # Arrays for holding gradients and velocities
    grad_X_temp = CUDA.zeros(Float32, (K_opt_X, M_comb))
    grad_X = CUDA.zeros(Float32, (K_opt_X, M_new))
    vel_X = CUDA.zeros(Float32, (K_opt_X, M_new))

    println("GRAD_X_TEMP: ", size(grad_X_temp))
    println("GRAD_X: ", size(grad_X))
    println("VEL_X: ", size(vel_X))

    while iter < max_iter

        ############################
        # Update X

        # take momentum "half-step"
        X_opt_view .+= (momentum.*vel_X)

        # compute gradient at half-step
        XY_d .= transpose(X_new_view)*Y_d
        XY_d .*= obs_mask
        XY_d .+= missing_mask
        compute_quadloss_delta!(XY_ql_view, A_ql_view)
        compute_logloss_delta!(XY_ll_view, A_ll_view)
        compute_poissonloss_delta!(XY_pl_view, A_pl_view)

        XY_d .*= feature_scales
        grad_X .= (Y_grad_view * transpose(XY_d)) ./ N

        grad_X_temp .= 0.0
        add_reg_grad!(grad_X_temp, X_comb_d, inst_reg_mats_d)
        grad_X .+= grad_X_temp[:,(M_old+1):(M_old+M_new)]

        # scale gradient by learning rate
        grad_X .*= lr

        # update velocity
        vel_X .*= momentum
        vel_X .-= grad_X

        # complete X update
        X_opt_view .-= grad_X

        iter += 1
        print_str = "Iteration: $iter"

        ############################
        # Every so often, compute the loss
        if (iter % loss_iter == 0)
            XY_d .= transpose(X_new_view)*Y_d
            XY_d .*= obs_mask

            new_loss = compute_quadloss!(XY_ql_view, A_ql_view)
            new_loss += compute_logloss!(XY_ll_view, A_ll_view)
            new_loss += compute_poissonloss!(XY_pl_view, A_pl_view)
            new_loss += compute_reg_loss(X_comb_d, inst_reg_mats_d)
            print_str = string(print_str, "\tLoss: $new_loss")
            println(print_str)

            if abs(new_loss - cur_loss) < abs_tol
                println(string("Absolute change <", abs_tol, ". Terminating."))
                break
            end
            if abs((new_loss - cur_loss)/cur_loss) < rel_tol
                println(string("Relative change <", rel_tol, ". Terminating."))
                break
            end
            cur_loss = new_loss
        end

    end # while

    # Move transformed X back to CPU
    X_new = Array{Float32}(X_new_view)

    return X_new
end


text:
const _residue_labels = map(string, reverse!(res"ARNDCQEGHILKMFPSTWYV-"))

@recipe function plot(msa::AbstractMatrix{Residue})
    seriestype --> :heatmap
    yflip --> true
    grid --> false
    foreground_color_border --> nothing
    foreground_color_axis --> nothing
    linewidth --> 0
    zdiscrete_values --> _residue_labels
    nseq = nsequences(msa)
    if nseq > 20
        step = div(nseq,20)
        yticks --> (1:step:nseq,names(msa))
        html_output_format :=  :png
    end
    1:ncolumns(msa), names(msa), map(string, msa)
end


text:
using NURBS.Bases

bs = BSplineBasis(0, 1, 1, 4)
pts = linspace(0, 1, 100000)

bs(pts)
@time bs(pts)

@profile bs(pts)
Profile.print()

Profile.clear_malloc_data()
bs(pts)


text:
mutable struct User
    id::Int
    name::String
    password::String
    email::String
    token::UInt32
end

const USERS = Dict{Int, User}()
const NEXT_USER_ID = Ref(0)

function getNextUserId()
    id = NEXT_USER_ID[]
    NEXT_USER_ID[] += 1
    return id
end

mutable struct Schedule
    id::Int
    uid::Int
    starttime::Int64
    endtime::Int64
    place::String
    prediction::Float64
    timestamp::Int64
end

const SCHEDULES = Dict{Int, Schedule}()
const NEXT_SCHEDULE_ID = Ref(0)

function getNextScheduleId()
    id = NEXT_SCHEDULE_ID[]
    NEXT_SCHEDULE_ID[] += 1
    return id
end


text:
module Example

f() = 42

end


text:
function solve()
    n, m, k = [parse(Int, x) for x in split(readline())]
    as = [parse(Int, x) for x in split(readline())]
    bs = [parse(Int, x) for x in split(readline())]

    as = cat(0, cumsum(as), dims=1)
    bs = cat(0, cumsum(bs), dims=1)

    ret = 0
    bi = m+1

    for ai=1:n+1
        while bi > 0 && as[ai] + bs[bi] > k
            bi -= 1
        end
        ret = max(ret, ai + bi - 2)
        if bi == 0
            break
        end
    end
    ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
PYDSTOOL_CURVE_CLASSES = Set(["EP-C","LP-C","H-C1","H-C2","FP-C","LC-C"])

CONT_BIF_POINTS = ["B", "SP"]
EQUILIBRIUM_BIF_POINTS = ["BP", "LP", "H"]
FOLD_BIF_POINTS = ["BT", "ZH", "CP"]
HOPF_BIF_POINTS = ["BT", "ZH", "GH", "DH"]
FIXEDPOINT_BIF_POINTS = ["BP", "PD", "LPC", "NS"]
FOLD_MAP_BIF_POINTS = ["CP"]
LIMITCYCLE_BIF_POINTS = ["PD", "LPC", "NS"]
OTHER_SPECIAL_POINTS = ["RG", "UZ", "P", "MX", "B"]

ALL_POINT_TYPES = union([CONT_BIF_POINTS;EQUILIBRIUM_BIF_POINTS;FOLD_BIF_POINTS;
                   HOPF_BIF_POINTS;FIXEDPOINT_BIF_POINTS;FOLD_MAP_BIF_POINTS;
                   LIMITCYCLE_BIF_POINTS;OTHER_SPECIAL_POINTS])

set_name(dsargs,name::String) = (dsargs[:name] = name; nothing)

set_ics(dsargs,icdict::PyDict) = (dsargs[:ics] = icdict; nothing)
set_ics(dsargs,icdict::Dict{String,T}) where {T} = set_ics(dsargs,PyDict(icdict))

set_pars(dsargs,pardict::PyDict) = (dsargs[:pars] = pardict; nothing)
set_pars(dsargs,pardict::Dict{String,T}) where {T} = set_pars(dsargs,PyDict(pardict))

set_vars(dsargs,vardict::PyDict) = (dsargs[:varspecs] = vardict; nothing)
set_vars(dsargs,vardict::Dict{String,T}) where {T} = set_vars(dsargs,PyDict(vardict))

set_fnspecs(dsargs,specsdict::PyDict) = (dsargs[:fnspecs] = specsdict; nothing)
set_fnspecs(dsargs,specsdict::Dict{String,T}) where {T} = set_vars(dsargs,PyDict(specsdict))

set_tdata(dsargs,tdata) = (dsargs[:tdata] = tdata; nothing)
set_tdomain(dsargs,tdomain) = (dsargs[:tdomain] = tdomain; nothing)


text:
"""
    AbstractBridgeOptimizer

A bridge optimizer applies given constraint bridges to a given optimizer thus
extending the types of supported constraints. The attributes of the inner
optimizer are automatically transformed to make the bridges transparent, e.g.
the variables and constraints created by the bridges are hidden.

By convention, the inner optimizer should be stored in a `model` field and
the dictionary mapping constraint indices to bridges should be stored in a
`bridges` field. If a bridge optimizer deviates from these conventions, it
should implement the functions `MOI.optimize!` and `bridge` respectively.
"""
abstract type AbstractBridgeOptimizer <: MOI.AbstractOptimizer end

# AbstractBridgeOptimizer interface

"""
    is_bridged(b::AbstractBridgeOptimizer, F::Type{<:MOI.AbstractFunction},
              S::Type{<:MOI.AbstractSet})::Bool

Return a `Bool` indicating whether `b` tries to bridge `F`-in-`S` constraints
instead of passing it as is to its internal model.
"""
function is_bridged end
# Syntactic sugar
function is_bridged(b::AbstractBridgeOptimizer, ::Type{CI{F, S}}) where {F, S}
    return is_bridged(b, F, S)
end
# We don't bridge variables.
is_bridged(b::AbstractBridgeOptimizer, ::Type{VI}) = false

"""
    supports_bridging_constraint(b::AbstractBridgeOptimizer,
                               F::Type{<:MOI.AbstractFunction},
                               S::Type{<:MOI.AbstractSet})::Bool

Return a `Bool` indicating whether `b` supports bridging `F`-in-`S` constraints.
"""
function supports_bridging_constraint(::AbstractBridgeOptimizer,
                                      ::Type{<:MOI.AbstractFunction},
                                      ::Type{<:MOI.AbstractSet})
    return false
end

"""
    bridge_type(b::AbstractBridgeOptimizer,
                F::Type{<:MOI.AbstractFunction},
                S::Type{<:MOI.AbstractSet})

Return the `AbstractBridge` type to be used to bridge `F`-in-`S` constraints.
This function should only be called if `is_bridged(b, F, S)`.
"""
function bridge_type end

function concrete_bridge_type(b::AbstractBridgeOptimizer,
                              F::Type{<:MOI.AbstractFunction},
                              S::Type{<:MOI.AbstractSet})
    return concrete_bridge_type(bridge_type(b, F, S), F, S)
end

"""
    bridge(b::AbstractBridgeOptimizer, ci::CI)

Return the `AbstractBridge` used to bridge the constraint with index `ci`.
"""
bridge(b::AbstractBridgeOptimizer, ci::CI) = b.bridges[ci.value]
# By convention, they should be stored in a `bridges` field using a
# dictionary-like object.
function bridge(b::AbstractBridgeOptimizer,
                ci::CI{MOI.SingleVariable, S}) where S
    return b.single_variable_constraints[(ci.value, S)]
end

# Implementation of the MOI interface for AbstractBridgeOptimizer

MOI.optimize!(b::AbstractBridgeOptimizer) = MOI.optimize!(b.model)
# By convention, the model should be stored in a `model` field

function MOI.is_empty(b::AbstractBridgeOptimizer)
    return isempty(b.bridges) && isempty(b.constraint_types) &&
           MOI.is_empty(b.model) && isempty(b.single_variable_constraints)
end
function MOI.empty!(b::AbstractBridgeOptimizer)
    MOI.empty!(b.model)
    empty!(b.bridges)
    empty!(b.constraint_types)
    empty!(b.single_variable_constraints)
end
function MOI.supports(b::AbstractBridgeOptimizer,
                      attr::Union{MOI.AbstractModelAttribute,
                                  MOI.AbstractOptimizerAttribute})
    return MOI.supports(b.model, attr)
end

function MOI.copy_to(mock::AbstractBridgeOptimizer, src::MOI.ModelLike; kws...)
    MOIU.automatic_copy_to(mock, src; kws...)
end
function MOIU.supports_default_copy_to(b::AbstractBridgeOptimizer,
                                       copy_names::Bool)
    return MOIU.supports_default_copy_to(b.model, copy_names)
end

# References
function _constraint_index(b::AbstractBridgeOptimizer, i::Integer)
    F, S = b.constraint_types[i]
    return MOI.ConstraintIndex{F, S}(i)
end
MOI.is_valid(b::AbstractBridgeOptimizer, vi::VI) = MOI.is_valid(b.model, vi)
function MOI.is_valid(b::AbstractBridgeOptimizer, ci::CI{F, S}) where {F, S}
    if is_bridged(b, typeof(ci))
        return 1 ≤ ci.value ≤ length(b.bridges) && b.bridges[ci.value] !== nothing &&
            (F, S) == b.constraint_types[ci.value]
    else
        return MOI.is_valid(b.model, ci)
    end
end
function MOI.is_valid(b::AbstractBridgeOptimizer,
                      ci::CI{MOI.SingleVariable, S}) where S
    if is_bridged(b, typeof(ci))
        return haskey(b.single_variable_constraints, (ci.value, S))
    else
        return MOI.is_valid(b.model, ci)
    end
end
function MOI.delete(b::AbstractBridgeOptimizer, vi::VI)
    # Delete all `MOI.SingleVariable` constraint of this variable
    for key in keys(b.single_variable_constraints)
        if key[1] == vi.value
            MOI.delete(b, MOI.ConstraintIndex{MOI.SingleVariable,
                                              key[2]}(vi.value))
        end
    end
    MOI.delete(b.model, vi)
end
function _remove_bridge(b, ci)
    b.bridges[ci.value] = nothing
end
function _remove_bridge(b, ci::MOI.ConstraintIndex{MOI.SingleVariable, S}) where S
    delete!(b.single_variable_constraints, (ci.value, S))
end
function MOI.delete(b::AbstractBridgeOptimizer, ci::CI)
    if is_bridged(b, typeof(ci))
        if !MOI.is_valid(b, ci)
            throw(MOI.InvalidIndex(ci))
        end
        MOI.delete(b, bridge(b, ci))
        _remove_bridge(b, ci)
    else
        MOI.delete(b.model, ci)
    end
end

# Attributes
# List of indices of all constraints, including those bridged
function get_all_including_bridged(
    b::AbstractBridgeOptimizer,
    attr::MOI.ListOfConstraintIndices{F, S}) where {F, S}
    if is_bridged(b, F, S)
        return MOI.ConstraintIndex{F, S}[MOI.ConstraintIndex{F, S}(i) for i in eachindex(b.bridges) if MOI.is_valid(b, MOI.ConstraintIndex{F, S}(i))]
    else
        return MOI.get(b.model, attr)
    end
end
function get_all_including_bridged(
    b::AbstractBridgeOptimizer,
    attr::MOI.ListOfConstraintIndices{MOI.SingleVariable, S}) where S
    F = MOI.SingleVariable
    if is_bridged(b, F, S)
        return MOI.ConstraintIndex{F, S}[MOI.ConstraintIndex{F, S}(key[1]) for key in keys(b.single_variable_constraints) if key[2] == S]
    else
        return MOI.get(b.model, attr)
    end
end
# Remove constraints bridged by `bridge` from `list`
function _remove_bridged(list, bridge, attr)
    for c in MOI.get(bridge, attr)
        i = something(findfirst(isequal(c), list), 0)
        if !iszero(i)
            MOI.deleteat!(list, i)
        end
    end
end
function MOI.get(b::AbstractBridgeOptimizer,
                 attr::MOI.ListOfConstraintIndices{F, S}) where {F, S}
    list = get_all_including_bridged(b, attr)
    for bridge in b.bridges
        if bridge !== nothing
            _remove_bridged(list, bridge, attr)
        end
    end
    for bridge in values(b.single_variable_constraints)
        _remove_bridged(list, bridge, attr)
    end
    return list
end
function _number_of(b::AbstractBridgeOptimizer, model::MOI.ModelLike,
                   attr::Union{MOI.NumberOfConstraints, MOI.NumberOfVariables})
    s = MOI.get(model, attr)
    for bridge in b.bridges
        if bridge !== nothing
            s -= MOI.get(bridge, attr)
        end
    end
    for bridge in values(b.single_variable_constraints)
        s -= MOI.get(bridge, attr)
    end
    return s
end
function MOI.get(b::AbstractBridgeOptimizer, attr::MOI.NumberOfVariables)
    return _number_of(b, b.model, attr)
end

# Number of all constraints, including those bridged
function get_all_including_bridged(
    b::AbstractBridgeOptimizer,
    attr::MOI.NumberOfConstraints{F, S}) where {F, S}
    if is_bridged(b, F, S)
        return count(i -> MOI.is_valid(b, MOI.ConstraintIndex{F, S}(i)), eachindex(b.bridges))
    else
        return MOI.get(b.model, attr)
    end
end
function get_all_including_bridged(
    b::AbstractBridgeOptimizer,
    attr::MOI.NumberOfConstraints{MOI.SingleVariable, S}) where S
    if is_bridged(b, MOI.SingleVariable, S)
        return count(key -> key[2] == S, keys(b.single_variable_constraints))
    else
        return MOI.get(b.model, attr)
    end
end
function MOI.get(b::AbstractBridgeOptimizer,
                 attr::MOI.NumberOfConstraints{F, S}) where {F, S}
    s = get_all_including_bridged(b, attr)
    # The constraints counted in `s` may have been added by bridges
    for bridge in b.bridges
        if bridge !== nothing
            s -= MOI.get(bridge, attr)
        end
    end
    for bridge in values(b.single_variable_constraints)
        s -= MOI.get(bridge, attr)
    end
    return s
end
function MOI.get(b::AbstractBridgeOptimizer, attr::MOI.ListOfConstraints)
    list_of_types = MOI.get(b.model, attr)
    list_of_bridged_types = Set{Tuple{DataType, DataType}}()
    for i in eachindex(b.bridges)
        if b.bridges[i] !== nothing
            push!(list_of_bridged_types, b.constraint_types[i])
        end
    end
    for key in keys(b.single_variable_constraints)
        push!(list_of_bridged_types, (MOI.SingleVariable, key[2]))
    end
    list_of_types = [list_of_types; collect(list_of_bridged_types)]
    # Some constraint types show up in `list_of_types` including when all the
    # constraints of that type have been created by bridges and not by the user.
    # The code in `NumberOfConstraints` takes care of removing these constraints
    # from the counter so we can rely on it to remove these constraint types.
    types_to_remove = findall(iszero.(
        map(FS -> MOI.get(b, MOI.NumberOfConstraints{FS...}()), list_of_types)))
    deleteat!(list_of_types, types_to_remove)
    return list_of_types
end

# Model an optimizer attributes
function MOI.get(b::AbstractBridgeOptimizer,
                 attr::Union{MOI.AbstractModelAttribute,
                             MOI.AbstractOptimizerAttribute})
    return MOI.get(b.model, attr)
end
function MOI.set(b::AbstractBridgeOptimizer,
                  attr::Union{MOI.AbstractModelAttribute,
                              MOI.AbstractOptimizerAttribute},
                  value)
    return MOI.set(b.model, attr, value)
end

# Variable attributes
function MOI.get(b::AbstractBridgeOptimizer,
                 attr::MOI.AbstractVariableAttribute,
                 index::VI)
    return MOI.get(b.model, attr, index)
end
function MOI.get(b::AbstractBridgeOptimizer,
                 attr::MOI.AbstractVariableAttribute,
                 indices::Vector{VI})
    return MOI.get(b.model, attr, indices)
end
function MOI.supports(b::AbstractBridgeOptimizer,
                      attr::MOI.AbstractVariableAttribute,
                      IndexType::Type{<:MOI.Index})
    return MOI.supports(b.model, attr, IndexType)
end
function MOI.set(b::AbstractBridgeOptimizer,
                  attr::MOI.AbstractVariableAttribute,
                  index::MOI.Index, value)
    return MOI.set(b.model, attr, index, value)
end
function MOI.set(b::AbstractBridgeOptimizer,
                  attr::MOI.AbstractVariableAttribute,
                  indices::Vector{<:MOI.Index}, values::Vector)
    return MOI.set(b.model, attr, indices, values)
end

# Constraint attributes
function MOI.get(b::AbstractBridgeOptimizer,
                 attr::MOI.AbstractConstraintAttribute, ci::CI)
    if is_bridged(b, typeof(ci))
        return MOI.get(b, attr, bridge(b, ci))
    else
        return MOI.get(b.model, attr, ci)
    end
end
function MOI.supports(b::AbstractBridgeOptimizer,
                      attr::MOI.AbstractConstraintAttribute,
                      IndexType::Type{CI{F, S}}) where {F, S}
    if is_bridged(b, IndexType)
        return MOI.supports(b, attr, concrete_bridge_type(b, F, S))
    else
        return MOI.supports(b.model, attr, IndexType)
    end
end

function MOI.set(b::AbstractBridgeOptimizer,
                 attr::MOI.AbstractConstraintAttribute,
                 index::CI, value)
    if is_bridged(b, typeof(index))
        return MOI.set(b, attr, bridge(b, index), value)
    else
        return MOI.set(b.model, attr, index, value)
    end
end
## Getting and Setting names
function MOI.get(b::AbstractBridgeOptimizer, attr::MOI.ConstraintName,
                 constraint_index::CI)
    if is_bridged(b, typeof(constraint_index))
        return get(b.con_to_name, constraint_index, MOIU.EMPTYSTRING)
    else
        return MOI.get(b.model, attr, constraint_index)
    end
end

function MOI.supports(b::AbstractBridgeOptimizer, attr::MOI.ConstraintName,
                      Index::Type{CI{F,S}}) where {F,S}
    if is_bridged(b, Index)
        return true
    else
        return MOI.supports(b.model, attr, Index)
    end
end
function MOI.set(b::AbstractBridgeOptimizer, attr::MOI.ConstraintName,
                  constraint_index::CI, name::String)
    if is_bridged(b, typeof(constraint_index))
        b.con_to_name[constraint_index] = name
        b.name_to_con = nothing # Invalidate the name map.
    else
        MOI.set(b.model, attr, constraint_index, name)
    end
end

# Name (similar to `UniversalFallback`)
function MOI.get(b::AbstractBridgeOptimizer, IdxT::Type{MOI.VariableIndex},
                 name::String)
    return MOI.get(b.model, IdxT, name)
end

function MOI.get(b::AbstractBridgeOptimizer, IdxT::Type{<:MOI.ConstraintIndex},
                 name::String)
    if b.name_to_con === nothing
        b.name_to_con = MOIU.build_name_to_con_map(b.con_to_name)
    end

    if is_bridged(b, IdxT)
        ci = get(b.name_to_con, name, nothing)
        if ci == CI{Nothing, Nothing}(-1)
            error("Multiple constraints have the name $name.")
        elseif ci isa IdxT
            return ci
        else
            return nothing
        end
    else
        ci = MOI.get(b.model, IdxT, name)
        if ci !== nothing && haskey(b.name_to_con, name)
            error("Multiple constraints have the name $name.")
        end
        return ci
    end
end

# We have no information as to whether the constraint is in the bridge or the
# model. Therefore, we try the model first, and then the bridge if that fails.
function MOI.get(b::AbstractBridgeOptimizer, IdxT::Type{CI},
                 name::String)
    if b.name_to_con === nothing
        b.name_to_con = MOIU.build_name_to_con_map(b.con_to_name)
    end

    ci = MOI.get(b.model, IdxT, name)
    if ci === nothing
        b_ci = get(b.name_to_con, name, nothing)
        if b_ci == CI{Nothing, Nothing}(-1)
            error("Multiple constraints have the name $name.")
        else
            return b_ci
        end
    else
        if haskey(b.name_to_con, name)
            error("Multiple constraints have the name $name.")
        end

        return ci
    end
end

# Constraints
function MOI.supports_constraint(b::AbstractBridgeOptimizer,
                                F::Type{<:MOI.AbstractFunction},
                                S::Type{<:MOI.AbstractSet})
    if is_bridged(b, F, S)
        return supports_bridging_constraint(b, F, S)
    else
        return MOI.supports_constraint(b.model, F, S)
    end
end
function store_bridge(b::AbstractBridgeOptimizer, func::MOI.SingleVariable,
                      set::MOI.AbstractSet, bridge)
    b.single_variable_constraints[(func.variable.value, typeof(set))] = bridge
    return MOI.ConstraintIndex{MOI.SingleVariable, typeof(set)}(func.variable.value)
end
function store_bridge(b::AbstractBridgeOptimizer, func::MOI.AbstractFunction,
                      set::MOI.AbstractSet, bridge)
    push!(b.bridges, bridge)
    push!(b.constraint_types, (typeof(func), typeof(set)))
    return MOI.ConstraintIndex{typeof(func), typeof(set)}(length(b.bridges))
end
function MOI.add_constraint(b::AbstractBridgeOptimizer, f::MOI.AbstractFunction,
                            s::MOI.AbstractSet)
    if is_bridged(b, typeof(f), typeof(s))
        # We compute `BridgeType` first as `concrete_bridge_type` calls
        # `bridge_type` which might throw an `UnsupportedConstraint` error in
        # which case, we do not want any modification to have been done
        BridgeType = concrete_bridge_type(b, typeof(f), typeof(s))
        # `add_constraint` might throw an `UnsupportedConstraint` but no
        # modification has been done in the previous line
        return store_bridge(b, f, s, bridge_constraint(BridgeType, b, f, s))
    else
        return MOI.add_constraint(b.model, f, s)
    end
end
function MOI.add_constraints(b::AbstractBridgeOptimizer, f::Vector{F},
                             s::Vector{S}) where { F <: MOI.AbstractFunction,
                             S <: MOI.AbstractSet}
    if is_bridged(b, F, S)
        return MOI.add_constraint.(b, f, s)
    else
        return MOI.add_constraints(b.model, f, s)
    end
end
function MOI.modify(b::AbstractBridgeOptimizer, ci::CI,
                     change::MOI.AbstractFunctionModification)
    if is_bridged(b, typeof(ci))
        MOI.modify(b, bridge(b, ci), change)
    else
        MOI.modify(b.model, ci, change)
    end
end

# Objective
function MOI.modify(b::AbstractBridgeOptimizer, obj::MOI.ObjectiveFunction,
                     change::MOI.AbstractFunctionModification)
    MOI.modify(b.model, obj, change)
end

# Variables
MOI.add_variable(b::AbstractBridgeOptimizer) = MOI.add_variable(b.model)
MOI.add_variables(b::AbstractBridgeOptimizer, n) = MOI.add_variables(b.model, n)

# TODO add transform


text:
function _is_cell_done(cell)
    if cell.running_disabled
        return true
    else
        return !cell.queued && !cell.running
    end
end

"""
    _is_notebook_done(notebook::Notebook)

Return whether all cells in the `notebook` have executed.
This method is more reliable than using `notebook.executetoken` because Pluto.jl drops that lock also after installing packages.
"""
function _notebook_done(notebook::Notebook)
    cells = [last(elem) for elem in notebook.cells_dict]
    return all(_is_cell_done, cells)
end

"""
    parallel_build!(
        dir,
        files;
        print_log=true,
        session=ServerSession()
    )

Build HTML files in parallel and write output to files with a ".html" extension.
This can be useful to speed up the build locally or in CI.
"""
function parallel_build!(
        dir,
        files;
        print_log=true,
        session=ServerSession()
    )

    # Start all the notebooks in parallel with async enabled.
    # This way, Pluto handles concurrency.
    notebooks = map(files) do in_file
        in_path = joinpath(dir, in_file)
        @assert isfile(in_path) "Expected .jl file at $in_path"

        @info "Starting evaluation of Pluto notebook at $in_file"
        notebook = SessionActions.open(session, in_path; run_async=true)
        return notebook
    end

    for (in_file, notebook) in zip(files, notebooks)
        while !_notebook_done(notebook)
            sleep(1)
        end

        without_extension, _ = splitext(in_file)
        out_file = "$(without_extension).html"
        out_path = joinpath(dir, out_file)

        html = notebook2html(notebook)
        SessionActions.shutdown(session, notebook)
        write(out_path, html)
    end

    return nothing
end

function parallel_build!(dir; print_log=true)
    files = filter(endswith(".jl"), readdir(dir))
    parallel_build!(dir, files; print_log)
    return nothing
end


text:
"""
     AbstractMutableComp

Abstract base type for computation units (read layers) which may mutate.
"""
abstract type AbstractMutableComp end
Base.Broadcast.broadcastable(m::AbstractMutableComp) = Ref(m)

# Generic forwarding methods. Just implement layer(t::Type) and wrapped(t::Type) to enable

# Not possible in julia <= 1.1. See #14919
# (m::AbstractMutableComp)(x...) = layer(m)(x...)
layertype(m::AbstractMutableComp) = layertype(layer(m))

NaiveNASlib.nin(m::AbstractMutableComp) = nin(wrapped(m))
NaiveNASlib.nout(m::AbstractMutableComp) = nout(wrapped(m))

function NaiveNASlib.Δsize!(m::AbstractMutableComp, inputs::AbstractVector, outputs::AbstractVector;kwargs...)
     NaiveNASlib.Δsize!(wrapped(m), inputs, outputs;kwargs...)
end

mutate_weights(m::AbstractMutableComp, w) = mutate_weights(wrapped(m), w)

function NaiveNASlib.compconstraint!(case, s::AbstractJuMPΔSizeStrategy, m::AbstractMutableComp, data)
     NaiveNASlib.compconstraint!(case, s, layertype(m), data)
end

function NaiveNASlib.defaultutility(m::AbstractMutableComp)
    util = neuronutility_safe(m)
    return util === 1 ? NaiveNASlib.NoDefaultUtility() : util
end

"""
    MutableLayer

Wraps a layer in order to allow for mutation as layer structures are typically immutable
"""
mutable struct MutableLayer <: AbstractMutableComp
    layer # Can't specialize as this might mutate into something else, e.g. Dense{Float32} -> Dense{Float64} thanks to mutable functor. Try to make functor rebuild everything?
end
(m::MutableLayer)(x...) = layer(m)(x...)
wrapped(m::MutableLayer) = m.layer
layer(m::MutableLayer) = wrapped(m)
layertype(m::MutableLayer) = layertype(layer(m))

@functor MutableLayer

function NaiveNASlib.Δsize!(m::MutableLayer, inputs::AbstractVector, outputs::AbstractVector; insert=neuroninsert, kwargs...)
    @assert length(inputs) == 1 "Only one input per layer!"
    mutate(m; inputs=inputs[1], outputs=outputs, insert=insert, kwargs...)
    nothing
end

mutate_weights(m::MutableLayer, w) = mutate(layertype(m), m, other=w)

function mutate(m::MutableLayer; inputs, outputs, other = l -> (), insert=neuroninsert)
    if inputs === missing
        mutate(layertype(m), m; outputs, other, insert)
    else
        mutate(layertype(m), m; inputs, outputs, other, insert)
    end
end

function mutate(lt::FluxParLayer, m::MutableLayer; inputs=1:nin(m)[], outputs=1:nout(m), other= l -> (), insert=neuroninsert)
    l = layer(m)
    otherdims = other(l)
    w = select(weights(l), indim(l) => inputs, outdim(l) => outputs, otherdims...; newfun=insert(lt, WeightParam()))
    b = select(bias(l), 1 => outputs; newfun=insert(lt, BiasParam()))
    newlayer(m, w, b, otherpars(other, l))
end
otherpars(o, l) = ()

function mutate(lt::FluxDepthwiseConv{N}, m::MutableLayer; inputs=1:nin(m)[], outputs=1:nout(m), other= l -> (), insert=neuroninsert) where N
    l = layer(m)
    otherdims = other(l)

    ngroups = div(length(outputs), length(inputs))

    # inputs and outputs are coupled through the constraints (which hopefully were enforced) so we only need to consider outputs
    currsize =size(weights(l))
    wo = select(reshape(weights(l), currsize[1:N]...,:), N+1 => outputs, otherdims...; newfun=insert(lt, WeightParam()))
    newks = size(wo)[1:N]
    w = collect(reshape(wo, newks...,ngroups, :))
    b = select(bias(l), 1 => outputs; newfun=insert(lt, BiasParam()))
    newlayer(m, w, b, otherpars(other, l))
end

function mutate(lt::FluxRecurrent, m::MutableLayer; inputs=1:nin(m)[], outputs=1:nout(m), other=missing, insert=neuroninsert)
    l = layer(m)
    outputs_scaled = mapfoldl(vcat, 1:outscale(l)) do i
        offs = (i-1) * nout(l)
        return map(x -> x > 0 ? x + offs : x, outputs)
    end

    wi = select(weights(l), indim(l) => inputs, outdim(l) => outputs_scaled, newfun=insert(lt, WeightParam()))
    wh = select(hiddenweights(l), 2 => outputs, 1 => outputs_scaled, newfun=insert(lt, RecurrentWeightParam()))
    b = select(bias(l), 1 => outputs_scaled, newfun=insert(lt, BiasParam()))
    mutate_recurrent_state(lt, m, outputs, wi, wh, b, insert)
end

function mutate_recurrent_state(lt::FluxRecurrent, m::MutableLayer, outputs, wi, wh, b, insert)
    l = layer(m)
    state0 = select(hiddenstate(l), 1 => outputs, newfun=insert(lt, RecurrentState()))
    s = select(state(l), 1 => outputs; newfun=insert(lt, RecurrentState()))

    cellnew = setproperties(layer(m).cell, (Wi=wi, Wh=wh, b = b, state0 = state0))
    lnew = setproperties(layer(m), (cell=cellnew, state = s))
    m.layer = lnew
end

function mutate_recurrent_state(lt::FluxLstm, m::MutableLayer, outputs, wi, wh, b, insert)
    l = layer(m)
    s0curr, scurr = hiddenstate(l), state(l)
    s0 = select.(s0curr, repeat([1 => outputs], length(s0curr)); newfun=insert(lt, RecurrentState()))
    s = select.(scurr, repeat([1 => outputs], length(scurr)); newfun=insert(lt, RecurrentState()))

    cellnew = setproperties(layer(m).cell, (Wi=wi, Wh=wh, b = b, state0 = Tuple(s0)))
    lnew = setproperties(layer(m), (cell=cellnew, state = Tuple(s)))
    m.layer = lnew

    return (;state0 = Tuple(s0)), Tuple(s)
end

function mutate(t::FluxParInvLayer, m::MutableLayer; inputs=missing, outputs=missing, other=missing, insert=neuroninsert)
    @assert any(ismissing.((inputs, outputs))) || inputs == outputs "Try to mutate $inputs and $outputs for invariant layer $(m)!"
    ismissing(inputs) || return mutate(t, m, inputs; insert=insert)
    ismissing(outputs) || return mutate(t, m, outputs; insert=insert)
end

function mutate(lt::FluxDiagonal, m::MutableLayer, inds; insert=neuroninsert)
    l = layer(m)
    w = select(weights(l), 1 => inds, newfun=insert(lt, WeightParam()))
    b = select(bias(l), 1 => inds; newfun=insert(lt, BiasParam()))
    newlayer(m, w, b)
end

function mutate(::FluxLayerNorm, m::MutableLayer, inds; insert=neuroninsert)
    # LayerNorm is only a wrapped Diagonal. Just mutate the Diagonal and make a new LayerNorm of it
    proxy = MutableLayer(layer(m).diag)
    mutate(proxy; inputs=inds, outputs=inds, other=l->(), insert=insert)

    updatelayer = layer(m)
    m.layer = @set updatelayer.diag = layer(proxy)
end

function mutate(lt::FluxParNorm, m::MutableLayer, inds; insert=neuroninsert)

    # Filter out the parameters which need to change and decide for each name (e.g. γ, β etc) what to do (typically insert 1 for scaling things and 0 for offset things)
    parselect(p::Pair) = parselect(p...)
    parselect(pname, x::AbstractArray) = select(x, 1 => inds; newfun = neuroninsert(lt, pname))
    parselect(pname, x) = x

    fs, re = Flux.functor(m.layer)
    newlayer = re(map(parselect, pairs(fs) |> collect))
    newlayer = @set newlayer.chs = length(inds)
    m.layer = newlayer
end

function mutate(lt::FluxGroupNorm, m::MutableLayer, inds; insert=neuroninsert)

    l = m.layer
    ngroups = l.G
    nchannels = length(inds)

    # Make a "best fit" to new group size, but don't allow 1 as group size
    # Step 1: Possible group sizes are all integers which divde the new number of channels evenly
    g_alts = filter(ga -> ga > 1, findall(r -> r == 0, nchannels .% (1:nchannels)))
    # Step 2: Select the size which is closest to the original number of groups
    ngroups = g_alts[argmin(abs.(g_alts .- ngroups))]

    nchannels_per_group = div(nchannels, ngroups)
    # Map indices to (new) groups
    inds_groups = ceil.(Int, map(ind -> ind > 0 ? ind / nchannels_per_group : ind, inds))
    # TODO: Select the most commonly occuring index in each column (except -1)
    inds_groups = reshape(inds_groups, nchannels_per_group, ngroups)[1,:]

    sizetoinds = Dict(nin(l)[] => inds, l.G => inds_groups)

    parselect(p::Pair) = parselect(p...)
    parselect(pname, x::AbstractArray) = select(x, 1 => sizetoinds[length(x)]; newfun = insert(lt, pname))
    parselect(pname, x) = x

    fs, re = Flux.functor(m.layer)
    newlayer = re(map(parselect, pairs(fs) |> collect))
    newlayer.G = ngroups
    newlayer = @set newlayer.chs = length(inds)
    m.layer = newlayer
    m.layer.G = ngroups
end

newlayer(m::MutableLayer, w, b, other=nothing) = m.layer = newlayer(layertype(m), m, w, b, other)

newlayer(::FluxDense, m::MutableLayer, w, b, other) = Dense(w, b, deepcopy(layer(m).σ))
newlayer(::FluxConvolutional, m::MutableLayer, w, b, other) = setproperties(layer(m), (weight=w, bias=b, σ=deepcopy(layer(m).σ), other...))
newlayer(::FluxDiagonal, m::MutableLayer, w, b, other) = Flux.Diagonal(w, b)

"""
    LazyMutable
    LazyMutable(m::AbstractMutableComp)

Lazy version of MutableLayer in the sense that it does not perform any mutations
until invoked to perform a computation.

This reduces the need to garbage collect when multiple mutations might be applied
to a vertex before evaluating the model.

Also useable for factory-like designs where the actual layers of a computation graph
are not instantiated until the graph is used.

# Examples
```jldoctest
julia> using NaiveNASflux, Flux

julia> struct DenseConfig end

julia> lazy = LazyMutable(DenseConfig(), 2, 3);

julia> layer(lazy)
DenseConfig()

julia> function NaiveNASflux.dispatch!(m::LazyMutable, ::DenseConfig, x)
       m.mutable = Dense(nin(m)[1], nout(m), relu)
       return m.mutable(x)
       end;

julia> lazy(ones(Float32, 2, 5)) |> size
(3, 5)

julia> layer(lazy)
Dense(2, 3, relu)   # 9 parameters
```
"""
mutable struct LazyMutable <: AbstractMutableComp
    mutable # May change at any time (see example)
    inputs::Vector{Vector{Int}}
    outputs::Vector{Int}
    other # May change at any time
    insert # May change at any time
end
LazyMutable(m::AbstractMutableComp) = LazyMutable(m, nin(m), nout(m))
LazyMutable(m, nin::Integer, nout::Integer) = LazyMutable(m, [nin], nout)
function LazyMutable(m, nins::AbstractVector{<:Integer}, nout::Integer)
    inputs = map(nin -> collect(Int, 1:nin), nins)
    outputs = collect(Int, 1:nout)
    LazyMutable(m, inputs, outputs, m -> (), neuroninsert)
end

function Base.show(io::IO, lm::LazyMutable)
    print(io, "LazyMutable(")
    show(io, lm.mutable)
    print(io, '[')
    for inpt in lm.inputs
        print(io, NaiveNASlib.compressed_string(inpt))
    end
    print(io, ']', NaiveNASlib.compressed_string(lm.outputs))
    show(io, lm.other)
    show(io, lm.insert)
end

wrapped(m::LazyMutable) = m.mutable
layer(m::LazyMutable) = layer(wrapped(m))

function Functors.functor(::Type{<:LazyMutable}, m)
    forcemutation(m)
    return (mutable=m.mutable,), y -> LazyMutable(y.mutable)
end

(m::LazyMutable)(x...) = dispatch!(m, m.mutable, x...)
dispatch!(::LazyMutable, m::AbstractMutableComp, x...) = m(x...)

NaiveNASlib.nin(m::LazyMutable) = length.(m.inputs)
NaiveNASlib.nout(m::LazyMutable) = length(m.outputs)

function NaiveNASlib.Δsize!(m::LazyMutable, inputs::AbstractVector, outputs::AbstractVector; insert=neuroninsert)
    (all((i,k)::Tuple -> ismissing(i) || i == 1:k, zip(inputs, nin(m))) && outputs == 1:nout(m)) && return
    m.insert = insert
    m.mutable = ResetLazyMutable(trigger_mutation(m.mutable))
    m.inputs = map(m.inputs, inputs) do currins, newins
        ismissing(newins) && return currins
        select(currins, 1 => newins, newfun = (args...) -> -1)
    end
    m.outputs = select(m.outputs, 1=>outputs, newfun = (args...) -> -1)
    nothing
end

function mutate_weights(m::LazyMutable, w)
    m.other == w && return

    m.mutable = ResetLazyMutable(trigger_mutation(m.mutable))
    m.other = w
end

function forcemutation(x) end
function forcemutation(::InputVertex) end
forcemutation(g::CompGraph) = forcemutation.(vertices(g::CompGraph))
forcemutation(v::AbstractVertex) = forcemutation(base(v))
forcemutation(v::CompVertex) = forcemutation(v.computation)
forcemutation(m::AbstractMutableComp) = forcemutation(NaiveNASflux.wrapped(m))
forcemutation(m::LazyMutable) = m(NoComp())

struct NoComp end
function (::MutableLayer)(::NoComp) end

"""
    MutationTriggered

Dispatches mutation for LazyMutable.
"""
struct MutationTriggered{T}
    wrapped::T
end
# Functionality is opt-in
trigger_mutation(m) = m
trigger_mutation(m::AbstractMutableComp) = MutationTriggered(m)

function dispatch!(lm::LazyMutable, m::MutationTriggered, x...)
    NaiveNASlib.Δsize!(m.wrapped, lm.inputs, lm.outputs; other=lm.other, insert=lm.insert)
    lm.mutable = m.wrapped
    return lm(x...)
end

layer(m::MutationTriggered) = layer(m.wrapped)
layertype(m::MutationTriggered) = layertype(layer(m))

@functor MutationTriggered

"""
    ResetLazyMutable

Reset a `LazyMutable` when dispatching.
"""
struct ResetLazyMutable{T}
    wrapped::T
end
ResetLazyMutable(r::ResetLazyMutable) = r

function dispatch!(lm::LazyMutable, m::ResetLazyMutable, x...)
    lm.mutable = m.wrapped
    output = lm(x...)
    lm.inputs = map(i -> collect(1:i), nin(lm))
    lm.outputs = 1:nout(lm)
    lm.other = m -> ()
    lm.insert = neuroninsert
    return output
end

layer(m::ResetLazyMutable) = layer(m.wrapped)
layertype(m::ResetLazyMutable) = layertype(layer(m))

@functor ResetLazyMutable

"""
    NoParams

Ignores size mutation.

Useful for layers which don't have parameters.
"""
struct NoParams{T}
    layer::T
end

(i::NoParams)(x...) = layer(i)(x...)
layer(i::NoParams) = i.layer
layertype(i::NoParams) = layertype(layer(i))

LazyMutable(m::NoParams) = m

function mutate_weights(::NoParams, w) end


text:
#===
ファイルの読み込み
===#
# include(".jl")

# export save

#===
共通型，関数
・全ての構造体はOptimizerの子型になる
===#

export SGD,Momentum,Adam,update!

#===SGD==============#
mutable struct SGD <: Optimizer
    learning_rate
    function SGD(;learning_rate=0.001)
        new(learning_rate)
    end
end
function fit!(this::SGD,model::Models)
    model.common.optimizer = this
end
function update!(this::SGD,model::Models)
    cnt = 1
    @inbounds for i in 1:length(model.common.layers)
        grads = model.common.layers[i].grads
        if grads === nothing
            continue
        end
        for j in 1:length(grads)
            grad = grads[j]
            model.common.layers[i].params[j] .-= this.learning_rate*grad
            cnt+=1
        end
    end
    nothing
end
function convert_to_cu!(this::SGD,model::Models)
    return nothing
end
function convert_to_array!(this::SGD,model::Models)
    return nothing
end

#===Momentum=========#
mutable struct Momentum <: Optimizer
    learning_rate
    vs #速度
    momentum
    function Momentum(;learning_rate=0.001,momentum=0.8)
        new(learning_rate,[],momentum)
    end
end
function fit!(this::Momentum,model::Models)
    @inbounds for i in 1:length(model.common.layers)
        params = model.common.layers[i].params
        if params === nothing
            continue
        end
        for j in 1:length(params)
            param = params[j]
            s = copy(param) * 0
            append!(this.vs,[s])
        end
    end
    model.common.optimizer = this
    nothing
end
function update!(this::Momentum,model::Models)
    cnt = 1
    @inbounds for i in 1:length(model.common.layers)
        grads = model.common.layers[i].grads
        if grads === nothing
            continue
        end
        for j in 1:length(grads)
            grad = grads[j]
            this.vs[cnt] = this.momentum * this.vs[cnt] - this.learning_rate * grad
            model.common.layers[i].params[j] .+= this.vs[cnt]
            cnt+=1
        end
    end
    nothing
end
function convert_to_cu!(this::Momentum,model::Models)
    this.vs = cu.(this.vs)
    return nothing
end
function convert_to_array!(this::Momentum,model::Models)
    this.vs = Array.(this.vs)
    return nothing
end

#====Adam============#
mutable struct Adam <: Optimizer
    learning_rate
    p1
    p2
    e
    vs #momentam
    ss #RMSProp
    function Adam(;learning_rate=0.001,p1=0.95,p2=0.99,e=10^(-12))
        new(learning_rate,p1,p2,e,[],[])
    end
end
function fit!(this::Adam,model::Models)
    @inbounds for i in 1:length(model.common.layers)
        params = model.common.layers[i].params
        if params === nothing
            continue
        end
        for j in 1:length(params)
            param = params[j]
            s = copy(param) * 0
            append!(this.vs,[s])
            append!(this.ss,[s])
        end
    end
    model.common.optimizer = this
    nothing
end
function update!(this::Adam,model::Models)
    cnt = 1
    @inbounds for i in 1:length(model.common.layers)
        grads = model.common.layers[i].grads
        if grads === nothing
            continue
        end
        for j in 1:length(grads)
            grad = grads[j]
            v = this.p1 .* this.vs[cnt] .+ (1-this.p1).*grad
            s = this.p2 .* this.ss[cnt] .+ (1-this.p2).*(grad .^2)
            model.common.layers[i].params[j] .-= this.learning_rate*(v./sqrt.(this.e .+s))
            this.vs[cnt] = v
            this.ss[cnt] = s
            cnt+=1
        end
    end
    nothing
end
function convert_to_cu!(this::Adam,model::Models)
    this.vs = cu.(this.vs)
    this.ss = cu.(this.ss)
    return nothing
end
function convert_to_array!(this::Adam,model::Models)
    this.vs = Array.(this.vs)
    this.ss = Array.(this.ss)
    return nothing
end

max_stars_repo_name: catmanium/Spike

max_issues_repo_name: catmanium/Spike

max_forks_repo_name: catmanium/Spike