text:
# Test: ODESystem built with ModelingToolkit + Ensemble Simulation with DifferentialEquations + Optimisation with DiffEqFlux
using DiffEqFlux, DifferentialEquations, Plots, Statistics, ModelingToolkit
ann = FastChain(FastDense(1,32,tanh), FastDense(32,32,tanh), FastDense(32,1))
dim_p = DiffEqFlux.paramlength(ann)

@variables t x1(t) x2(t)
@parameters p_NN[1:dim_p]
D = Differential(t)

deqs = [D(x1) ~ x2
        D(x2) ~ ann([t],p_NN)[1]^3]
sys = ODESystem(deqs, t, [x1, x2], p_NN; name = :sys)
# sys = ODESystem(deqs, iv, dvs, ps; )
# sys_simplified = structural_simplify(sys) # ---> not working

u0 = [x1 => -4f0,
      x2 => 0f0]
p_NN_0 = initial_params(ann)
p = [p_NN[i] => p_NN_0[i] for i in 1:dim_p]
tspan = (0.0f0, 8.0f0)
prob = ODEProblem(sys, u0, tspan, p; jac = true)
@time sol = solve(prob, Tsit5(), saveat = 0.1f0) # ---> not working

##--------
using DiffEqFlux, DifferentialEquations, Plots, Statistics, ModelingToolkit
dim_h = 32
tspan = (0.0f0,8.0f0)
ann = FastChain(FastDense(1,32,tanh), FastDense(32,32,tanh), FastDense(32,1))
θ = initial_params(ann)
function dxdt_(dx,x,p,t)
    x1, x2 = x
    dx[1] = x[2]
    dx[2] = ann([t],p)[1]^3
end
x0 = [-4f0, 0f0]
prob = ODEProblem(dxdt_,x0,tspan,θ)
sys = structural_simplify(modelingtoolkitize(prob))
prob_mtk = ODEProblem(sys, [], tspan, jac = true)

@time sol = solve(prob_mtk, Tsit5(), saveat = 0.1f0) # ---> working

#------------
ts = Float32.(collect(0.0:0.01:tspan[2]))
function predict_adjoint(θ)
  Array(solve(prob_mtk,Tsit5(),p=θ,saveat=ts,sensealg=InterpolatingAdjoint(autojacvec=ZygoteVJP())))
end
function loss_adjoint(θ)
  x = predict_adjoint(θ)
  mean(abs2,4.0 .- x[1,:]) + 2mean(abs2,x[2,:]) + mean(abs2,[first(ann([t],θ)) for t in ts])/10
end
l = loss_adjoint(θ)
cb = function (θ,l)
  println(l)
  # p = plot(solve(remake(prob,p=θ),Tsit5(),saveat=0.01),ylim=(-6,6),lw=3)
  # plot!(p,ts,[first(ann([t],θ)) for t in ts],label="u(t)",lw=3)
  # display(p)
  return false
end
# Display the ODE with the current parameter values.
cb(θ,l)
loss1 = loss_adjoint(θ)
res1 = DiffEqFlux.sciml_train(loss_adjoint, θ, ADAM(0.005), cb = cb,maxiters=100)
res2 = DiffEqFlux.sciml_train(loss_adjoint, res1.u,
                              BFGS(initial_stepnorm=0.01), cb = cb,maxiters=100,
                              allow_f_increases = false)

## ------------------------------------------------------------------------------------
using ModelingToolkit

function decay(;name)
  @parameters t a
  @variables x(t) f(t)
  D = Differential(t)
  ODESystem([
      D(x) ~ -a*x + f
    ];
    name=name)
end

@named decay1 = decay()
@named decay2 = decay()

@parameters t
D = Differential(t)
connected = compose(ODESystem([
                        decay2.f ~ decay1.x
                        D(decay1.f) ~ 0
                      ], t; name=:connected), decay1, decay2)

equations(connected)

#4-element Vector{Equation}:
# Differential(t)(decay1₊f(t)) ~ 0
# decay2₊f(t) ~ decay1₊x(t)
# Differential(t)(decay1₊x(t)) ~ decay1₊f(t) - (decay1₊a*(decay1₊x(t)))
# Differential(t)(decay2₊x(t)) ~ decay2₊f(t) - (decay2₊a*(decay2₊x(t)))

simplified_sys = structural_simplify(connected)

equations(simplified_sys)

#3-element Vector{Equation}:
# Differential(t)(decay1₊f(t)) ~ 0
# Differential(t)(decay1₊x(t)) ~ decay1₊f(t) - (decay1₊a*(decay1₊x(t)))
# Differential(t)(decay2₊x(t)) ~ decay1₊x(t) - (decay2₊a*(decay2₊x(t)))

x0 = [
  decay1.x => 1.0
  decay1.f => 0.0
  decay2.x => 1.0
]
p = [
  decay1.a => 0.1
  decay2.a => 0.2
]

using DifferentialEquations
prob = ODEProblem(simplified_sys, x0, (0.0, 100.0), p)
sol = solve(prob, Tsit5())
sol[decay2.f]

## ------------------------------------------------------------------------------------
using DifferentialEquations, DiffEqFlux
pa = [1.0]
u0 = [3.0]
θ = [u0;pa]

function model1(θ,ensemble)
  prob = ODEProblem((u, p, t) -> 1.01u .* p, [θ[1]], (0.0, 1.0), [θ[2]])

  function prob_func(prob, i, repeat)
    remake(prob, u0 = 0.5 .+ i/100 .* prob.u0)
  end

  ensemble_prob = EnsembleProblem(prob, prob_func = prob_func)
  sim = solve(ensemble_prob, Tsit5(), ensemble, saveat = 0.1, trajectories = 100)
end

# loss function
loss_serial(θ)   = sum(abs2,1.0.-Array(model1(θ,EnsembleSerial())))
loss_threaded(θ) = sum(abs2,1.0.-Array(model1(θ,EnsembleThreads())))

cb = function (θ,l) # callback function to observe training
  @show l
  false
end

opt = ADAM(0.1)
l1 = loss_serial(θ)
res_serial = DiffEqFlux.sciml_train(loss_serial, θ, opt; cb = cb, maxiters=100)
res_threads = DiffEqFlux.sciml_train(loss_threaded, θ, opt; cb = cb, maxiters=100)

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

using DifferentialEquations
prob = ODEProblem((u,p,t)->1.01*[u[2];u[1]],[0.5, 1.0],(0.0,1.0))
function prob_func(prob,i,repeat)
  remake(prob,u0=rand()*prob.u0)
end
ensemble_prob = EnsembleProblem(prob,prob_func=prob_func)
sim = solve(ensemble_prob,Tsit5(),EnsembleThreads(),trajectories=100)

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

using DifferentialEquations, Plots, GalacticOptim

function lotka_volterra!(du,u,p,t)
    rab, wol = u
    α,β,γ,δ=p
    du[1] = drab = α*rab - β*rab*wol
    du[2] = dwol = γ*rab*wol - δ*wol
    nothing
end

u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0]
prob = ODEProblem(lotka_volterra!,u0,tspan,p)
sol = solve(prob,saveat=0.1)
plot(sol)

dataset = Array(sol)
scatter!(sol.t,dataset')

tmp_prob = remake(prob, p=[1.2,0.8,2.5,0.8])
tmp_sol = solve(tmp_prob)
plot(tmp_sol)
scatter!(sol.t,dataset')

function loss(p)
  tmp_prob = remake(prob, p=p)
  tmp_sol = solve(tmp_prob,Tsit5(),saveat=0.1)
  if tmp_sol.retcode == :Success
    return sum(abs2,Array(tmp_sol) - dataset)
  else
    return 1E10
  end
end

# function loss(p)
#   tmp_prob = remake(prob, p=p)
#   tmp_sol = solve(tmp_prob,Tsit5(),saveat=0.1)
#   if size(tmp_sol) == size(dataset)
#     return sum(abs2,Array(tmp_sol) .- dataset)
#   else
#     return Inf
#   end
# end

using DiffEqFlux

pinit = [1.2,0.8,2.5,0.8]
res = DiffEqFlux.sciml_train(loss,pinit,ADAM(), maxiters = 100)

# res = DiffEqFlux.sciml_train(loss,pinit,BFGS(), maxiters = 100) ### errors!

#try Newton method of optimization
res = DiffEqFlux.sciml_train(loss,pinit,BFGS())

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

# compatibility of modelingtoolkitize and ensemble parellel simulation
using DifferentialEquations, ModelingToolkit
function rober(du,u,p,t)
  y₁,y₂,y₃ = u
  k₁,k₂,k₃ = p
  du[1] = -k₁*y₁+k₃*y₂*y₃
  du[2] =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
  du[3] =  k₂*y₂^2
  nothing
end

function prob_func(prob,i,repeat)
    p_i = collect(prob.p)
    p_i[1] = 1E-2 * i
    remake(prob,p=p_i)
    # remake(prob,u0=rand()*prob.u0)
end

prob = ODEProblem(rober,[1.0,0.0,0.0],(0.0,1e1),(0.04,3e7,1e4))
sys = modelingtoolkitize(prob)
prob_jac = ODEProblem(sys,[],(0.0,1e1),jac=true)

ensemble_prob = EnsembleProblem(prob,prob_func=prob_func)
ensemble_prob_jac = EnsembleProblem(prob_jac,prob_func=prob_func)

@time sim = solve(ensemble_prob,Tsit5(),EnsembleSerial(),trajectories=10)
@time sim_jac = solve(ensemble_prob_jac,Tsit5(),EnsembleSerial(),trajectories=10)

# conclusion 1) sim_jac is faster than sim if ensemble algorithm is EnsembleSerial().
# ---> sim_jac is even slower than sim when EnsembleThreads() is used.
# conclusion 2) sim_jac[i].prob.p is determined as intended by prob_func even though the modelingtoolkitized system is used.

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

using Flux

d = Dense(5,2)
p1 = d.weight, d.bias
p2 = Flux.params(d)
p1[1] == p2[1]
p2[1] == p2[2]

m = Chain(Dense(10,5), Dense(5,2))
x = rand(10)
# m(x) == m[2](m[1](x))

# loss(ŷ, y)

p, re = Flux.destructure(m)

m(x) - re(p)(x)

## ------------------------------------------------------------------------------------
using Flux

W = rand(2, 5)
b = rand(2)

predict(x) = (W * x) .+ b
loss(x, y) = sum((predict(x) .- y).^2)

x, y = rand(5), rand(2) # Dummy data
l = loss(x, y) # ~ 3

θ = params(W, b)
grads = gradient(() -> loss(x, y), θ)

using Flux.Optimise: update!

η = 0.1 # Learning Rate
for p in (W, b)
  update!(p, η * grads[p])
end

opt = Descent(0.1) # Gradient descent with learning rate 0.1

for p in (W, b)
  update!(opt, p, grads[p])
end

## ------------------------------------------------------------------------------------
# Comparison between the automatic differentiation results for ForwardDiffSensitivity and InterpolatingAdjoint
using DiffEqSensitivity, OrdinaryDiffEq, Zygote, LinearAlgebra, QuadGK

function fiip(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -p[3]*u[2] + p[4]*u[1]*u[2]
  du[3] = u[2]^2
end

p       = [1.5,1.0,3.0,1.0]
u0      = [1.0;1.0;0.0]
prob    = ODEProblem(fiip,u0,(0.0,10.0),p)
# y_sense = solve(prob,Tsit5(),u0=u0,p=p,saveat=0.1;sensealg = InterpolatingAdjoint(autojacvec=ZygoteVJP()))

# loss(u0,p) = sum(solve(prob,Tsit5(),u0=u0,p=p,saveat=0.1,abstol=1e-14,reltol=1e-14))
# loss_adjoint(u0,p) = sum(solve(prob,Tsit5(),u0=u0,p=p,saveat=0.1,abstol=1e-14,reltol=1e-14;sensealg = InterpolatingAdjoint()))
function loss(u0,p)
  y = Array(solve(prob,Tsit5(),u0=u0,p=p,saveat=0.1,abstol=1e-14,reltol=1e-14))
  return y[3,end]
end
function loss_adjoint(u0,p)
  y = Array(solve(prob,Tsit5(),u0=u0,p=p,saveat=0.1,abstol=1e-14,reltol=1e-14;sensealg = InterpolatingAdjoint()))
  return y[3,end]
end
function loss_discrete(u0,p)
  y = Array(solve(prob,Tsit5(),u0=u0,p=p,saveat=0.1,abstol=1e-14,reltol=1e-14;sensealg = InterpolatingAdjoint()))
  # return quadgk((t) -> (y(t)[2])^2, 0.0, 10.0)[1] # atol=1e-14,rtol=1e-10)  # this is incorrect and not working
  return 0.1*sum(y[2,:].^2)
end
function loss_adjoint_Bolza(u0,p)
  y = Array(solve(prob,Tsit5(),u0=u0,p=p,abstol=1e-14,reltol=1e-14;sensealg = InterpolatingAdjoint()))
  return y[3,end] + 100.0*y[1:2,end]'*y[1:2,end]
end

du01,dp1 = Zygote.gradient(loss,u0,p)
du02,dp2 = Zygote.gradient(loss_adjoint,u0,p)
du03,dp3 = Zygote.gradient(loss_discrete,u0,p)
du04,dp4 = Zygote.gradient(loss_adjoint_Bolza,u0,p)

# Comparison between the results of automatic differentiation and continuous adjoint sensitivity analysis method
y = solve(prob,Tsit5(),u0=u0,p=p,abstol=1e-14,reltol=1e-14;sensealg = InterpolatingAdjoint())
bwd_sol = solve(ODEAdjointProblem(y, InterpolatingAdjoint(), (out, x, p, t, i) -> (out[:] = [200.0*y[end][1:2]; 1]'), [10.0]), Tsit5(),
                # dense = false, save_everystep = false, save_start = false, # saveat = 0.1,
                abstol=1e-14,reltol=1e-14)
du05, dp5 = bwd_sol[end][1:length(u0)], -bwd_sol[end][1+length(u0):end]

isapprox(du04, du05)
isapprox(dp4, dp5)

## ------------------------------------------------------------------------------------
# Comparison between forward and adjoint sensitivity calculation results
using DiffEqSensitivity, OrdinaryDiffEq, Zygote

function f(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -p[3]*u[2] + u[1]*u[2]
end

p = [1.5,1.0,3.0]
u0 = [1.0;1.0]
prob = ODEProblem(f,u0,(0.0,10.0),p)
sol = solve(prob,Vern9(),abstol=1e-10,reltol=1e-10)

g(u,p,t) = (sum(u).^2) ./ 2 + p'*p / 2

function dgdu(out,u,p,t)
  out[1]= u[1] + u[2]
  out[2]= u[1] + u[2]
end

function dgdp_s(out,u,p,t)
  out[1] = p[1]
  out[2] = p[2]
  out[3] = p[3]
end

function dgdp_v1(out,u,p,t)
  out .= p
end

function dgdp_v2(out,u,p,t)
  out = p
end

function dgdp_v3(out,u,p,t)
  out = p'
end

res1 = adjoint_sensitivities(sol,Vern9(),g,nothing,abstol=1e-8, reltol=1e-8,iabstol=1e-8,ireltol=1e-8)                # this is incorrect
res2 = adjoint_sensitivities(sol,Vern9(),g,nothing,dgdu,abstol=1e-8, reltol=1e-8,iabstol=1e-8,ireltol=1e-8)           # this is incorrect
res3 = adjoint_sensitivities(sol,Vern9(),g,nothing,(dgdu,dgdp_s),abstol=1e-8, reltol=1e-8,iabstol=1e-8,ireltol=1e-8)  # this is correct
res4 = adjoint_sensitivities(sol,Vern9(),g,nothing,(dgdu,dgdp_v1),abstol=1e-8, reltol=1e-8,iabstol=1e-8,ireltol=1e-8) # this is correct
res5 = adjoint_sensitivities(sol,Vern9(),g,nothing,(dgdu,dgdp_v2),abstol=1e-8, reltol=1e-8,iabstol=1e-8,ireltol=1e-8) # this is incorrect
res6 = adjoint_sensitivities(sol,Vern9(),g,nothing,(dgdu,dgdp_v3),abstol=1e-8, reltol=1e-8,iabstol=1e-8,ireltol=1e-8) # this is incorrect

using ForwardDiff, Calculus
using QuadGK
function G(p)
  tmp_prob = remake(prob,p=p)
  sol = solve(tmp_prob,Vern9(),abstol=1e-14,reltol=1e-14)
  res,err = quadgk((t)-> (sum(sol(t)).^2)./2 + p'*p/2,0.0,10.0,atol=1e-14,rtol=1e-10)
  res
end
res7 = ForwardDiff.gradient(G,p)
res8 = Calculus.gradient(G,p)
# res9 = Zygote.gradient(G,p) # not working

## ------------------------------------------------------------------------------------
## ------------------------------------------------------------------------------------
## ------------------------------------------------------------------------------------
# Samuel Aintworth's code

import DiffEqBase
import DiffEqSensitivity:
    solve,
    ODEProblem,
    ODEAdjointProblem,
    InterpolatingAdjoint
import Zygote
import Statistics: mean

function extract_loss_and_xT(fwd_sol)
    fwd_sol[end][1], fwd_sol[end][2:end]
end

"""Returns a differentiable loss function that rolls out a policy in an
environment and calculates its cost."""
function ppg_goodies(dynamics, cost, policy, T)

    function aug_dynamics(z, policy_params, t)
        x = @view z[2:end]
        u = policy(x, t, policy_params)
        [cost(x, u); dynamics(x, u)]
    end

    function loss_pullback(x0, policy_params, solvealg, solve_kwargs)
        z0 = vcat(0.0, x0)
        fwd_sol = solve(
            ODEProblem(aug_dynamics, z0, (0, T), policy_params),
            solvealg,
            u0 = z0,
            p = policy_params;
            solve_kwargs...,
        )

        function _adjoint_solve(g_zT, sensealg; kwargs...)
            # See https://diffeq.sciml.ai/stable/analysis/sensitivity/#Syntax-1
            # and https://github.com/SciML/DiffEqSensitivity.jl/blob/master/src/local_sensitivity/sensitivity_interface.jl#L9.
            solve(
                ODEAdjointProblem(
                    fwd_sol,
                    sensealg,
                    (out, x, p, t, i) -> (out[:] = g_zT),
                    [T],
                ),
                solvealg;
                kwargs...,
            )
        end

        # This is the pullback using the augmented system and a discrete
        # gradient input at time T. Alternatively one could use the continuous
        # adjoints on the non-augmented system although this seems to be slower
        # and a less stable feature.

        function pullback(g_zT, sensealg::InterpolatingAdjoint)
            bwd_sol = _adjoint_solve(
                g_zT,
                sensealg,
                dense = false,
                save_everystep = false,
                save_start = false,
                # reltol = 1e-3,
                # abstol = 1e-3,
            )

            # The first z_dim elements of bwd_sol.u are the gradient wrt z0,
            # next however many are the gradient wrt policy_params.
            p = fwd_sol.prob.p
            l = p === nothing || p === DiffEqBase.NullParameters() ? 0 :
                length(fwd_sol.prob.p)
            g_x0 = bwd_sol[end][1:length(fwd_sol.prob.u0)]

            # We do exactly as many f calls as there are function calls in the
            # forward pass, and in the backward pass we don't need to call f,
            # but instead we call ∇f.
            (
                g = -bwd_sol[end][(1:l).+length(fwd_sol.prob.u0)],
                nf = 0,
                n∇ₓf = bwd_sol.destats.nf,
                n∇ᵤf = bwd_sol.destats.nf,
            )
        end

        fwd_sol, pullback
    end

    function ez_loss_and_grad(
        x0,
        policy_params,
        solvealg,
        sensealg;
        fwd_solve_kwargs = Dict(),
    )
        # @info "fwd"
        fwd_sol, vjp = loss_pullback(x0, policy_params, solvealg, fwd_solve_kwargs)
        # @info "bwd"
        bwd = vjp(vcat(1, zero(x0)), sensealg)
        loss, _ = extract_loss_and_xT(fwd_sol)
        # @info "fin"
        loss, bwd.g, (nf = fwd_sol.destats.nf + bwd.nf, n∇ₓf = bwd.n∇ₓf, n∇ᵤf = bwd.n∇ᵤf)
    end

    function _aggregate_batch_results(res)
        (
            mean(loss for (loss, _, _) in res),
            mean(g for (_, g, _) in res),
            (
                nf = sum(info.nf for (_, _, info) in res),
                n∇ₓf = sum(info.n∇ₓf for (_, _, info) in res),
                n∇ᵤf = sum(info.n∇ᵤf for (_, _, info) in res),
            ),
        )
    end

    function ez_loss_and_grad_many(
        x0_batch,
        policy_params,
        solvealg,
        sensealg;
        fwd_solve_kwargs = Dict(),
    )
        # Using tmap here gives a segfault. See https://github.com/tro3/ThreadPools.jl/issues/18.
        _aggregate_batch_results(
            map(x0_batch) do x0
                ez_loss_and_grad(
                    x0,
                    policy_params,
                    solvealg,
                    sensealg,
                    fwd_solve_kwargs = fwd_solve_kwargs,
                )
            end,
        )
    end

    (
        aug_dynamics = aug_dynamics,
        loss_pullback = loss_pullback,
        ez_loss_and_grad = ez_loss_and_grad,
        ez_loss_and_grad_many = ez_loss_and_grad_many,
        ez_euler_bptt = ez_euler_bptt,
        ez_euler_loss_and_grad_many = ez_euler_loss_and_grad_many,
    )
end

function policy_dynamics!(dx, x, policy_params, t)
    u = policy(x, policy_params)
    dx .= dynamics(x, u)
end

function cost_functional(x, policy_params, t)
    cost(x, policy(x, policy_params))
end

# See https://github.com/SciML/DiffEqSensitivity.jl/issues/302 for context.
dcost_tuple = (
    (out, u, p, t) -> begin
        ū, _, _ = Zygote.gradient(cost_functional, u, p, t)
        out .= ū
    end,
    (out, u, p, t) -> begin
        _, p̄, _ = Zygote.gradient(cost_functional, u, p, t)
        out .= p̄
    end,
)

function gold_standard_gradient(x0, policy_params)
    # Actual/gold standard evaluation. Using high-fidelity Vern9 method with
    # small tolerances. We want to use Float64s for maximum accuracy. Also 1e-14
    # is recommended as the minimum allowable tolerance here: https://docs.sciml.ai/stable/basics/faq/#How-to-get-to-zero-error-1.
    x0_f64 = convert(Array{Float64}, x0)
    policy_params_f64 = convert(Array{Float64}, policy_params)
    fwd_sol = solve(
        ODEProblem(policy_dynamics!, x0, (0, T), policy_params),
        Vern9(),
        u0 = x0_f64,
        p = policy_params_f64,
        abstol = 1e-14,
        reltol = 1e-14,
    )
    # Note that specifying dense = false is essential for getting acceptable
    # performance. save_everystep = false is another small win.
    bwd_sol = solve(
        ODEAdjointProblem(
            fwd_sol,
            InterpolatingAdjoint(),
            cost_functional,
            nothing,
            dcost_tuple,
        ),
        Vern9(),
        dense = false,
        save_everystep = false,
        abstol = 1e-14,
        reltol = 1e-14,
    )
    @assert typeof(fwd_sol.u) == Array{Array{Float64,1},1}
    @assert typeof(bwd_sol.u) == Array{Array{Float64,1},1}

    # Note that the backwards solution includes the gradient on x0, as well as
    # policy_params. The full ODESolution things can't be serialized easily
    # since they `policy_dynamics!` and shit...
    (xT = fwd_sol.u[end], g = bwd_sol.u[end])
end

function eval_interp(x0, policy_params, abstol, reltol)
    fwd_sol = solve(
        ODEProblem(policy_dynamics!, x0, (0, T), policy_params),
        Tsit5(),
        u0 = x0,
        p = policy_params,
        abstol = abstol,
        reltol = reltol,
    )
    bwd_sol = solve(
        ODEAdjointProblem(
            fwd_sol,
            InterpolatingAdjoint(),
            cost_functional,
            nothing,
            dcost_tuple,
        ),
        Tsit5(),
        dense = false,
        save_everystep = false,
        abstol = abstol,
        reltol = reltol,
    )
    @assert typeof(fwd_sol.u) == Array{Array{floatT,1},1}
    @assert typeof(bwd_sol.u) == Array{Array{floatT,1},1}

    # Note that g includes the x0 gradient and the gradient on parameters.
    # We do exactly as many f calls as there are function calls in the forward
    # pass, and in the backward pass we don't need to call f, but instead we
    # call ∇f.
    (
        xT = fwd_sol.u[end],
        g = bwd_sol.u[end],
        nf = fwd_sol.destats.nf,
        n∇f = bwd_sol.destats.nf,
    )
end

## ------------------------------------------------------------------------------------
# using DiffEqFlux, DifferentialEquations, Plots
# using ComponentArrays
# using FSimBase  # for apply_inputs

# u0 = Float32[2.; 0.]
# u0_aug = vcat(u0, zeros(1))
# datasize = 30
# tspan = (0.0f0,1.5f0)

# L(u) = u' * u  # running cost

# function trueODEfunc(dx,x,p,t; a)
#     true_A = [-0.1 2.0; -2.0 -0.1]
#     dx[1:2, :] .= ((x[1:2, :] .^ 3)'true_A)' + a
#     dx[3, :] = L(x[1:2, :])
# end
# t = range(tspan[1],tspan[2],length=datasize)
# # prob = ODEProblem(trueODEfunc,u0,tspan)
# # prob = ODEProblem(trueODEfunc,u0_aug,tspan)
# # ode_data = Array(solve(prob,Tsit5(),saveat=t))

# # dudt2 = Chain(x -> x.^3,
# #              Dense(2,50,tanh),
# #              Dense(50,2))
# # p,re = Flux.destructure(dudt2) # use this p as the initial condition!
# # dudt(u,p,t) = re(p)(u) # need to restructure for backprop!
# # prob = ODEProblem(dudt,u0_aug,tspan)
# controller2 = Chain(u -> u.^3,
#                    Dense(2,50,tanh),
#                    Dense(50,2))
# p,re = Flux.destructure(controller2)
# controller(u, p, t) = re(p)(u)
# prob = ODEProblem(apply_inputs(trueODEfunc; a=(x, p, t) -> controller(x[1:2], p, t)),u0_aug,tspan)

# function predict_n_ode()
#   Array(solve(prob,Tsit5(),u0=u0_aug,p=p,saveat=t))  # it seems that the output should be an Array
# end

# function loss_n_ode()
#     pred = predict_n_ode()
#     # loss = sum(abs2,ode_data .- pred)
#     loss = pred[end][end]  # to make it scalar for gradient calculation
#     loss
# end

# loss_n_ode() # n_ode.p stores the initial parameters of the neural ODE

# cb = function (;doplot=false) # callback function to observe training
#     pred = predict_n_ode()
#     @show pred[end][end]
#     # display(sum(abs2,ode_data .- pred))
#     # # plot current prediction against data
#     # pl = scatter(t,ode_data[1,:],label="data")
#     # scatter!(pl,t,pred[1,:],label="prediction")
#     # display(plot(pl))
#     return false
# end

# # Display the ODE with the initial parameter values.
# cb()

# data = Iterators.repeated((), 1000)
# # Flux.train!(loss_n_ode, Flux.params(u0,p), data, ADAM(0.05), cb = cb)
# Flux.train!(loss_n_ode, Flux.params(p), data, ADAM(0.05), cb = cb)

#
# #-------------------------------------------------------------------------------
# using DiffEqFlux, DifferentialEquations, Plots, Statistics
# tspan = (0.0f0,8.0f0)
# ann = FastChain(FastDense(1,32,tanh), FastDense(32,32,tanh), FastDense(32,1))
# θ = initial_params(ann)
# function dxdt_(dx,x,p,t)
#     x1, x2 = x
#     dx[1] = x[2]
#     dx[2] = ann([t],p)[1]^3
# end
# x0 = [-4f0,0f0]
# ts = Float32.(collect(0.0:0.01:tspan[2]))
# prob = ODEProblem(dxdt_,x0,tspan,θ)
# solve(prob,Vern9(),abstol=1e-10,reltol=1e-10)
# ##

# function predict_adjoint(θ)
#   Array(solve(prob,Vern9(),p=θ,saveat=ts,sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true))))
# end
# function loss_adjoint(θ)
#   x = predict_adjoint(θ)
#   mean(abs2,4.0 .- x[1,:]) + 2mean(abs2,x[2,:]) + mean(abs2,[first(ann([t],θ)) for t in ts])/10
# end
# l = loss_adjoint(θ)
# cb = function (θ,l)
#   println(l)
#   p = plot(solve(remake(prob,p=θ),Tsit5(),saveat=0.01),ylim=(-6,6),lw=3)
#   plot!(p,ts,[first(ann([t],θ)) for t in ts],label="u(t)",lw=3)
#   display(p)
#   return false
# end
# # Display the ODE with the current parameter values.
# cb(θ,l)
# loss1 = loss_adjoint(θ)
# res1 = DiffEqFlux.sciml_train(loss_adjoint, θ, ADAM(0.005), cb = cb,maxiters=100)
# res2 = DiffEqFlux.sciml_train(loss_adjoint, res1.u,
#                               BFGS(initial_stepnorm=0.01), cb = cb,maxiters=100,
#                               allow_f_increases = false)

##-------------------------------------------------------------------------------
using DiffEqFlux, DiffEqSensitivity, DifferentialEquations, Plots, Statistics
tspan = (0.0f0,8.0f0)
ann = FastChain(FastDense(1,32,tanh), FastDense(32,32,tanh), FastDense(32,1))
θ = initial_params(ann)
function dxdt_(x,p,t)
    x1, x2 = x
    return dx = [x2; ann([t],p)[1]^3]
end
x0 = [-4f0,0f0]
ts = Float32.(collect(0.0:0.01:tspan[2]))
prob = ODEProblem(dxdt_,x0,tspan,θ)
# solve(prob,Vern9(),abstol=1e-10,reltol=1e-10)
# function predict_adjoint(θ)
#   Array(solve(prob,Vern9(),p=θ,saveat=0.01f0,sensealg=InterpolatingAdjoint()))
# end
function loss_adjoint_test(θ)
#   x = predict_adjoint(θ)
  x = Array(solve(prob,Tsit5(),p=θ,saveat=0.01f0,sensealg=InterpolatingAdjoint( autojacvec = ZygoteVJP() )))
  mean(abs2,4.0f0 .- x[1,:]) + 2mean(abs2,x[2,:]) + mean(abs2,[first(ann([t],θ)) for t in ts])/10.0f0
end
l = loss_adjoint_test(θ)
cb = function (θ,l)
  println(l)
  p = plot(solve(remake(prob,p=θ),Tsit5(),saveat=0.01),ylim=(-6,6),lw=3)
  plot!(p,ts,[first(ann([t],θ)) for t in ts],label="u(t)",lw=3)
  display(p)
  return false
end
# Display the ODE with the current parameter values.

cb(θ,l)
loss1 = loss_adjoint_test(θ)
res1 = DiffEqFlux.sciml_train(loss_adjoint_test, θ, ADAM(0.005); cb = cb, maxiters=100)
res2 = DiffEqFlux.sciml_train(loss_adjoint_test, res1.u,
                              BFGS(initial_stepnorm=0.01), cb = cb,maxiters=100,
                              allow_f_increases = false)

## -------------------------------------------------------------------------
using OrdinaryDiffEq, DiffEqFlux
pa = [1.0]
u0 = [3.0]
θ = [u0;pa]

function model1(θ,ensemble)
  prob = ODEProblem((u, p, t) -> 1.01u .* p, [θ[1]], (0.0, 1.0), [θ[2]])

  function prob_func(prob, i, repeat)
    remake(prob, u0 = 0.5 .+ i/100 .* prob.u0)
  end

  ensemble_prob = EnsembleProblem(prob, prob_func = prob_func)
  sim = solve(ensemble_prob, Tsit5(), ensemble, saveat = 0.1, trajectories = 100, sensealg = InterpolatingAdjoint(autojacvec = ZygoteVJP()))
end

# loss function
loss_serial(θ)   = sum(abs2,1.0.-Array(model1(θ,EnsembleSerial())))
loss_threaded(θ) = sum(abs2,1.0.-Array(model1(θ,EnsembleThreads())))

prob_base = ODEProblem((u, p, t) -> 1.01u .* p, [θ[1]], (0.0, 1.0), [θ[2]])
function prob_gen(prob, i, repeat)
  remake(prob, u0 = 0.5 .+ i/100 .* prob.u0)
end
function loss_test(θ)
    prob = remake(prob_base, u0 = [θ[1]], p = [θ[2]])
    fwd_ensemble_sol = Array( solve(EnsembleProblem(prob, prob_func = prob_gen), Tsit5(), EnsembleThreads(), saveat = 0.1, trajectories = 100, sensealg = InterpolatingAdjoint(autojacvec = ZygoteVJP())))
    loss_val = sum(abs2, 1.0 .- fwd_ensemble_sol)
    return loss_val
end

cb = function (θ,l) # callback function to observe training
  @show l
  false
end

opt = ADAM(0.1)
l1 = loss_serial(θ)
# res_serial = DiffEqFlux.sciml_train(loss_serial, θ, opt; cb = cb, maxiters=100)
# res_threads = DiffEqFlux.sciml_train(loss_threaded, θ, opt; cb = cb, maxiters=100)
res_test = DiffEqFlux.sciml_train(loss_test, θ, opt; cb = cb, maxiters = 100)

## ----------------
using DiffEqGPU, OrdinaryDiffEq
function lorenz(du,u,p,t)
    du[1] = p[1]*(u[2]-u[1])
    du[2] = u[1]*(p[2]-u[3]) - u[2]
    du[3] = u[1]*u[2] - p[3]*u[3]
end

u0 = Float32[1.0;0.0;0.0]
tspan = (0.0f0,100.0f0)
p = [10.0f0,28.0f0,8/3f0]
prob = ODEProblem(lorenz,u0,tspan,p)
prob_func = (prob,i,repeat) -> remake(prob,p=rand(Float32,3).*p)
monteprob = EnsembleProblem(prob, prob_func = prob_func, safetycopy=false)
@time sol = solve(monteprob,Tsit5(),EnsembleGPUArray(),trajectories=10000,saveat=1.0f0)

## -------------------
using OrdinaryDiffEq, Flux, Optim, DiffEqFlux, DiffEqSensitivity

model_gpu = Chain(Dense(2, 50, tanh), Dense(50, 2)) |> gpu
p, re = Flux.destructure(model_gpu)
dudt!(u, p, t) = re(p)(u)

# Simulation interval and intermediary points
tspan = (0.0, 10.0)
tsteps = 0.0:0.1:10.0

u0 = Float32[2.0; 0.0] |> gpu
prob_gpu = ODEProblem(dudt!, u0, tspan, p)

# Runs on a GPU
sol_gpu = solve(prob_gpu, Tsit5(), saveat = tsteps)


text:
abstract type AbstractOptimizerState end

"""
    value(state)

返回当前状态的极小值
"""
value(state::AbstractOptimizerState) = error("`value` is not implemented for $(state).")

"""
    minimizer(state)

返回当前状态的决策变量
"""
minimizer(state::AbstractOptimizerState) = error("`minimizer` is not implemented for $(state).")

"""
    terminate(state)

判断是否提前终止
"""
terminate(state::AbstractOptimizerState) = false

ext: jl


text:
using Dates: DateTime
using Registrator.WebUI: isauthorized, AuthFailure, AuthSuccess, User
using GitForge: GitForge, GitHub, GitLab
using HTTP: stacktrace

using Mocking

Mocking.activate()

function patch_gitforge(body::Function; is_collaborator=false, is_member=false)
    patches = [
        @patch GitForge.is_collaborator(args...) =
            GitForge.Result{Bool}(is_collaborator, nothing, nothing, stacktrace())
        @patch GitForge.is_member(args...) =
            GitForge.Result{Bool}(is_member, nothing, nothing, stacktrace())
    ]

    apply(patches) do
        return body()
    end
end

@testset "gitutils" begin

@testset "isauthorized" begin
    @test isauthorized("username", "reponame") == AuthFailure("Unkown user type or repo type")

    @testset "GitHub" begin

        user = GitHub.User(login="user123")
        org = GitHub.User(login="JuliaLang")
        private_repo = GitHub.Repo(name="Example.jl", private=true, owner=user)
        public_repo_of_user = GitHub.Repo(name="Example.jl", private=false, owner=user, organization=nothing)
        public_repo_of_org = GitHub.Repo(name="Example.jl", private=false, owner=org, organization=org)
        u = User(user, GitHub.GitHubAPI())

        @testset "private repo" begin
            # Assuming CONFIG["allow_private"] is false
            @test isauthorized(u, private_repo) == AuthFailure("Repo Example.jl is private")
        end

        @testset "public repo of user" begin
            # authorized if user is a collaborator on the repo
            patch_gitforge(is_collaborator=true) do
                @test isauthorized(u, public_repo_of_user) == AuthSuccess()
            end
            patch_gitforge(is_collaborator=false) do
                @test isauthorized(u, public_repo_of_user) == AuthFailure("User user123 is not a collaborator on repo Example.jl")
            end
        end

        @testset "public repo of org" begin
            # authorized if user is either a collaborator on the repo or member of the org
            patch_gitforge(is_collaborator=true, is_member=true) do
                @test isauthorized(u, public_repo_of_org) == AuthSuccess()
            end
            patch_gitforge(is_collaborator=true, is_member=false) do
                @test isauthorized(u, public_repo_of_org) == AuthSuccess()
            end
            patch_gitforge(is_collaborator=false, is_member=true) do
                @test isauthorized(u, public_repo_of_org) == AuthSuccess()
            end
            patch_gitforge(is_collaborator=false, is_member=false) do
                @test isauthorized(u, public_repo_of_org) == AuthFailure("User user123 is not a member of the org JuliaLang and not a collaborator on repo Example.jl")
            end
        end
    end

    @testset "GitLab" begin

        user = GitLab.User(name="user123", username="user123", id=111)
        org = GitLab.User(name="org123", username="org123", id=222)
        private_project = GitLab.Project(name="Example.jl", visibility="private", owner=user)
        public_project_of_user = GitLab.Project(name="Example.jl", visibility="public", owner=user, namespace=GitLab.Namespace(kind="user"))
        public_project_of_group = GitLab.Project(name="Example.jl", visibility="public", owner=org, namespace=GitLab.Namespace(kind="group", full_path="org123/subgroup/Example.jl"))
        u = User(user, GitLab.GitLabAPI())

        @testset "private project" begin
            # Assuming CONFIG["allow_private"] is false
            @test isauthorized(u, private_project) == AuthFailure("Project Example.jl is private")
        end

        @testset "public project of user" begin
            # authorized if user is a collaborator on the project
            patch_gitforge(is_collaborator=true) do
                @test isauthorized(u, public_project_of_user) == AuthSuccess()
            end
            patch_gitforge(is_collaborator=false) do
                @test isauthorized(u, public_project_of_user) == AuthFailure("User user123 is not a member of project Example.jl")
            end
        end

        @testset "public project of group" begin
            # authorized if user is a collaborator on the project or member of the group/subgroups
            patch_gitforge(is_collaborator=true, is_member=true) do
                @test isauthorized(u, public_project_of_group) == AuthSuccess()
            end
            patch_gitforge(is_collaborator=false, is_member=true) do
                @test isauthorized(u, public_project_of_group) == AuthSuccess()
            end
            patch_gitforge(is_collaborator=true, is_member=false) do
                @test isauthorized(u, public_project_of_group) == AuthSuccess()
            end
            patch_gitforge(is_collaborator=false, is_member=false) do
                @test isauthorized(u, public_project_of_group) == AuthFailure("Project Example.jl belongs to the group org123/subgroup/Example.jl, and user user123 is not a member of that group or its parent group(s)")
            end
        end
    end
end

end


text:
function union_constraint_types(MCT)
    UMCT = SOS.union_constraint_indices_types(MCT)
    return Union{UMCT, Vector{UMCT}}
    # `Vector{UMCT}` is for sparse SOS
    # As basis can have different sizes, we need `Vector{Union{...}}` instead of
    # `Union{Vector, ...}`.
end
function union_set_types(MCT)
    return Union{typeof(SOS.matrix_cone(MCT, 0)),
                 typeof(SOS.matrix_cone(MCT, 1)),
                 typeof(SOS.matrix_cone(MCT, 2)),
                 typeof(SOS.matrix_cone(MCT, 3))}
end
function constrained_variable_types(MCT)
    return [
        (typeof(SOS.matrix_cone(MCT, 0)),),
        (typeof(SOS.matrix_cone(MCT, 1)),),
        (typeof(SOS.matrix_cone(MCT, 2)),),
        (typeof(SOS.matrix_cone(MCT, 3)),)
   ]
end


text:
#===============================================================================

                                Cloning

                This is quite inefficient, but has a nice syntax

===============================================================================#

"""
    clone(P::T, θ::AbstractDict) where T <: DiffusionProcess

Simplified cloning of diffusion law `P`. Substitute relevant parameters with
new values. `θ` must be a dict corresponding to parameters returned after a call
to `var_parameters`.
"""
function clone(P::T, θ::AbstractDict) where T <: DiffusionProcess
    cp = const_parameters(P)
    remove_curly(T)(;cp..., θ...)
end

#===============================================================================

                        in-place parameter setting

            Efficient way of re-parameterizing the diffusion law

===============================================================================#
"""
    set_parameters!(P::DiffusionProcess, θ, entries)

Set parameters of a diffusion law `P` in-place. `entries` should be a collection
of pairs Pair{Int64,Symbol} that list the relevant entries in `θ` for
reparameterization, together with the corresponding parameter names.
"""
function set_parameters!(P::DiffusionProcess, θ, entries)
    for (i,s) in entries
        setfield!(P, s, θ[i])
    end
end

function set_parameters!(P::DiffusionProcess, P°::DiffusionProcess, entries)
    for e in entries
        setfield!(P°, e, getfield(P, e))
    end
end

function set_parameters!(P::T, θ::Dict) where T <: DiffusionProcess
    for k in keys(θ)
        hasfield(T, k) && setfield!(P, k, θ[k])
    end
end

function same_entries(P::DiffusionProcess, P°::DiffusionProcess, entries)
    for e in entries
        getfield(P, e) == getfield(P°, e) || return false
    end
    true
end

function has_any(P::T, entries) where T <: DiffusionProcess
    for (i,s) in entries
        hasfield(T, s) && return true
    end
    false
end


text:
# abstract types
abstract type PointProcess end;

"""
    DeterminantalPointProcess(L::AbstractMatrix; parallelize=false)
    DeterminantalPointProcess(L::Eigen; parallelize=false)

Given a symmetric matrix `L`, creates a `DeterminantalPointProcess` (DPP). One can also pass the eigen object directly

---

    DeterminantalPointProcess(kernel::Kernel, X::AbstractVector; parallelize=false)
    DeterminantalPointProcess(kernel::Kernel, X::AbstractMatrix; obsdim=1; parallelize=false)

Similar to the basic constructor, will first build the kernel matrix with `kernel` on observations `X`.
If your input is a `AbstractMatrix` you can pass the `obsdim` argument to indicate if rows or columns represent the samples (see docs of [KernelFunctions](https://juliagaussianprocesses.github.io/KernelFunctions.jl/stable/))
"""
struct DeterminantalPointProcess{T,parallelize,TL<:AbstractMatrix{T},Tfact} <: PointProcess
    L::TL
    Lfact::Tfact
    size::Int
end

function DeterminantalPointProcess(L::AbstractMatrix; parallelize=false)
    issymmetric(L) || error("Given matrix is not symmetric")
    Lfact = eigen(L)
    return DeterminantalPointProcess{eltype(L),parallelize,typeof(L),typeof(Lfact)}(L, Lfact, length(Lfact.values))
end

function DeterminantalPointProcess(Lfact::Eigen; parallelize=false)
    L = Symmetric((Lfact.vectors .* Lfact.values') * Lfact.vectors')
    return DeterminantalPointProcess{eltype(L),parallelize,typeof(L),typeof(Lfact)}(L, Lfact, length(Lfact.values))
end

"""
    kDeterminantalPointProcess(k::Int, dpp::DeterminantalPointProcess)

Create a k-DPP where the size of the subsets is fixed to k.
You can also create a k-DPP by calling `dpp(k)`
"""
struct kDeterminantalPointProcess{T,parallelize,Tdpp<:DeterminantalPointProcess{T,parallelize}} <: PointProcess
    k::Int
    dpp::Tdpp
end

(dpp::DeterminantalPointProcess{T,p})(k::Int) where {T,p} = kDPP{T,p,typeof(dpp)}(k, dpp)

# aliases
const DPP = DeterminantalPointProcess
const kDPP = kDeterminantalPointProcess
# const KDPP = KroneckerDeterminantalPointProcess
const MCMCState = Tuple{BitArray{1},Matrix{Float64}}


text:
"""
    beamforming(csm,v)

Calculate frequency-domain beamforming using cross-spectral matrix `csm` and steering vector `v`.
`csm` must be a square (Hermitian) matrix optionally with a third (frequency) dimension.
First dimension of `v` and `csm` must be equal.
"""
function beamforming(csm::A,v::B) where {T1<:AbstractFloat,T2<:AbstractFloat, N,  A<:AbstractArray{Complex{T1},N},B<:AbstractArray{Complex{T2},N}}
    #@assert size(csm,1) == size(csm,2) "csm must be square with dimensions M x M (x Nf)!"
    #@assert size(v,1) == size(csm,1) "First dimension of v and csm must be equal!"
    b = Array{T2, 2}(undef, size(v,2), size(csm,3))
    for j in axes(csm,3)
        csmd = @view csm[:,:,j]
        for i in axes(v,2)
            vd = @view v[:,i,j]
            b[i,j] = real(vd'*csmd*vd)
        end
    end
    return b
end

"""
    beamforming(Environment)

Calculate frequency-domain beamforming using cross-spectral matrix `csm` and steering vector `v`
stored in Environment struct
"""
function beamforming(E::Environment)
    @unpack CSM_s,w,steeringvec,M,N,Nf,fn = E
    b = Array{Float64}(undef, N, Nf)
    @views @inbounds for j in 1:Nf
        st = E.steeringvec.arr[:,:,j]
        @. st = complex(w[:,j].*real(st), imag(st))
        bf_col!(b[:,j],st,selectdim(CSM_s.arr, 3, j))
    end
    return FreqArray(b,fn)
end

function bf_col!(b::A,st::B,csm::C) where {T1<:AbstractFloat,T2<:AbstractFloat, N, A<:AbstractArray{T1,1}, B<:AbstractArray{Complex{T1},N},C<:AbstractArray{Complex{T2},N}}
    @views @inbounds for i in 1:length(b)
        b[i] = real(mydot(st[:,i]',csm,st[:,i]))
    end
end

function mydot(xt, A, x)
    s = zero(promote_type(eltype(x), eltype(A)))
    @inbounds @simd for m in 1:size(A,1)
        for n in 1:size(A,2)
            s += xt[m] * A[m,n] * x[n]
        end
    end
    s
end


text:
##k-mean++

#Square Euclidean Function to use (internal sql2() or SqEuclidean() from Distances.jl)
SqEuDist = SqEuclidean
# SqEuDist = sql2

#Return distance to closest centre for each data point (helper function for k-mean++)
"""
    mindist2cntrs(centres, data)

Computes the square of the L2 distance to closest centre for each data point.
* `data` must be a matrix of `Float` with each column representing a point
* `centres` must be a matrix of `Float` with each column representing a centre

!!! note
    `mindist2cntrs(centres, data)` does not check the whether the dimensions of the input match.
    Users must ensure that `size(data,1) = size(centres,1)`.
"""
function mindist2cntrs(cntrs::AbstractMatrix{T}, data::AbstractMatrix{T}) where T<:AbstractFloat
    mdist = Vector{T}(undef,size(data,2))
    temp = Vector{T}(undef,size(cntrs,2))
    @inbounds @views for i ∈ eachindex(mdist)
        colwise!(temp,SqEuDist(), data[:,i], cntrs)
        mdist[i] = minimum(temp)
    end
    return mdist
end

#k-means++ algorithm for initialising centres
"""
    kmeanspp(data, number_of_centres)

Indentfies initial cluster centres for k-means clustering using the k-means++ alogrithm for given `data`.
* `data` must be a matrix of `Float` with each column vector representing a point
* Function will return a vector of indices corresponding to columns of `data` identified as cluster centres
"""
function kmeanspp(data::AbstractMatrix{T}, ncntrs::Int) where T<:AbstractFloat
    cin = [rand(1:size(data,2))]
    din = collect(1:size(data,2))
    deleteat!(din,cin)
    while(length(cin)<ncntrs)
        mindist = mindist2cntrs(data[:,cin], data[:,din])
        cdist = pdist2cdist(mindist)
        ncin = findclosest(rand(),cdist)
        push!(cin, din[ncin])
        deleteat!(din,ncin)
    end
    return cin
end

## K-means

#Identify closest cluster centre
"""
    closest_centre!(closest_centre, centres, data)

Identifies the centre closest to each data point based on L2 distance.
* `data` must be a matrix of `Float` with each column representing a point
* `centres` must be a matrix of `Float` with each column representing a centre
* `closest_centre` must be a vector of `Int` with each element holding the column index of `centres` that corresponds to the closest centre

!!! note
    `closest_centre!(closest_centre, centres, data)` does not check the whether the dimensions of the input match.
    Users must ensure that `length(closest_centre) = size(data,2)` and `size(data,1) = size(centres,1)`.
"""
function closest_centre!(ccntr::AbstractVector{Int}, cntrs::AbstractMatrix{T}, data::AbstractMatrix{T}) where T<:AbstractFloat
    temp = Vector{T}(undef,size(cntrs,2))
    @inbounds @views for i ∈ eachindex(ccntr)
        colwise!(temp,SqEuDist(), data[:,i], cntrs)
        ccntr[i] = argmin(temp)
    end
end

"""
    closest_centre(centres, data)

Identifies the centre closest to each data point based on L2 distance.
* `data` must be a matrix of `Float` with each column representing a point
* `centres` must be a matrix of `Float` with each column representing a centre
* Output is a vector of `Int` with each element holding the column index of `centres` that corresponds to the closest centre

!!! note
    `closest_centre(centres, data)` does not check the whether the dimensions of the input match.
    Users must ensure that `size(data,1) = size(centres,1)`.
"""
function closest_centre(cntrs::AbstractMatrix{T}, data::AbstractMatrix{T}) where T<:AbstractFloat
    ccntr = Vector{Int}(undef,size(data,2))
    closest_centre!(ccntr, cntrs, data)
    return ccntr
end

#Maximum L2 distance by which a centre has shifted
function max_centre_shift(new_cntrs::AbstractMatrix{T}, old_cntrs::AbstractMatrix{T}) where T<:AbstractFloat
    s = zero(T)
    @views for i ∈ axes(old_cntrs,2)
        temp = evaluate(SqEuDist(),new_cntrs[:,i],old_cntrs[:,i])
        if(temp>s)
            s = temp
        end
    end
    return sqrt(s)
end

function max_centre_shift2(new_cntrs::AbstractMatrix{T}, old_cntrs::AbstractMatrix{T}) where T<:AbstractFloat
    temp = colwise(SqEuDist(),new_cntrs,old_cntrs)
    return sqrt(maximum(temp))
end

#k-means initialisation
abstract type KMinit end
struct kmpp <: KMinit end
struct kmrand <: KMinit end

"""
    kmeans_init([initialise_method], data, number_of_centres)

Pick a set of points from `data` to serve as the initial cluster centres for k-means clustering. Method for initialisation can random or k-means++.
* `data` must be a matrix of `Float` with each column representing a point
* Use `kmpp()` as `initialise_method` to use k-mean++ as initialisation method [DEFAULT if `initialise_method` is not specified]
* Use `kmrand()` as `intialise_method` to randomly select points as initial centres. This can be much faster than k-means++ for large datasets
"""
kmeans_init(::kmpp, data::AbstractMatrix{<:AbstractFloat}, ncntrs::Int) = data[:,kmeanspp(data,ncntrs)]

kmeans_init(::kmrand, data::AbstractMatrix{<:AbstractFloat}, ncntrs::Int) = data[:,rand(1:size(data,2),ncntrs)]

kmeans_init(data::AbstractMatrix{<:AbstractFloat}, ncntrs::Int) = kmeans_init(kmpp(), data, ncntrs)

#One iteration of updating k-means cluster centres
function kmeans_update!(ccntr::AbstractVector{Int}, cntrs::AbstractMatrix{T}, data::AbstractMatrix{T}) where T<:AbstractFloat
    closest_centre!(ccntr,cntrs,data)
    for i ∈ axes(cntrs,2)
        if(sum(ccntr.==i)>0)
            cntrs[:,i] = running_mean(data[:,ccntr.==i])
        end
    end
end

"""
    kmeans!(centres, data; [update_threshold = 0.0])

Updates (inplace) cluster centres `centres` using the k-means clustering algorithm based on `data`.
Cluster centres are iteratively updated until the largest shift in cluster centres,
measured in terms of square of L2 distance, falls below `update_threshold`.
* `centres` must be a matrix of `Float` with each column representing a cluster centre
* `data` must be a matrix of `Float` with each column representing a point
* Default value of optional named parameter is `update_threshold = 0.0` and corresponds to the setting that algorithm will converge when cluster centres do not shift
"""
function kmeans!(cntrs::AbstractMatrix{T}, data::AbstractMatrix{T}; update_threshold = 0.0) where T<:AbstractFloat
    ndims,ncntrs = size(cntrs)
    old_cntrs = Matrix{T}(undef,ndims,ncntrs)
    ccntr = Vector{Int}(undef,size(data,2))
    max_shift::T = Inf
    while(max_shift > update_threshold)
        copy!(old_cntrs,cntrs)
        kmeans_update!(ccntr, cntrs, data)
        max_shift = max_centre_shift(cntrs, old_cntrs)
    end
end

"""
    kmeans(data, number_of_centres; [update_threshold = 0.0])
    kmeans(data, number_of_centres, initialise_method; [update_threshold = 0.0])

Estimate cluster centres from `data` using the k-means clustering algorithm.
Cluster centres are iteratively updated until the largest shift in cluster centres,
measured in terms of square of L2 distance, falls below `update_threshold`.
* `data` must be a matrix of `Float` with each column representing a point
* Default value of optional named parameter is `update_threshold = 0.0` and corresponds to the setting that algorithm will converge when cluster centres do not shift
* Use `kmpp()` as `initialise_method` to use k-mean++ as initialisation method [DEFAULT if `initialise_method` is not specified]
* Use `kmrand()` as `intialise_method` to randomly select points as initial centres. This can be faster than k-means++ for large datasets
"""
function kmeans(data::AbstractMatrix{T}, ncntrs::Int, init_method::KMinit; update_threshold = 0.0) where T<:AbstractFloat
    cntrs = kmeans_init(init_method, data, ncntrs)
    kmeans!(cntrs, data; update_threshold)
    return cntrs
end

kmeans(data::AbstractMatrix{T}, ncntrs::Int; update_threshold = 0.0) where T<:AbstractFloat = kmeans(data, ncntrs, kmpp(); update_threshold)


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

function _enforce_transmission(
    model::JuMP.Model,
    violations::Vector{_Violation},
)::Nothing
    for v in violations
        _enforce_transmission(
            model = model,
            violation = v,
            isf = model[:isf],
            lodf = model[:lodf],
        )
    end
    return
end

function _enforce_transmission(;
    model::JuMP.Model,
    violation::_Violation,
    isf::Matrix{Float64},
    lodf::Matrix{Float64},
)::Nothing
    instance = model[:instance]
    limit::Float64 = 0.0
    overflow = model[:overflow]
    net_injection = model[:net_injection]

    if violation.outage_line === nothing
        limit = violation.monitored_line.normal_flow_limit[violation.time]
        @info @sprintf(
            "    %8.3f MW overflow in %-5s time %3d (pre-contingency)",
            violation.amount,
            violation.monitored_line.name,
            violation.time,
        )
    else
        limit = violation.monitored_line.emergency_flow_limit[violation.time]
        @info @sprintf(
            "    %8.3f MW overflow in %-5s time %3d (outage: line %s)",
            violation.amount,
            violation.monitored_line.name,
            violation.time,
            violation.outage_line.name,
        )
    end

    fm = violation.monitored_line.name
    t = violation.time
    flow = @variable(model, base_name = "flow[$fm,$t]")

    v = overflow[violation.monitored_line.name, violation.time]
    @constraint(model, flow <= limit + v)
    @constraint(model, -flow <= limit + v)

    if violation.outage_line === nothing
        @constraint(
            model,
            flow == sum(
                net_injection[b.name, violation.time] *
                isf[violation.monitored_line.offset, b.offset] for
                b in instance.buses if b.offset > 0
            )
        )
    else
        @constraint(
            model,
            flow == sum(
                net_injection[b.name, violation.time] * (
                    isf[violation.monitored_line.offset, b.offset] + (
                        lodf[
                            violation.monitored_line.offset,
                            violation.outage_line.offset,
                        ] * isf[violation.outage_line.offset, b.offset]
                    )
                ) for b in instance.buses if b.offset > 0
            )
        )
    end
    return nothing
end


text:
#ДАНО: Робот находится в произвольной клетке ограниченного прямоугольного поля с внутренними перегородками.
#РЕЗУЛЬТАТ: Робот — в исходном положении в центре прямого креста из маркеров, расставленных вплоть до внешней рамки.
include("library.jl")
function krest_with_obstacles(r, dx::Integer = 0, dy::Integer = 0)

    R = CrossRobot(r, Coords(0, 0), true, dx, dy)
    back_path = BackPath(R)
    marks!(R)
    BackPath(r)
    back!(R, back_path)
end

max_stars_repo_name: listova-nastya/mirea-progs

max_issues_repo_name: listova-nastya/mirea-progs

max_forks_repo_name: listova-nastya/mirea-progs


text:
using MultivariateOrthogonalPolynomials, Test

include("test_modalinterlace.jl")
include("test_disk.jl")
include("test_rectdisk.jl")
include("test_triangle.jl")
# include("test_dirichlettriangle.jl")


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export TropicalWeight

"""
`TropicalWeight(x)`

| Set                                 | ``\\oplus``|   ``\\otimes`` | ``\\bar{0}`` | ``\\bar{1}`` |
|:-----------------------------------:|:----------:|:--------------:|:------------:|:------------:|
|``\\mathbb{R}\\cup\\{\\pm\\infty\\}``|  ``\\min`` |     ``+``      |``\\infty``   |   ``0``      |
"""
struct TropicalWeight{T <: AbstractFloat} <: Semiring
  x::T
end

zero(::Type{TropicalWeight{T}}) where T = TropicalWeight{T}(T(Inf))
one(::Type{TropicalWeight{T}}) where T = TropicalWeight{T}(zero(T))

*(a::TropicalWeight{T}, b::TropicalWeight{T}) where {T <: AbstractFloat} = TropicalWeight{T}(a.x + b.x)
+(a::TropicalWeight{T}, b::TropicalWeight{T}) where {T <: AbstractFloat} = TropicalWeight{T}( min(a.x,b.x) )
/(a::TropicalWeight{T}, b::TropicalWeight{T}) where {T <: AbstractFloat} = TropicalWeight{T}( a.x-b.x )

# parsing
parse(::Type{S},str) where {T, S <: TropicalWeight{T}} = S(parse(T,str))

#properties
isidempotent(::Type{W}) where {W <: TropicalWeight} = true
iscommulative(::Type{W}) where {W <: TropicalWeight} = true
isleft(::Type{W}) where {W <: TropicalWeight}= true
isright(::Type{W}) where {W <: TropicalWeight}= true
isweaklydivisible(::Type{W}) where {W <: TropicalWeight}= true
iscomplete(::Type{W}) where {W <: TropicalWeight}= true
ispath(::Type{W}) where {W <: TropicalWeight}= true


text:
"""
    IntDict{T}

A dictionary to convert key::T into integer id.

## 👉 Example
```julia
dict = IntDict{String}()
push!(dict, "abc") == 1
push!(dict, "def") == 2
push!(dict, "abc") == 1

dict["abc"] == 1
getkey(dict, id1) == "abc"
```
"""
struct IntDict{T}
    key2id::Dict{T,Int}
    id2key::Vector{T}
    id2count::Vector{Int}
    default::Int
end
IntDict{T}(; default=0) where T = IntDict(Dict{T,Int}(), T[], default)

Base.count(dict::IntDict, id::Int) = dict.id2count[id]
Base.getkey(dict::IntDict, id::Int) = dict.id2key[id]
Base.getindex(dict::IntDict, key) = dict.key2id[key]
Base.get(dict::IntDict, key) = get(dict.key2id, key, dict.default)
Base.length(d::IntDict) = length(d.key2id)

function Base.push!(dict::IntDict, key)
    if haskey(dict.key2id, key)
        id = dict.key2id[key]
        dict.id2count[id] += 1
    else
        id = length(dict.id2key) + 1
        dict.key2id[key] = id
        push!(dict.id2key, key)
        push!(dict.id2count, 1)
    end
    id
end
Base.append!(dict::IntDict, keys::Vector) = map(k -> push!(dict,k), keys)

function save(dict::IntDict)
    throw("Not implemented yet.")
end

function load(::Type{IntDict{T}}, path::String) where T
    throw("Not implemented yet.")
end


text:
using Printf

function variable_examples()
    println("--------------------------------------")
    println("variable_examples")
    println("Variables types can be dynamically assigned")
    s = 0
    println(s)
    s = true
    println(s)
    s = "Dog"
    println(s)

    println("We can also specify variable types")
    x::Int8 = 10
    println(x)
    y::Int16 = -10
    println(y)

    # If we try to assign a different type to a typed variable, we get an exception
    #y = 3.14
    #println(y)

    # Int8 : -128 - 127
    # Int16 : −32,768 - 32,767
    # Int32 : −2,147,483,648 - 2,147,483,647
    # Int64 : −2^63 - 2^63 - 1
    # Float32, Float64, UInt8, UInt16, etc.
    # UInt8, UInt16, etc.

    println("We can also use booleans")
    somebool = true
    println(somebool)

    println("Chars have to be surrounded by single quotes")
    c = 'c'
    println(c)

    println("Cast from Int to Char")
    c2 = Char(120)
    println(c2)

    println("Cast from Float to Int")
    i1 = UInt8(trunc(3.14))
    println(i1)

    println("Cast string to Float")
    f1 = parse(Float64,"1")
    println(f1)

    println("Cast string to Int")
    i2 = parse(Int8,"1")
    println(i2)

    println("We can also create BigInts")
    someBigInt::BigInt = 1234 * 1234 * 1234 * 1234
    println(someBigInt)

    println("And we can create BigFloats")
    someBigFloat::BigFloat = 1234 / 1234 / 1234 / 1234
    println(someBigFloat)

end

max_stars_repo_name: James-P-D/JuliaDump

max_issues_repo_name: James-P-D/JuliaDump

max_forks_repo_name: James-P-D/JuliaDump


text:
using .InteractiveViz
using Statistics
using DSP.Periodograms

export iplot, iplot!, ispecgram

"""
$(SIGNATURES)
Plots interactive timeseries of the signal.
"""
function InteractiveViz.iplot(s::SampledSignal, args...; kwargs...)
  s1 = samples(s)
  if isanalytic(s1)
    @warn "Plotting only real part of complex signal"
    s1 = real.(s1)
  end
  t = domain(s)
  if maximum(t) <= 1.0
    t *= 1000.0
    xlabel = "Time (ms)"
  else
    xlabel = "Time (s)"
  end
  iplot(t, s1, args...; xlabel=xlabel, kwargs...)
end

"""
$(SIGNATURES)
Plots interactive timeseries of the signal over a previous plot.
"""
function InteractiveViz.iplot!(s::SampledSignal; kwargs...)
  s1 = samples(s)
  if isanalytic(s1)
    @warn "Plotting only real part of complex signal"
    s1 = real.(s1)
  end
  t = domain(s)
  maximum(t) <= 1.0 && (t *= 1000.0)
  iplot!(t, s1; kwargs...)
end

"""
$(SIGNATURES)
Plots interactive spectrogram of the signal.
"""
function ispecgram(s; fs=framerate(s), nfft=min(div(length(s),8),256), noverlap=div(nfft,2), window=nothing, pooling=mean, kwargs...)
  @assert Base.size(s,2) == 1 "ispecgram only works with vectors"
  s1 = samples(s)[:,1]
  if isanalytic(s1)
    @warn "Using only real part of complex signal"
    s1 = real.(s1)
  end
  p = spectrogram(s1, nfft, noverlap; fs=inHz(fs), window=window)
  t = time(p)
  f = freq(p)
  tunit = "s"
  funit = "Hz"
  if maximum(f) >= 10000
    f /= 1000.0
    funit = "kHz"
  end
  if maximum(t) <= 1.0
    t *= 1000.0
    tunit = "ms"
  end
  xlabel = "Time ("*tunit*")"
  ylabel = "Frequency ("*funit*")"
  z = pow2db.(power(p)')
  cmax = ceil(Int, maximum(z)/5)*5
  cmin = floor(Int, max(cmax-50, minimum(z))/5)*5
  iheatmap(z, t[1], t[end], f[1], f[end]; clim=(cmin, cmax), xlabel=xlabel, ylabel=ylabel, pooling=pooling, kwargs...)
end


text:
# Tests that the weighting scheme works for multiclass to balance skewed data

# Get margin for each sample in votes, according to its true label in labels
function margin{T1 <: Real, T2 <: Real}( votes::Array{T1,2}, labels::Vector{T2} )
    const numClasses = size(votes,1)
    const numSamples = size(votes,2)

    votes = copy(votes)

    mg = zeros(eltype(votes), numSamples)
    for i=1:numSamples
        local v = votes[Int(labels[i]) + 1, i]  # copy score of true label

        # compute margin
        votes[Int(labels[i]) + 1, i] = -Inf
        mg[i] = v - maximum(votes[:,i])
    end

    return mg
end

# Create a skewed 3-class problem
N1 = 100
N2 = 100000
N3 = 1000000
N = N1 + N2 + N3

numClasses = 3
numFeats = 5

X = rand((N, numFeats))
y = vcat(zeros(N1), ones(N2), 2*ones(N3))

# one feature is exactly the label
X[:,1] = y

estimator = LightGBM.LGBMMulticlass(num_iterations = 10,
                                    learning_rate = .5,
                                    feature_fraction = 0.5,
                                    bagging_fraction = 1.0,
                                    bagging_freq = 1,
                                    num_leaves = 5,
                                    metric = ["multi_logloss"],
                                    num_class = numClasses,
                                    min_data_in_leaf=1,
                                    min_sum_hessian_in_leaf=1);

# Compute required weights
classWeights = [N / (numClasses * sum(y .== n)) for n in 0:(numClasses - 1)]
weights = [classWeights[Int(n+1)] for n in y]

LightGBM.fit(estimator, X, y; weights=weights);

pred = LightGBM.predict(estimator, X; predict_type=1)
pred = reshape(pred, ( numClasses, size(X,1)))

# Average margins for each prediction should be similar if weighting works
margins = margin(pred, y)
margins = [mean(margins[y .== label]) for label in 0:(numClasses-1)]

meanMargin = mean(margins)
marginDeviations = abs.(margins - meanMargin)
maxMarginDeviation = maximum(marginDeviations)

@test maxMarginDeviation / meanMargin < 0.1


text:
######################################
# Example of rksm function
# Based on example_rksm.m on V. Simoncini's Website
######################################
workspace()
include("../src/LargeMatrixEquations.jl")
using LargeMatrixEquations

nh = 20
n = nh^2
T = diagm(ones(nh)*2)-diagm(ones(nh-1),-1)-diagm(ones(nh-1),1)
I = eye(nh,nh)
A = sparse(-(kron(T,I)+kron(I,T)))

srand(123)

E = eye(n)#diagm(rand(n),0)

B=randn(n,2)

m=100
tol=1e-9
tolY=1e-12
ch=true
Z,resnorm=rksm(A,B)
print(norm(A*Z*Z'*E+E*Z*Z'*A'+B*B'))

max_stars_repo_name: garrettthomaskth/LargeMatrixEquations

max_issues_repo_name: garrettthomaskth/LargeMatrixEquations

max_forks_repo_name: garrettthomaskth/LargeMatrixEquations


text:
# some implementations for convenience
# maintained by Zach Sunberg

rand(rng::AbstractRNG, t::Tuple{Bool, Bool}) = rand(rng, Bool)
rand(t::Tuple{Bool, Bool}) = rand(Bool)

support(s::AbstractVector) = s
support(s::Tuple) = s
support(r::AbstractRange) = r
support(g::Base.Generator) = g


text:
# move
"""
    move(data::Data, move_vec::Array)

Do this will move the model uniformly respect to `move_vec` in unit of Angstrom

# Example
```julia-repl
data_cell = genr_cell([1 1 1])
str = Si3N4()
data_atom = genr_atom(data_cell, str)
data_move = move(data_atom, [10 10 10])
```
"""
function move(data::Data_Unit, move_vec)
    # Motion of box
    data_new = Data_Unit(data)
    if typeof(move_vec) == Array{Int64,2}
        move_vec = conv(move_vec, 1)
    end
    data_new.data_basic.box_size = broadcast(+, data_new.data_basic.box_size, move_vec)
    # Motion of atoms
    for atom = 1 : data_new.data_basic.num_atoms
        data_new.vec_atom[atom].coord += move_vec
    end
    data_new
end

function move(data::Data_Sum, move_vec)
    # Motion of box
    data_new = Data_Sum(data)
    if typeof(move_vec) == Array{Int64,2}
        move_vec = conv(move_vec, 1)
    end
    data_new.data_basic.box_size = broadcast(+, data_new.data_basic.box_size, move_vec)
    # Motion of atoms
    for atom = 1 : data_new.data_basic.num_atoms
        data_new.vec_atom[atom].coord += move_vec
    end
    data_new
end

# select and remove

function find(vec_unit::Union{Vector{T}, Int64}, list_atom) where T <: Unit
    len = length(vec_unit)
    if typeof(vec_unit) == Int64
        return 0
    else
        atom = [vec_unit[n].atom for n = 1:len]
        list = 0
        id = 0
        for i in atom
            id += 1
            for j in i
                if sum(list_atom .== j) == 1
                    list = vcat(list, id)
                    break
                end
            end
        end
        list[2:end]
    end
end

"""
    select(data_cell::Data_Cell; mode::String, para)

Do this will return a list of id of cells in specific region controled by variables `mode` and `para`

# Arguments
- `data_cell`:
- `mode`: "cylinder" and "sphere" are supported currently.
- `para`: parameters needed to describe a region correspond to each modes

# Parameters List

- `"cylinder"`: `[radius dim]`, `dim` represents which axies cylinder will be prallel to
- `"sphere"`: `[radius]`
- `"plane"`: `[vec_plane, radius, origin]`, `vec_plane` is a vector prependicular to plane that you want to cut, `radius` is the cutoff radius cells whose projection to `vec_plane` over it will be selected, `origin` is the origin point of `vec_plane`, `[0, 0, 0]` is set defaultly

# Example
```julia-repl
data_cell = genr_cell([10 10 2])
data_select = select(data_cell, mode="cylinder", para=[3, 3])
```
"""
function select(data_cell::Data_Cell; mode::String, para)
    list_mode = split("cylinder sphere plane")
    if !in(mode, list_mode)
        error(join(["Error, mode: ",mode," is not supported."]))
    end

    center = data_cell.cell_vec ./ 2

    cell_list = 0

    if mode == "cylinder"
        for cell in 1:data_cell.num_cells
            if dist(data_cell.cell_mat[cell, :], center, para[2]) <= para[1]
                cell_list = vcat(cell_list, cell)
            end
        end
    elseif mode == "sphere"
        for cell in 1:data_cell.num_cells
            if dist(data_cell.cell_mat[cell, :], center) <= para[1]
                cell_list = vcat(cell_list, cell)
            end
        end
    elseif mode == "plane"
        vec = norm_vec(para[1:3])
        dist_proj = para[4]
        if length(para) == 4
            org = [0, 0, 0]
        else
            org = para[5:7]
        end
        for cell = 1:data_cell.num_cells
            dist_now = (data_cell.cell_mat[cell, :] .- org)' * vec
            if dist_now >= dist_proj
                cell_list = vcat(cell_list, cell)
            end
        end
    end
    if length(cell_list) == 1
        error("No cell has been selected")
    end
    return cell_list[2:end]
end

"""
    select(data_atom::Data; mode::String, para)

Do this will return a list of id of atoms in specific region controled by variables `mode` and `para`

# Arguments
- `data_atom`:
- `mode`: "cylinder" and "sphere" are supported currently.
- `para`: parameters needed to describe a region correspond to each modes

# Parameters List

- `"cylinder"`: `[radius, dim]`, `dim` represents which axies cylinder will be prallel to
- `"sphere"`: `[radius]`
- `"plane"`: `[vec_plane, radius, origin]`, `vec_plane` is a vector prependicular to plane that you want to cut, `radius` is the cutoff radius atoms whose projection to `vec_plane` over it will be selected, `origin` is the origin point of `vec_plane`, `[0, 0, 0]` is set defaultly

# Example
```julia-repl
data_cell = genr_cell([10 10 2])
str = Si3N4()
data_atom = genr_atom(data_cell, str)
data_select = select(data_atom, mode="cylinder", para=[10, 3])
```
"""
function select(data_unit::Data_Unit; mode::String, para)
    list_mode = split("cylinder sphere plane")
    if !in(mode, list_mode)
        error(join(["Error, mode: ",mode," is not supported."]))
    end

    box_size = data_unit.data_basic.box_size
    box_tilt = data_unit.data_basic.box_tilt
    center = (max(data_unit.vec_atom, "coord")-min(data_unit.vec_atom, "coord")) ./ 2
    #center[1] += box_tilt[1] + box_tilt[2]
    #center[2] += box_tilt[3]

    atom_list = 0

    if mode == "cylinder"
        for atom = 1:data_unit.data_basic.num_atoms
            if dist(data_unit.vec_atom[atom], center, para[2]) <= para[1]
                atom_list = vcat(atom_list, atom)
            end
        end
    elseif mode == "sphere"
        for atom = 1:data_unit.data_basic.num_atoms
            if dist(data_unit.vec_atom[atom], center) <= para[1]
                atom_list = vcat(atom_list, atom)
            end
        end
    elseif mode == "plane"
        vec = norm_vec(para[1:3])
        dist_proj = para[4]
        if length(para) == 4
            org = [0, 0, 0]
        else
            org = para[5:7]
        end
        for atom = 1:data_unit.data_basic.num_atoms
            dist_now = (data_unit.vec_atom[atom].coord .- org)' * vec
            #dist_now = sqrt(dist_now' * dist_now)
            if dist_now >= dist_proj
                atom_list = vcat(atom_list, atom)
            end
        end
    end

    if length(atom_list) == 1
        error("No atom has been selected")
    end
    return atom_list[2:end]
end

"""
    remove!(data_cell::Data_Cell, list_cell::Array)

Do this will remove all cells in `list_cell` from `data_cell`

# Example
```julia-repl
data_cell = genr_cell([10 10 2])
data_select = select(data_cell, mode="cylinder", para=[3, 3])
data_new = remove!(data_cell, data_select)
```
"""
function remove!(data_cell::Data_Cell, list_cell)
    # Reading Input
    cell_mat = data_cell.cell_mat
    len = size(cell_mat)[1]
    # Deleting Cells
    judge = trues(len)
    for id in list_cell
        judge .&= 1:len .!= id
    end
    data_cell.cell_mat = cell_mat[judge, :]
    data_cell.num_cells = size(data_cell.cell_mat)[1]

    # Output
    data_cell
end

"""
    remove!(data::Data, list_atom::Array)

Do this will remove all atoms in `list_atom` themselves and other infomation related to them  from `data`

# Example
```julia-repl
data_cell = genr_cell([10 10 2])
str = Si3N4()
data_atom = genr_atom(data_cell, str)
data_select = select(data_atom, mode="cylinder", para=[3, 3])
data_new = remove!(data_atom, data_select)
```
"""
function remove!(data::Data, list_atom::Array)
    fields = fieldnames(typeof(data))
    name_fields = [string(fields[n]) for n = 1:length(fields)]
    num_fields = length(fields)
    list_fields = findall(x->occursin("vec", x), name_fields)

    for field in  list_fields
        # Find all elements that need to be removed
        list_id = find(getfield(data, fields[field]), list_atom)

        if list_id == 0
            continue
        else
            # Changing # of specifc field
            num_ids = length(list_id)
            para_now = name_fields[field][findall(x->in('_', x), name_fields[field])[1]+1 : end]
            para_now = Meta.parse(join(["num_", para_now, "s"]))
            para_result = getfield(data.data_basic, para_now) - num_ids
            setfield!(data.data_basic, para_now, para_result)

            # Changing Vector of each field
            setfield!(data, fields[field], remove!(getfield(data, fields[field]), list_id)) # remove!(vec::Vector{T} , id::Array) where T<:Data_Unit
        end
    end
    sort_data!(data, list_atom)
end

function remove!(vec::Union{Vector{T}, Int64}, id::Array) where T <: Unit
    len = length(vec)
    judge = trues(len)
    for i in id
        judge .&= 1 : len .!= i
    end
    vec[judge]
end

function remove!(mat; id=1, dim=1)
    if dim == 1
        result = mat[1:end .!= id, :]
    elseif dim == 2
        result = mat[:, 1:end .!= id]
    else
        error("Error, dim should be 1 or 2, representing row or column respectively!")
    end
    result
end

"""
    remove_else!(data::Data, list_atom::Array)

Do this will remove all atoms but atoms in `list_atom` themselves and other infomation related to them  from `data`

# Example
```julia-repl
data_cell = genr_cell([10 10 2])
str = Si3N4()
data_atom = genr_atom(data_cell, str)
data_select = select(data_atom, mode="cylinder", para=[3, 3])
data_new = remove_else!(data_atom, data_select)
```
"""
function remove_else!(data::Data, list_atom::Array)
    atom_all = get_data(data.vec_atom, "atom")
    atom_else = setdiff(atom_all, list_atom) # Notice: order of `atom_all` and `atom_else` is vital
    remove!(data, atom_else)
end


text:
module EMLPIDs
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...types.EMLtypes
	import EMLtypes.length
	include("PIDs/PID.jl")
	include("PIDs/FirstOrder.jl")
	include("PIDs/StepSignal.jl")
	include("PIDs/ConstantSignal.jl")
end


text:
using StaticArrays
using Base: ReinterpretArray, ReshapedArray
using LinearAlgebra: norm

import Base.convert

export mat, vecs, xyz
export SVec, SMat, JVec, JVecF, JMat, JMatF
export maxdist, maxnorm

"typealias fora static vector type; currently `StaticArrays.SVector`"
const SVec = SVector
"typealias fora static matrix type; currently `StaticArrays.SMatrix`"
const SMat = SMatrix

"`JVec{T}` : 3-dimensional immutable vector"
const JVec{T} = SVec{3,T}
const JVecF = JVec{Float64}
const JVecI = JVec{Int}

"`JMat{T}` : 3 × 3 immutable matrix"
const JMat{T} = SMatrix{3,3,T,9}
const JMatF = SMatrix{3,3,Float64,9}
const JMatI = JMat{Int}

#

"""
`vecs(V)` : convert (as reference) a `3 x N` matrix or `3*N` vector representing
N vectors (e.g. positions or forces) in R³ to a vector of `JVec` vectors.

If `V` is obtained from a call to  `mat(X)` where `X::Vector{JVec{T}}` then
`vecs(V) === X`.
"""
vecs(V::AbstractMatrix{T}) where {T} = (@assert size(V,1) == 3;
                                        reinterpret(JVec{T}, vec(V)))
vecs(V::AbstractVector{T}) where {T} = reinterpret(JVec{T}, V)
vecs(M::ReshapedArray{T,2,ReinterpretArray{T,1,JVec{T},Array{JVec{T},1}},Tuple{}}
    ) where {T} = M.parent.parent

"""
`mat(X)`: convert (as reference) a vector of
`SVec` to a 3 x N matrix representing those.

### Usage:
```
X = positions(at)          # returns a Vector{<: JVec}
M = positions(at) |> mat   # returns an AbstractMatrix
X === vecs(M)
```
"""
mat(V::AbstractVector{SVec{N,T}}) where {N,T} =
      reshape( reinterpret(T, V), (N, length(V)) )
# mat(X::Base.ReinterpretArray) = reshape(X.parent, 3, :)

"""
convert a Vector{JVec{T}} or a 3 x N Matrix{T} into a 3-tuple
(x, y, z). E.g.,
```
x, y, z = positions(at) |> xyz
```
a short-cut is
```
x, y, z = xyz(at)
```
Conversely, `set_positions!(at, x, y, z)` is also allowed.
"""
function xyz(V::AbstractMatrix)
   @assert size(V, 1) == 3
   return (view(V, 1, :), view(V, 2, :), view(V, 3, :))
end

xyz(V::AbstractVector{JVec{T}}) where {T} = xyz(mat(V))

"""
`maxdist(A, B): ` maximum of distances between two ordered sets of objects,
typically positions. `maximum(norm(a - b) for (a,b) in zip(A,B))`
"""
function maxdist(x::AbstractArray, y::AbstractArray)
   @assert length(x) == length(y)
   return maximum( norm(a-b)  for (a,b) in zip(x,y) )
end

"""
`maxdist(X, y): ` maximum of distance of `y` to all items of `X`;
`maximum( norm(x - y) for x in X)`
"""
maxdist(X::AbstractArray{T}, y::T) where {T} = maximum(norm(x - y) for x in X)
maxdist(x::T, X::AbstractArray{T}) where {T} = maxdist(X, x)

"`maximum(norm(y) for y in x);` typically, x is a vector of forces"
maxnorm(X::AbstractVector) = maximum( norm(x) for x in X )


text:
struct IcwbScoreTable
        TOTAL_INSERTIONS
        TOTAL_DELETIONS
        TOTAL_SUBSTITUTIONS
        TOTAL_NCHANGE
        TOTAL_TRUE_WORD_COUNT
        TOTAL_TEST_WORD_COUNT
        TOTAL_TRUE_WORDS_RECALL
        TOTAL_TEST_WORDS_PRECISION
        F_MEASURE
        OOV_Rate
        OOV_Recall_Rate
        IV_Recall_Rate

        n::Int
        source::String
###	D:\KongYiji.jl\src\KongYiji.jl\..\..\test\as_kongyiji_tokens.txt	32323	2437	38777	73537	122610	152496	0.664	0.534	0.592	0.043	0.113	0.689

end

function show(io::IO, tb::IcwbScoreTable)
        fnames = fieldnames(IcwbScoreTable)
        nf = length(fnames)
        mat = Matrix{Any}(undef, (nf - 2, 2))
        for i = 1:nf - 2
                val = getfield(tb, fnames[i]) / tb.n
                val = i <= 6 ? Int(round(val)) : trunc(val; digits=3)
                mat[i,:] = [fnames[i], val]
        end
        show(io, UselessTable(mat; heads=["icwb $(tb.source) score table $(tb.n) combined"], cnames=["measure", "score"]))
end

+(a::IcwbScoreTable, b::IcwbScoreTable) = IcwbScoreTable(map(x->getfield(a, x) + getfield(b, x), fieldnames(IcwbScoreTable)[:end-1])..., a.source)

function IcwbScoreTable(dict::Set{String}, truths::Vector{Vector{SubString{String}}}, outputs::Vector{Vector{String}}, source::String)
        n = length(truths)
        maxl = 2000 #todo
        dp = Matrix{Int}(undef, (maxl, maxl))
        prev = Matrix{Tuple{Int, Int}}(undef, (maxl, maxl))
        nins, ndel, nsub, ntruew, ntestw, nhits, noov, oov_hits, niv, iv_hits = repeat([0], 10)
        for (a, b) in zip(truths, outputs)
                na, nb = length(a), length(b)
                for i in 1:na, j in 1:nb
                        dp[i,j] = 0
                        if a[i] == b[j]
                                dp[i,j] = get(dp, (i - 1,j - 1), 0) + 1
                                prev[i,j] = (i - 1,j - 1)
                        end
                        if dp[i,j] < get(dp, (i - 1,j), 0)
                                dp[i,j] = get(dp, (i - 1,j), 0)
                                prev[i,j] = get(prev, (i - 1,j), (0,0))
                        end
                        if dp[i,j] < get(dp, (i,j - 1), 0)
                                dp[i,j] = get(dp, (i,j - 1), 0)
                                prev[i,j] = get(prev, (i,j - 1), (0,0))
                        end
                end
                vs = [(max(na, nb), max(na, nb)), (na, nb)]
                ia, ib = na, nb
                while true
                        ia, ib = vs[end]
                        if ia == 0 || ib == 0 break end
                        push!(vs, prev[ia,ib])
                end
                push!(vs, (0,0))
                reverse!(vs)
                for i = 2:length(vs)
                        da, db = vs[i] .- vs[i - 1]
                        if da == db nsub += da - 1
                        elseif da < db nins += db - da; nsub += db - da
                        else ndel += da - db; nsub += da - db
                        end
                end
                ntruew += na
                ntestw += nb
                nhits += dp[na,nb]
                for wa in a !in(wa, dict) ? noov += 1 : niv += 1 end
                for v in vs[2:end-1] all(v .> 0) && !in(b[v[2]], dict) ? oov_hits += 1 : iv_hits += 1 end
        end
        recall, precision = nhits / ntruew, nhits / ntestw
        F1 = f1(precision, recall)
        oov_rate = noov / ntruew
        oov_recall = oov_hits / noov
        iv_recall = iv_hits / niv
        return IcwbScoreTable(nins, ndel, nsub, (nins+ndel+nsub), ntruew, ntestw, recall, precision, F1, oov_rate, oov_recall, iv_recall, 1, source)
end


text:
##################################
## capture un-implemented cases ##
##################################

function pdf(cop::PairCop,
             u1::FloatVec, u2::FloatVec)
    error("pdf function is not yet implemented for this type of copula")
end

function cdf(cop::PairCop,
             u1::FloatVec, u2::FloatVec)
    error("cdf function is not yet implemented for this type of copula")
end

function hfun(cop::PairCop,
              u1::FloatVec, u2::FloatVec)
    error("h-function is not yet implemented for this type of copula")
end

function vfun(cop::PairCop,
              u2::FloatVec, u1::FloatVec)
    error("v-function is not yet implemented for this type of copula")
end

function hinv(cop::PairCop,
              u1::FloatVec, u2::FloatVec)
    error("inverse h-function is not yet implemented for this type of copula")
end

function vinv(cop::PairCop,
              u2::FloatVec, u1::FloatVec)
    error("inverse v-function is not yet implemented for this type of copula")
end

################################
## handle single observations ##
################################

funcs = [:pdf,
         :cdf,
         :hfun,
         :vfun,
         :hinv,
         :vinv]

macro treatSingleObs(f)
    esc(quote
        function $(f)(cop::PairCop, u1::Float64, u2::Float64)
            return $(f)(cop, [u1], [u2])
        end
        function $(f)(cop::PairCop, u1::Float64, u2::FloatVec)
            return $(f)(cop, u1*ones(length(u2)), u2)
        end
        function $(f)(cop::PairCop, u1::FloatVec, u2::Float64)
            return $(f)(cop, u1, u2*ones(length(u1)))
        end
    end)
end

for f in funcs
        eval(macroexpand(:(@treatSingleObs($f))))
end


text:
function gradfunc{T <: Number}(
    f::Function,
    x::T;
    method::Symbol = :finite,
    direction::Symbol = :central,
) # -> Function
    if method == :finite
        function g(x::T)
            return finitegrad(f, x, direction)
        end
    elseif method == :ad
        function g(x::T)
            return autograd(f, x)
        end
    else
        throw(ArgumentError("method must be :finite or :ad"))
    end
    return g
end

function gradfunc!{T <: Number}(
    f::Function,
    x::Vector{T};
    method::Symbol = :finite,
    direction::Symbol = :central,
) # -> Function
    if method == :finite
        function g!(gr::Vector{T}, x::Vector{T})
            return finitegrad!(gr, f, x, direction)
        end
    elseif method == :ad
        x_d = Array(Dual{T}, length(x))
        function g!(gr::Vector{T}, x::Vector{T})
            return autograd!(gr, x_d, f, x)
        end
    else
        throw(ArgumentError("method must be :finite or :ad"))
    end
    return g!
end

function gradfunc{T <: Number}(
    f::Function,
    x::Vector{T};
    method::Symbol = :finite,
    direction::Symbol = :central,
) # -> Function
    if method == :finite
        function g(x::Vector{T})
            gr = similar(x)
            finitegrad!(gr, f, x, direction)
            return gr
        end
    elseif method == :ad
        x_d = Array(Dual{T}, length(x))
        function g(x::Vector{T})
            gr = similar(x)
            autograd!(gr, x_d, f, x)
            return gr
        end
    else
        throw(ArgumentError("method must be :finite or :ad"))
    end
    return g
end


text:
# Robertson chemical reaction:
#
# Adapted from
# http://www.unige.ch/~hairer/testset/testset.html
# reference solution from
# http://www.dm.uniba.it/~testset/problems/rober.php
#
# See also:
# http://octave.sourceforge.net/odepkg/function/odepkg_testsuite_robertson.html
#
# Note that this runs over large tspan = [0,1e11] with the first step being < 1!

export rober

rober = let
    tcname = :rober
    T = Float64 # the datatype used
    Tarr = Matrix
    dof = 3 # degrees of freedom
    dae=0
    # stiffness of system, one of the three constants
    stiffness = [nonstiff, mildlystiff, stiff][3]
    function fn!(t,y,dydt)
        # the ode function
        dydt[1] = -0.04*y[1] + 1.0e4*y[2]*y[3]
        dydt[3] = 3.0e7 *y[2]*y[2]
        dydt[2] = -dydt[1]-dydt[3]
        nothing
    end
    # initializes storage for y:
    fn!( ; T_::Type=T, dof_=dof) = zeros(T_,dof_)

    function jac!(t,y,dfdy)
        # the Jacobian of f
        prod1 = 1.0e4*y[2]
        prod2 = 1.0e4*y[3]
        prod3 = 6.0e7*y[2]
        dfdy[1,1] = -0.04
        dfdy[1,2] = prod2
        dfdy[1,3] = prod1
        dfdy[2,1] = 0.04
        dfdy[2,2] = -prod2-prod3
        dfdy[2,3] = -prod1
        dfdy[3,1] = 0.0
        dfdy[3,2] = prod3
        dfdy[3,3] = 0.0
        return nothing
    end
    jac!( ; T_::Type=T, dof_=dof) = zeros(T_,dof_,dof_)  # if the problem is large better return an appropriate sparse matrix

    # mass matrix M(t,y) dydt = f(t,y)
    mass! = nothing

    ic = T[1,0,0] # vector of initial conditions
    tspan = T[0, 1e11] # integration interval
    refsol = T[0.2083340149701255e-07,
               0.8333360770334713e-13,
               0.9999999791665050] # reference solution at tspan[2]
    refsolinds = trues(dof)
    scd_absinds = falses(dof)

    tc = TestCaseExplicit{tcname, T, Tarr}(
                             stiffness,
                             dae,
                             dof,
                             fn!,
                             jac!,
                             mass!,
                             ic,
                             tspan,
                             refsol,
                             refsolinds,
                             scd_absinds)

    tc_all[tcname] = tc
    tc_stiff[tcname] = tc
end


text:
ndiff(A::Array, n::Integer) = n < 1 ? A : diff(ndiff(A, n-1))

s = [90, 47, 58, 29, 22, 32, 55, 5, 55, 73]
println.(collect(ndiff(s, i) for i in 0:9))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using LinearAlgebra, StatsBase, CoupledFields
# using Base.Test

# write your own tests here
# @test 1 == 1

function simfields(t::Vector{Float64}, p::Int64, σₑ::Float64)
    n = length(t)
    Xa = sin.(2π*repeat(t, inner=(1,p))/5.0) * diagm(0=>rand([-1, 1], p))
    Xb = sin.(2π*repeat(t, inner=(1,p))/2.5 .+ 1.5π) * diagm(0=>rand([-1, 1], p))
    Xa = StatsBase.zscore(Xa,1) .+ 0.05*randn(n, p)
    Xb = StatsBase.zscore(Xb,1) .+ σₑ*randn(n, p)
    return Xa, Xb
end

t1 = collect((2001+1/24):(1/12):2015)
lpfield, hpfield = simfields(t1, 10, 1.0)
Z = InputSpace(lpfield, hpfield, [1.1,1.1])
kpars = GaussianKP(Z.X)
model = gKCCA([1.0 -5.0 2], Z.X, Z.Y, kpars)


text:
import UncertainData: RandomSequences

"""
    RandomSequencesTest(test::CausalityTest, sequences_resampling::RandomSequences)

A causality test where the `test` is applied to multiple independent draws of
the datasets, where each draw from a randomly selected consecutive chunk of
points.

## Examples

Say we want do do a cross mapping test on a set of time series with `N = 200`
points. To get some ensemble statistics, we want to do an ensemble analysis
on 150 randomly selected chunks of the time series with lengths ranging from
`N-30` to `N-1`.

```julia
# A cross mapping test with default parameters (don't use default parameters in a real analysis!)
cm_test = CrossMappingTest()

# Combine with the random sequence subsamling
rtest = RandomSequencesTest(, RandomSequences(150, N-30:N-1))
```
"""
struct RandomSequencesTest{CT <: CausalityTest, RS <: RandomSequences}
    test::CT
    sequences_resampling::RS
end

export RandomSequencesTest


text:
#=
External Sampling Counterfactual Regret Minimization
- "sample only the actions of the opponent and chance (those choices external to the player)"
=#
using Random

# try `Ref{Int}` for `a_idx` to keep immutable
mutable struct MCInfoState <: AbstractInfoState
    σ::Vector{Float64}
    r::Vector{Float64}
    s::Vector{Float64}
    _tmp_σ::Vector{Float64}
    a_idx::Int
end

function MCInfoState(L::Int)
    return MCInfoState(
        fill(1/L, L),
        zeros(L),
        fill(1/L,L),
        fill(1/L,L),
        0
    )
end

mutable struct DebugMCInfoState <: AbstractInfoState
    σ::Vector{Float64}
    r::Vector{Float64}
    s::Vector{Float64}
    _tmp_σ::Vector{Float64}
    hist::Vector{Vector{Float64}}
    a_idx::Int
end

function DebugMCInfoState(L::Int)
    return DebugMCInfoState(
        fill(1/L, L),
        zeros(L),
        fill(1/L, L),
        fill(1/L, L),
        Vector{Float64}[],
        0
    )
end

struct ESCFRSolver{K,G,I} <: AbstractCFRSolver{K,G,I}
    I::Dict{K, I}
    game::G
end

"""
Generate random action index from info state according to strategy σ
"""
function Random.rand(rng::AbstractRNG, I::AbstractInfoState)
    σ = I.σ
    t = rand(rng)
    i = 1
    cw = σ[1]
    while cw < t && i < length(σ)
        i += 1
        @inbounds cw += σ[i]
    end
    return i
end

Random.rand(I::AbstractInfoState) = rand(Random.GLOBAL_RNG, I)

"""
    `ESCFRSolver(game::Game{H,K}; debug::Bool=false)`

Instantiate external sampling CFR solver with some `game`.

If `debug=true`, record history of strategies over training period, allowing
for training history of individual information states to be plotted with
`Plots.plot(is::DebugInfoState)`

"""
function ESCFRSolver(game::Game{H,K}; debug::Bool=false) where {H,K}
    if debug
        return ESCFRSolver(Dict{K, DebugMCInfoState}(), game)
    else
        return ESCFRSolver(Dict{K, MCInfoState}(), game)
    end
end

# TODO: Validate regret update. Not entirely confident in code nor MCCFR paper eqn.
function CFR(solver::ESCFRSolver, h, i, t, π_1, π_2)
    game = solver.game
    if isterminal(game, h)
        return utility(game, i, h)
    elseif player(game, h) === 0 # chance player
        a = chance_action(game, h)
        h′ = next_hist(game,h,a)
        return CFR(solver, h′, i, t, π_1, π_2)
    end

    I = infoset(solver, h)
    A = actions(game, h)

    v_σ = 0.0

    if player(game, h) === i
        v_σ_Ia = I._tmp_σ
        π_i = i == 1 ? π_1 : π_2
        for (k,a) in enumerate(A)
            h′ = next_hist(game, h, a)
            if i === 1
                v_σ_Ia[k] = CFR(solver, h′, i, t, I.σ[k]*π_1, π_2)
            else
                v_σ_Ia[k] = CFR(solver, h′, i, t, π_1, I.σ[k]*π_2)
            end
            v_σ += I.σ[k]*v_σ_Ia[k]
        end
        for (k,a) in enumerate(A)
            I.r[k] += (1 - I.σ[k])*(v_σ_Ia[k] - v_σ)
            I.s[k] += π_i*I.σ[k]
        end
    else
        a_idx = I.a_idx
        iszero(a_idx) && (a_idx = rand(I))
        I.a_idx = a_idx
        a = A[a_idx]
        h′ = next_hist(game, h, a)
        v_σ = CFR(solver, h′, i, t, π_1, π_2)
    end

    return v_σ
end

function train!(solver::ESCFRSolver{K,G,INFO}, N::Int; show_progress::Bool=false, cb=()->()) where {K,G,INFO<:MCInfoState}
    regret_match!(solver)
    ih = initialhist(solver.game)
    prog = Progress(N; enabled=show_progress)
    for t in 1:N
        for i in 1:2
            CFR(solver, ih, i, t, 1.0, 1.0)
        end
        for I in values(solver.I)
            regret_match!(I)
            I.a_idx = 0
        end
        cb()
        next!(prog)
    end
    finalize_strategies!(solver)
end

function train!(solver::ESCFRSolver{K,G,INFO}, N::Int; show_progress::Bool=false, cb=()->()) where {K,G,INFO<:DebugMCInfoState}
    regret_match!(solver)
    ih = initialhist(solver.game)
    prog = Progress(N; enabled=show_progress)
    for t in 1:N
        for i in 1:2
            CFR(solver, ih, i, t, 1.0, 1.0)
        end
        for I in values(solver.I)
            regret_match!(I)
            push!(I.hist, copy(I.s) ./ sum(I.s))
            I.a_idx = 0
        end
        cb()
        next!(prog)
    end
    finalize_strategies!(solver)
end

max_stars_repo_name: WhiffleFish/HelloCFR

max_issues_repo_name: WhiffleFish/HelloCFR

max_forks_repo_name: WhiffleFish/HelloCFR


text:
@testset "1502.can-make-arithmetic-progression-from-sequence.jl" begin
    @test can_make_arith([3, 5, 1])
    @test !can_make_arith([1, 2, 4])
end


text:
import Pkg
import UUIDs

@inline function get_all_versions(pkg_path::String)::Vector{VersionNumber}
    versions_toml = Pkg.TOML.parsefile(joinpath(pkg_path, "Versions.toml"))
    all_versions = VersionNumber.(collect(keys(versions_toml)))
    unique!(all_versions)
    sort!(all_versions)
    return all_versions
end

@inline function _get_latest_version(versions::AbstractVector{VersionNumber})
    return maximum(versions)
end

@inline function _get_latest_zero_version(versions::AbstractVector{VersionNumber})
    all_zero_versions = filter((x) -> (x < v"1"), versions)
    if isempty(all_zero_versions)
        return nothing
    else
        latest_zero_version = maximum(all_zero_versions)
        always_assert(latest_zero_version < v"1", "latest_zero_version < v\"1\"")
        return latest_zero_version
    end
end

@inline function parse_registry(registry_paths::Vector{String})
    pkg_to_path = Dict{Package, String}()
    pkg_to_num_versions = Dict{Package, Int}()
    pkg_to_latest_version = Dict{Package, VersionNumber}()
    pkg_to_latest_zero_version = Dict{Package, Union{VersionNumber, Nothing}}()
    parse_registry!(pkg_to_path,
                    pkg_to_num_versions,
                    pkg_to_latest_version,
                    pkg_to_latest_zero_version,
                    registry_paths)
    return pkg_to_path,
           pkg_to_num_versions,
           pkg_to_latest_version,
           pkg_to_latest_zero_version
end

@inline function parse_registry!(pkg_to_path::AbstractDict{Package, String},
                                 pkg_to_num_versions::AbstractDict{Package, Int},
                                 pkg_to_latest_version::AbstractDict{Package, VersionNumber},
                                 pkg_to_latest_zero_version::AbstractDict{Package, <:Union{VersionNumber, Nothing}},
                                 registry_paths::Vector{String})
    for path in registry_paths
        parse_registry!(pkg_to_path,
                        pkg_to_num_versions,
                        pkg_to_latest_version,
                        pkg_to_latest_zero_version,
                        path)
    end
    return pkg_to_path,
           pkg_to_num_versions,
           pkg_to_latest_version,
           pkg_to_latest_zero_version
end

@inline function parse_registry!(pkg_to_path::AbstractDict{Package, String},
                                 pkg_to_num_versions::AbstractDict{Package, Int},
                                 pkg_to_latest_version::AbstractDict{Package, VersionNumber},
                                 pkg_to_latest_zero_version::AbstractDict{Package, <:Union{VersionNumber, Nothing}},
                                 registry_path::String)
    registry = Pkg.TOML.parsefile(joinpath(registry_path, "Registry.toml"))
    packages = registry["packages"]
    for p in packages
        name = p[2]["name"]
        pkg_path = joinpath(registry_path, p[2]["path"])
        if !is_jll(name)
            all_versions = get_all_versions(pkg_path)
            num_versions = length(all_versions)
            latest_version = _get_latest_version(all_versions)
            latest_zero_version = _get_latest_zero_version(all_versions)
            pkg = Package(name)
            pkg_to_path[pkg] = pkg_path
            pkg_to_num_versions[pkg] = num_versions
            pkg_to_latest_version[pkg] = latest_version
            pkg_to_latest_zero_version[pkg] = latest_zero_version
        end
    end
    return nothing
end