text:
import ...timeevolution: dschroedinger!, dschroedinger_dynamic!, check_schroedinger

"""
    stochastic.schroedinger(tspan, state0, H, Hs[; fout, ...])

Integrate stochastic Schrödinger equation.

# Arguments
* `tspan`: Vector specifying the points of time for which the output should
        be displayed.
* `psi0`: Initial state as Ket.
* `H`: Deterministic part of the Hamiltonian.
* `Hs`: Stochastic part(s) of the Hamiltonian (either an operator or a vector
        of operators).
* `fout=nothing`: If given, this function `fout(t, state)` is called every time
        an output should be displayed. ATTENTION: The given state is neither
        normalized nor permanent!
* `normalize_state=false`: Specify whether or not to normalize the state after
        each time step taken by the solver.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function schroedinger(tspan, psi0::T, H::AbstractOperator{B,B}, Hs::Vector;
                fout=nothing,
                normalize_state=false,
                calback=nothing,
                kwargs...) where {B,T<:Ket{B}}
    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)

    n = length(Hs)
    dstate = copy(psi0)
    x0 = psi0.data
    state = copy(psi0)

    # TODO: replace checks by dispatch
    for h=Hs
        @assert isa(h, AbstractOperator{B,B})
    end

    dschroedinger_determ(t, psi, dpsi) = dschroedinger!(dpsi, H, psi)
    dschroedinger_stoch(dx, t, psi, dpsi, n) = dschroedinger_stochastic(dx, psi, Hs, dpsi, n)

    if normalize_state
        norm_func(u, t, integrator) = normalize!(u)
        ncb = DiffEqCallbacks.FunctionCallingCallback(norm_func;
                 func_everystep=true,
                 func_start=false)
    else
        ncb = nothing
    end

    integrate_stoch(tspan_, dschroedinger_determ, dschroedinger_stoch, x0, state, dstate, fout, n;
                    ncb=ncb,
                    kwargs...)
end
schroedinger(tspan, psi0::Ket{B}, H::AbstractOperator{B,B}, Hs::AbstractOperator{B,B}; kwargs...) where B = schroedinger(tspan, psi0, H, [Hs]; kwargs...)

"""
    stochastic.schroedinger_dynamic(tspan, state0, fdeterm, fstoch[; fout, ...])

Integrate stochastic Schrödinger equation with dynamic Hamiltonian.

# Arguments
* `tspan`: Vector specifying the points of time for which the output should
        be displayed.
* `psi0`: Initial state.
* `fdeterm`: Function `f(t, psi, u) -> H` returning the deterministic
        (time- or state-dependent) part of the Hamiltonian.
* `fstoch`: Function `f(t, psi, u, du) -> Hs` returning a vector that
        contains the stochastic terms of the Hamiltonian.
* `fout=nothing`: If given, this function `fout(t, state)` is called every time
        an output should be displayed. ATTENTION: The given state is neither
        normalized nor permanent!
* `noise_processes=0`: Number of distinct white-noise processes in the equation.
        This number has to be equal to the total number of noise operators
        returned by `fstoch`. If unset, the number is calculated automatically
        from the function output.
        NOTE: Set this number if you want to avoid an initial calculation of
        the function output!
* `normalize_state=false`: Specify whether or not to normalize the state after
        each time step taken by the solver.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function schroedinger_dynamic(tspan, psi0::Ket, fdeterm, fstoch;
                fout=nothing, noise_processes::Int=0,
                normalize_state=false,
                kwargs...)
    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)

    if noise_processes == 0
        fs_out = fstoch(0.0, psi0)
        n = length(fs_out)
    else
        n = noise_processes
    end

    dstate = copy(psi0)
    x0 = psi0.data
    state = copy(psi0)

    dschroedinger_determ(t, psi, dpsi) = dschroedinger_dynamic!(dpsi, fdeterm, psi, t)
    dschroedinger_stoch(dx, t, psi, dpsi, n) = dschroedinger_stochastic(dx, t, psi, fstoch, dpsi, n)

    if normalize_state
        norm_func(u, t, integrator) = normalize!(u)
        ncb = DiffEqCallbacks.FunctionCallingCallback(norm_func;
                 func_everystep=true,
                 func_start=false)
    else
        ncb = nothing
    end

    integrate_stoch(tspan, dschroedinger_determ, dschroedinger_stoch, x0, state,
            dstate, fout, n;
            ncb=ncb,
            kwargs...)
end

function dschroedinger_stochastic(dx::AbstractVector, psi, Hs, dpsi, index)
    recast!(dpsi,dx)
    dschroedinger!(dpsi, Hs[index], psi)
end
function dschroedinger_stochastic(dx::AbstractMatrix, psi, Hs, dpsi, n)
    for i=1:n
        dx_i = @view dx[:, i]
        recast!(dpsi,dx_i)
        dschroedinger!(dpsi, Hs[i], psi)
        recast!(dx_i,dpsi)
    end
end
function dschroedinger_stochastic(dx, t, psi, f, dpsi, n)
    ops = f(t, psi)
    if QO_CHECKS[]
        @inbounds for h=ops
            check_schroedinger(psi, h)
        end
    end
    dschroedinger_stochastic(dx, psi, ops, dpsi, n)
end

recast!(x::SubArray,psi::Ket) = (x .= psi.data)
recast!(psi::Ket,x::SubArray) = (psi.data = x)


text:
function get_axiom_network()
    network = Network()
    update_inputs_outputs!(network, [get_axiom_node()])
    return network
end
function get_axiom_node()
    n1 = Node(AXIOM_NODES[1], [], [], [], [])
    n2 = Node(AXIOM_NODES[2], [], [], [], [])
    n3 = Node(AXIOM_NODES[3], [], [], [], [])

    n1.type = InputNode()
    n3.type = OutputNode()

    push!(n1.out_excitatory, n2)
    push!(n2.in_excitatory, n1)

    push!(n2.out_excitatory, n3)
    push!(n3.in_excitatory, n2)

    return n1
end

ext: jl

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
module ForestFire

include("Forest.jl")
include("SpreadFire.jl")
include("Animation.jl")

export createforest
export burn
export plotforest

end

max_stars_repo_name: TheSquareRoot/ForestFire

max_issues_repo_name: TheSquareRoot/ForestFire

max_forks_repo_name: TheSquareRoot/ForestFire


text:
module SSLTestTiling
using SimpleSDMLayers
using Test

# hcat / vcat
l1 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=0.0, top=10.0)
l2 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=10.0, top=20.0)
l3 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=10.0, right=20.0, bottom=0.0, top=10.0)
l4 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=10.0, right=20.0, bottom=10.0, top=20.0)
l5 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=20.0, bottom=0.0, top=20.0)

ml1 = hcat(l1, l3)
vl1 = vcat(l1, l2)
ml2 = hcat(l2, l4)
vl2 = vcat(l3, l4)

vml = vcat(ml1, ml2)
mvl = hcat(vl1, vl2)

@test all(vml.grid == mvl.grid)

for l in (vml, mvl)
   @test all(l.grid == l5.grid)
   @test size(l) == size(l5)
   @test stride(l) == stride(l5)
   @test longitudes(l) == longitudes(l5)
   @test latitudes(l) == latitudes(l5)
   @test l.left == l5.left
   @test l.right == l5.right
   @test l.bottom == l5.bottom
   @test l.top == l5.top
end

end


text:
# Functions to parse strings to intervals

"""
    parse{T}(DecoratedInterval{T}, s::AbstractString)

Parse a string of the form `"[a, b]_dec"` as a `DecoratedInterval`
with decoration `dec`.
"""
function parse(::Type{DecoratedInterval{T}}, s::AbstractString) where T
    m = match(r"(\[.*\])(\_.*)?", s)

    if m == nothing  # matched
        throw(ArgumentError("Unable to process string $x as decorated interval"))

    end

    interval_string, decoration_string = m.captures
    interval = parse(Interval{T}, interval_string)

    # type unstable:
    if decoration_string == nothing
        decoration_string = "_com"
    end

    decoration_symbol = Symbol(decoration_string[2:end])
    decoration = getfield(IntervalArithmetic, decoration_symbol)

    return DecoratedInterval(interval, decoration)

end

"""
    parse{T}(Interval{T}, s::AbstractString)

Parse a string as an interval. Formats allowed include:
- "1"
- "[1]"
- "[3.5, 7.2]"
- "[-0x1.3p-1, 2/3]"  # use numerical expressions
"""
function parse(::Type{Interval{T}}, s::AbstractString) where T

    # Check version!
    if !(occursin("[", s))  # string like "3.1"

        m = match(r"(.*)±(.*)", s)
        if m != nothing
            a = parse(T, strip(m.captures[1]))
            b = parse(T, strip(m.captures[2]))

            return a ± b
        end

        a = parse(T, s, RoundDown)
        b = parse(T, s, RoundUp)

        return Interval(a, b)

    end

    # match string of form [a, b]_dec:
    m = match(r"\[(.*),(.*)\]", s)

    if m != nothing  # matched
        lo, hi = m.captures

    else

        m = match(r"\[(.*)\]", s)  # string like "[1]"

        if m == nothing
            throw(ArgumentError("Unable to process string $s as interval"))
        end

        if m.captures[1] == "Empty"
            return emptyinterval(T)
        end

        lo = m.captures[1]
        hi = lo

    end

    expr1 = Meta.parse(lo)
    expr2 = Meta.parse(hi)

    interval = eval(make_interval(T, expr1, [expr2]))

    return interval

end


text:
using jcontrol
using jcontrol.visuals
import jcontrol: Track

plt = draw(:arena)

track = Track(;start_waypoint=4, keep_n_waypoints=-1)
draw!(track)

display(plt)

max_stars_repo_name: FedeClaudi/LocomotionControl

max_issues_repo_name: FedeClaudi/LocomotionControl

max_forks_repo_name: FedeClaudi/LocomotionControl


text:
let
    a,b = Mem.info()
    # NOTE: actually testing this is pretty fragile on CI
    #=@test a == =# VectorEngine.available_memory()
    #=@test b == =# VectorEngine.total_memory()
end

# dummy data
T = UInt32
N = 5
data = rand(T, N)
nb = sizeof(data)

# buffers are untyped, so we use a convenience function to get a typed pointer
# we prefer to return a device pointer (for managed buffers) to maximize VectorEngine coverage
typed_pointer(buf::Mem.Device, T) = convert(VEPtr{T}, buf)
typed_pointer(buf::Mem.Host, T)   = convert(Ptr{T},   buf)

# allocations and copies
for srcTy in [Mem.Device, Mem.Host],
    dstTy in [Mem.Device, Mem.Host]

    local dummy = Mem.alloc(srcTy, 0)
    Mem.free(dummy)

    src = Mem.alloc(srcTy, nb)
    unsafe_copyto!(typed_pointer(src, T), pointer(data), N)

    dst = Mem.alloc(dstTy, nb)
    unsafe_copyto!(typed_pointer(dst, T), typed_pointer(src, T), N)

    ref = Array{T}(undef, N)
    unsafe_copyto!(pointer(ref), typed_pointer(dst, T), N)

    @test data == ref

    if isa(src, Mem.Device)
        Mem.set!(typed_pointer(src, T), zero(T), N)
    end

    Mem.free(src)
    Mem.free(dst)
end

# asynchronous operations
let
    src = Mem.alloc(Mem.Device, nb)

    @test_throws ArgumentError unsafe_copyto!(typed_pointer(src, T), pointer(data), N; async=true)
    unsafe_copyto!(typed_pointer(src, T), pointer(data), N; async=true, stream=VEDefaultStream())

    Mem.set!(typed_pointer(src, T), zero(T), N; async=true, stream=VEDefaultStream())

    Mem.free(src)
end


text:
## Environment and packages
cd(@__DIR__)
using Pkg; Pkg.activate("."); Pkg.instantiate()

using OrdinaryDiffEq
using ModelingToolkit
using DataDrivenDiffEq
using LinearAlgebra, DiffEqSensitivity, Optim
using DiffEqFlux, Flux
using Plots
gr()
using JLD2, FileIO
using Statistics
# Set a random seed for reproduceable behaviour
using Random

Random.seed!(2345)

# Create a name for saving ( basically a prefix )
svname = "Scenario_2_"

## Data generation
## Lotka-Volterra equation
function lotka!(du, u, p, t)

    α, β, γ, δ = p
    du[1] = α*u[1] - β*u[2]*u[1]
    du[2] = γ*u[1]*u[2]  - δ*u[2]
end

# Define the experimental parameter
tspan = (0.0f0,6.0f0)
u0 = Float32[0.44249296,4.6280594]
p_ = Float32[1.3, 0.9, 0.8, 1.8]
prob = ODEProblem(lotka!, u0,tspan, p_)
solution = solve(prob, Vern7(), abstol=1e-6, reltol=1e-6, saveat = 0.1)

scatter(solution, alpha = 0.25)
plot!(solution, alpha = 0.5)

# Ideal data
X = Array(solution)
t = solution.t
# Add noise in terms of the mean
x̄ = mean(X, dims = 2)
noise_magnitude = Float32(1e-2)
Xₙ = X .+ (noise_magnitude*x̄) .* randn(eltype(X), size(X))

# Subsample the data in y
# We assume we have only 5 measurements in y, evenly distributed
ty = collect(t[1]:Float32(6/5):t[end])
# Create datasets for the different measurements
round(Int64, mean(diff(ty))/mean(diff(t)))
XS = zeros(eltype(X), length(ty)-1, floor(Int64, mean(diff(ty))/mean(diff(t)))+1) # All x data
TS = zeros(eltype(t), length(ty)-1, floor(Int64, mean(diff(ty))/mean(diff(t)))+1) # Time data
YS = zeros(eltype(X), length(ty)-1, 2) # Just two measurements in y

for i in 1:length(ty)-1
    idxs = ty[i].<= t .<= ty[i+1]
    XS[i, :] = Xₙ[1, idxs]
    TS[i, :] = t[idxs]
    YS[i, :] = [Xₙ[2, t .== ty[i]]'; Xₙ[2, t .== ty[i+1]]]
end

XS

scatter!(t, transpose(Xₙ))

## Define the network
# Gaussian RBF as activation
rbf(x) = exp.(-(x.^2))

# Define the network 2->5->5->5->2
U = FastChain(
    FastDense(2,5,rbf), FastDense(5,5, rbf), FastDense(5,5, rbf), FastDense(5,2)
)
# Get the initial parameters, first is linear decay
p = [rand(Float32); initial_params(U)]

# Define the hybrid model
function ude_dynamics!(du,u, p, t, p_true)
    û = U(u, p[2:end]) # Network prediction
    du[1] = p_true[1]*u[1] + û[1]
    # We assume a linear decay rate for the predator
    du[2] = -p[1]*u[2] + û[2]
end

# Closure with the known parameter
nn_dynamics!(du,u,p,t) = ude_dynamics!(du,u,p,t,p_)
# Define the problem
prob_nn = ODEProblem(nn_dynamics!,Xₙ[:, 1], tspan, p)

## Function to train the network
# Define a predictor
function predict(θ, X = Xₙ[:,1], T = t)
    Array(solve(prob_nn, Vern7(), u0 = X, p=θ,
                tspan = (T[1], T[end]), saveat = T,
                abstol=1e-6, reltol=1e-6,
                sensealg = ForwardDiffSensitivity()
                ))
end

# Multiple shooting like loss
function loss(θ)
    # Start with a regularization on the network
    l = convert(eltype(θ), 1e-3)*sum(abs2, θ[2:end]) ./ length(θ[2:end])
    for i in 1:size(XS,1)
        X̂ = predict(θ, [XS[i,1], YS[i,1]], TS[i, :])
        # Full prediction in x
        l += sum(abs2, XS[i,:] .- X̂[1,:])
        # Add the boundary condition in y
        l += abs(YS[i, 2] .- X̂[2, end])
    end

    return l
end

# Container to track the losses
losses = Float32[]

# Callback to show the loss during training
callback(θ,l) = begin
    push!(losses, l)
    if length(losses)%50==0
        println("Current loss after $(length(losses)) iterations: $(losses[end])")
    end
    false
end

## Training

# First train with ADAM for better convergence -> move the parameters into a
# favourable starting positing for BFGS
res1 = DiffEqFlux.sciml_train(loss, p, ADAM(0.1f0), cb=callback, maxiters = 200)
println("Training loss after $(length(losses)) iterations: $(losses[end])")
# Train with BFGS
res2 = DiffEqFlux.sciml_train(loss, res1.minimizer, BFGS(initial_stepnorm=0.01f0), cb=callback, maxiters = 10000)
println("Final training loss after $(length(losses)) iterations: $(losses[end])")

# Plot the losses
pl_losses = plot(1:200, losses[1:200], yaxis = :log10, xaxis = :log10, xlabel = "Iterations", ylabel = "Loss", label = "ADAM", color = :blue)
plot!(201:length(losses), losses[201:end], yaxis = :log10, xaxis = :log10, xlabel = "Iterations", ylabel = "Loss", label = "BFGS", color = :red)
savefig(pl_losses, joinpath(pwd(), "plots", "$(svname)_losses.pdf"))

# Rename the best candidate
p_trained = res2.minimizer

## Analysis of the trained network
# Plot the data and the approximation
X̂ = predict(p_trained)
# Trained on noisy data vs real solution
pl_trajectory = plot(t, transpose(X̂), ylabel = "t", xlabel ="x(t), y(t)", color = :red, label = ["UDE Approximation" nothing])
scatter!(t, X[1,:], color = :black, label = "Measurements")
ymeasurements = unique!(vcat(YS...))
tmeasurements = unique!(vcat([[ts[1], ts[end]] for ts in eachrow(TS)]...))
scatter!(tmeasurements, ymeasurements, color = :black, label = nothing, legend = :topleft)
savefig(pl_trajectory, joinpath(pwd(), "plots", "$(svname)_trajectory_reconstruction.pdf"))

# Ideal unknown interactions of the predictor
Ȳ = [-p_[2]*(X̂[1,:].*X̂[2,:])';p_[3]*(X̂[1,:].*X̂[2,:])']
# Neural network guess
Ŷ = U(X̂,p_trained[2:end])

pl_reconstruction = plot(t, transpose(Ŷ), xlabel = "t", ylabel ="U(x,y)", color = :red, label = ["UDE Approximation" nothing])
plot!(t, transpose(Ȳ), color = :black, label = ["True Interaction" nothing], legend = :topleft)
savefig(pl_reconstruction, joinpath(pwd(), "plots", "$(svname)_missingterm_reconstruction.pdf"))

# Plot the error
pl_reconstruction_error = plot(t, norm.(eachcol(Ȳ-Ŷ)), yaxis = :log, xlabel = "t", ylabel = "L2-Error", color = :red, label = nothing)
pl_missing = plot(pl_reconstruction, pl_reconstruction_error, layout = (2,1))
savefig(pl_missing, joinpath(pwd(), "plots", "$(svname)_missingterm_reconstruction_and_error.pdf"))
pl_overall = plot(pl_trajectory, pl_missing)
savefig(pl_overall, joinpath(pwd(), "plots", "$(svname)_reconstruction.pdf"))
## Symbolic regression via sparse regression ( SINDy based )

# Create a Basis
@variables u[1:2]

# Generate the basis functions, multivariate polynomials up to deg 5
# and sine
b = [polynomial_basis(u, 5); sin.(u)]
basis = Basis(b, u)

# Create an optimizer for the SINDy problem
opt = SR3(Float32(1e-2), Float32(1e-2))
# Create the thresholds which should be used in the search process
λ = Float32.(exp10.(-7:0.1:3))
# Target function to choose the results from; x = L0 of coefficients and L2-Error of the model
g(x) = x[1] < 1 ? Inf : norm(x, 2)

# Test on ideal derivative data for unknown function ( not available )
println("SINDy on partial ideal, unavailable data")
Ψ = SINDy(X̂, Ȳ, basis, λ, opt, g = g, maxiter = 10000) # Succeed
println(Ψ)
print_equations(Ψ)

# Test on uode derivative data
println("SINDy on learned, partial, available data")
Ψ = SINDy(X̂, Ŷ, basis, λ,  opt, g = g, maxiter = 10000, normalize = true, denoise = true) # Succeed
println(Ψ)
print_equations(Ψ)

# Extract the parameter
p̂ = parameters(Ψ)
println("First parameter guess : $(p̂)")

# Just the equations
b = Basis((u, p, t)->Ψ(u, [1.; 1.], t), u)

# Retune for better parameters -> we could also use DiffEqFlux or other parameter estimation tools here.
Ψf = SINDy(X̂, Ŷ, b, STRRidge(0.01), maxiter = 100, convergence_error = 1e-18) # Succeed
println(Ψf)
p̂ = parameters(Ψf)
println("Second parameter guess : $(p̂)")
println("Overall parameter guess : $(abs.([p̂; p_trained[1]]))")
println("True paramter : $(p_[2:end])")

# Define the recovered, hyrid model
function recovered_dynamics!(du,u, p, t, p_true)
    û = Ψf(u, p) # Network prediction
    du[1] = p_true[1]*u[1] + û[1]
    du[2] = -p_trained[1]*u[2] + û[2] # Learned paramter + prediction
end

# Closure with the known parameter
estimated_dynamics!(du,u,p,t) = recovered_dynamics!(du,u,p,t,p_)

estimation_prob = ODEProblem(estimated_dynamics!, u0, tspan, p̂)
estimate = solve(estimation_prob, Tsit5(), saveat = t)

# Plot
plot(solution)
plot!(estimate)

## Simulation

# Look at long term prediction
t_long = (0.0f0, 50.0f0)
estimation_prob = ODEProblem(estimated_dynamics!, u0, t_long, p̂)
estimate_long = solve(estimation_prob, Tsit5(), saveat = 0.1) # Using higher tolerances here results in exit of julia
plot(estimate_long)

true_prob = ODEProblem(lotka!, u0, t_long, p_)
true_solution_long = solve(true_prob, Tsit5(), saveat = estimate_long.t)
plot!(true_solution_long)

## Save the results
save(joinpath(pwd(), "results" ,"$(svname)recovery_$(noise_magnitude).jld2"),
    "solution", solution, "X", Xₙ, "t" , t, "neural_network" , U, "initial_parameters", p, "trained_parameters" , p_trained, # Training
    "losses", losses, "result", Ψf, "recovered_parameters", p̂, # Recovery
    "long_solution", true_solution_long, "long_estimate", estimate_long) # Estimation

## Post Processing and Plots

c1 = 3 # RGBA(174/255,192/255,201/255,1) # Maroon
c2 = :orange # RGBA(132/255,159/255,173/255,1) # Red
c3 = :blue # RGBA(255/255,90/255,0,1) # Orange
c4 = :purple # RGBA(153/255,50/255,204/255,1) # Purple

p1 = plot(t,abs.(Array(solution) .- estimate)' .+ eps(Float32),
          lw = 3, yaxis = :log, title = "Timeseries of UODE Error",
          color = [3 :orange], xlabel = "t",
          label = ["x(t)" "y(t)"],
          titlefont = "Helvetica", legendfont = "Helvetica",
          legend = :topright)

# Plot L₂
p2 = plot(X[1,:], X[2,:], Ŷ[2,:], lw = 3,
     title = "Neural Network Fit of U2(t)", color = c1,
     label = "Neural Network", xaxis = "x", yaxis="y",
     titlefont = "Helvetica", legendfont = "Helvetica",
     legend = :bottomright)
plot!(X[1,:], X[2,:], Ȳ[2,:], lw = 3, label = "True Missing Term", color=c2)

p3 = scatter(t, Xₙ[1,:], color = [c1], label = "x data",
             title = "Extrapolated Fit From Incomplete Training Data",
             titlefont = "Helvetica", legendfont = "Helvetica",
             markersize = 5)

scatter!(ty, [YS[:, 1]; YS[end, end]], color = c2, label = "y data", markersize = 5)
plot!(p3,true_solution_long, color = [c1 c2], linestyle = :dot, lw=5, label = ["True x(t)" "True y(t)"])
plot!(p3,estimate_long, color = [c3 c4], lw=2, label = ["Estimated x(t)" "Estimated y(t)"])
plot!(p3,[5.99,6.01],[0.0,maximum(hcat(Array(solution),Array(estimate_long)))],lw=1,color=:black, label = nothing)
annotate!([(3,12,text("Training \nData", 10, :center, :top, :black, "Helvetica"))])
l = @layout [grid(1,2)
             grid(1,1)]
plot(p1,p2,p3,layout = l)

savefig(joinpath(pwd(),"plots","$(svname)full_plot.pdf"))


text:
"""
    QuantExQASM.jl - implementation of circuit generator for
    different quantum algorithms.
"""
module QuantExQASM

include("GateOps.jl")
include("CircuitList.jl")
include("Circuit.jl")
include("algs/NCU.jl")
include("algs/Oracle.jl")
include("algs/Diffusion.jl")
include("algs/Grover.jl")
include("Algorithms.jl")

end


text:
typealias IntMatrix{w,T} IntArray{w,T,2}

function call{w,T}(::Type{IntMatrix{w,T}}, m::Integer, n::Integer, mmap::Bool=false)
    return IntArray{w,T}((m, n), mmap)
end

function call{w,T}(::Type{IntMatrix{w,T}}, mmap::Bool=false)
    return IntArray{w,T}((0, 0), mmap)
end

function convert{w,T}(::Type{IntMatrix{w}}, matrix::AbstractMatrix{T})
    return convert(IntArray{w,T,2}, matrix)
end


text:
"""
    Module for functions shared by different parameterizations.
"""
module Common

import Thermodynamics
import CLIMAParameters
import CLIMAParameters.Planet
import CLIMAParameters.Atmos.Microphysics

const TD = Thermodynamics
const CP = CLIMAParameters
const CP_planet = CLIMAParameters.Planet
const CP_micro = CLIMAParameters.Atmos.Microphysics
const APS = CP.AbstractParameterSet

export G_func

"""
    G_func(param_set, T, Liquid())
    G_func(param_set, T, Ice())

 - `param_set` - abstract set with earth parameters
 - `T` - air temperature
 - `Liquid()`, `Ice()` - liquid or ice phase to dispatch over.

Utility function combining thermal conductivity and vapor diffusivity effects.
"""
function G_func(param_set::APS, T::FT, ::TD.Liquid) where {FT <: Real}

    _K_therm::FT = CP_micro.K_therm(param_set)
    _R_v::FT = CP_planet.R_v(param_set)
    _D_vapor::FT = CP_micro.D_vapor(param_set)

    L = TD.latent_heat_vapor(param_set, T)
    p_vs = TD.saturation_vapor_pressure(param_set, T, TD.Liquid())

    return FT(1) / (
        L / _K_therm / T * (L / _R_v / T - FT(1)) + _R_v * T / _D_vapor / p_vs
    )
end
function G_func(param_set::APS, T::FT, ::TD.Ice) where {FT <: Real}

    _K_therm::FT = CP_micro.K_therm(param_set)
    _R_v::FT = CP_planet.R_v(param_set)
    _D_vapor::FT = CP_micro.D_vapor(param_set)

    L = TD.latent_heat_sublim(param_set, T)
    p_vs = TD.saturation_vapor_pressure(param_set, T, TD.Ice())

    return FT(1) / (
        L / _K_therm / T * (L / _R_v / T - FT(1)) + _R_v * T / _D_vapor / p_vs
    )
end

end


text:
#   Unit examples of the FourierAnalysis Package for julia language
#   v 0.0.1 - last update 24th of September 2019
#
#   MIT License
#   Copyright (c) 2019, Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This example shows how to compute univariate measures in
#   the time-frequency domain.

#   ~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~  #
#                                                                             #
#   ~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~  #

using FourierAnalysis, FFTW, LinearAlgebra, Statistics, Plots

# add module for reading the two EEG text files to be used ater
push!(LOAD_PATH, @__DIR__)
using IOtxt

# Get EEG file names with complete path (they have extension .txt)
S=getFilesInDir(@__DIR__; ext=(".txt",))

# read the two EEG data files and put them in a Matrix object
X1=readEEG(S[1])
X2=readEEG(S[2])

#############################################################################
# Univariate measures
# see Congedo, 2018: https://hal.archives-ouvertes.fr/hal-01868538/document

sr, wl, bandwidth=128, 512, 2
Pz=15
𝐱=[X1[:, Pz], X2[:, Pz]] # get the two times-series at electrode Pz
𝐘=TFanalyticsignal(𝐱, sr, wl, bandwidth; fmax=32, nonlinear=false)
𝐀=TFamplitude(𝐘)

# Mean Amplitude Eq. 0.1 (MAmp)
# compute the MAmp averaging in a TF region from a TFAnalyticSignalVector object
MAmp=meanAmplitude(𝐘, (8, 12), (1, 512); mode=mean)
# compute the MAmp averaging in a TF region from a TFAmplitudeVector object
MAmp=meanAmplitude(𝐀, (8, 12), (1, 512); mode=mean)
# compute the MAmp averaging in a TF region directly from data
MAmp=meanAmplitude(𝐱, sr, wl, (8, 12), (1, 512), bandwidth; mode=mean)

# compute the MAmp in a TF region from a TFAnalyticSignalVector object
MAmp=meanAmplitude(𝐘, (8, 12), (1, 512); mode=extract)
# compute the MAmp in a TF region from a TFAmplitudeVector object
MAmp=meanAmplitude(𝐀, (8, 12), (1, 512); mode=extract)
# compute the MAmp in a TF region directly from data
MAmp=meanAmplitude(𝐱, sr, wl, (8, 12), (1, 512), bandwidth; mode=extract)
# NB the Analytic Signal or the amplitude objects must be linear (note the nonlinear=false above)

# All these operation can be obtained on smoothed Amplitude, e.g.,
MAmp=meanAmplitude(𝐱, sr, wl, (8, 12), (1, 512), bandwidth;
                   mode=extract,
                   fsmoothing=hannSmoother,
                   tsmoothing=hannSmoother)

# The same operations can be done with the other univariate measures. For example:
# Concentration Eq. 0.2 (Con)
# compute the Con averaging in a TF region from a TFAnalyticSignalVector object
Con=concentration(𝐘, (8, 12), (1, 512); mode=mean)
# compute the Con in a TF region directly from data
Con=concentration(𝐱, sr, wl, (8, 12), (1, 512), bandwidth; mode=extract)
# Mean Direction Eq. 0.3 (MDir)
# compute the MDir averaging in a TF region directly from data
MDir=meanDirection(𝐱, sr, wl, (8, 12), (1, 512), bandwidth; mode=mean)
# compute the MDir in a TF region from a TFAnalyticSignalVector object
MDir=meanDirection(𝐘, (8, 12), (1, 512); mode=extract)

# and for the non-linear counterpart:
# Phase Concentration Eq. 0.4 (PCon)
# compute the Con in a TF region directly from data
Con=concentration(𝐱, sr, wl, (8, 12), (1, 512), bandwidth; mode=extract, nonlinear=true)
# Phase Mean Direction Eq. 0.5 (PMDir)
# compute the MDir averaging in a TF region directly from data
MDir=meanDirection(𝐱, sr, wl, (8, 12), (1, 512), bandwidth; mode=mean, nonlinear=true)

# If you try to compute a non-linear measure from a linear AnalyticSignal object
# you will get en error (see the REPL)
Con=concentration(𝐘, (8, 12), (1, 512); mode=mean, nonlinear=true)

# In order to compute those quantities from Analytic Signal objects
# first we need to compute a non-linear Analytic Signal objects:
𝐘=TFanalyticsignal(𝐱, sr, wl, bandwidth; fmax=32, nonlinear=true)
# then
# Phase Concentration Eq. 0.4 (PCon)
Con=concentration(𝐘, (8, 12), (1, 512); mode=mean, nonlinear=true)
# compute the MDir in a TF region from a TFAnalyticSignalVector object
MDir=meanDirection(𝐘, (8, 12), (1, 512); mode=extract, nonlinear=true)
# NB the Analytic Signal or the amplitude objects must be non-linear (note the nonlinear=true above)

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


text:
###################################################################
############## Structures Common to All equipment #################
###################################################################
#an object that contains all cost factors used in the calculations
mutable struct cost_ks
   FC_ac::Float32
   FC_dc::Float32
   dc::Float32
   f_ct::Float32
   p_ct::Float32
   c_ct::Float32
   Qc_oss::Float32
   Qc_pcc::Float32
   life::Float32
   T_op::Float32
   E_op::Float32
   cf::Float32
   FC_bld::Float32
   p2e::Float32
   ldc::Float32
   lac::Float32
   lmac::Float32
   npv::Float32
   opx_c::Float32
   opx_pl::Float32
   opx_x::Float32
   opx_co::Float32
   opx_cp::Float32
   conv_c::Float32#incremental value of converter cost
   conv_d::Float32#fixed value of converter cost
   pac_e::Float32#incremental value of ac platform cost
   pac_f::Float32#fixed value of ac platform cost
   pdc_g::Float32#incremental value of dc platform cost
   pdc_h::Float32#fixed value of dc platform cost
end
cost_ks()=cost_ks(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0)

#Contains all cable data sourced from www.abb.com/cables
mutable struct cbls_data
   cbls33kV::Array{Array{Float32, 1}, 1}
   cbls66kV::Array{Array{Float32, 1}, 1}
   cbls132kV::Array{Array{Float32, 1}, 1}
   cbls220kV::Array{Array{Float32, 1}, 1}
   cbls400kV::Array{Array{Float32, 1}, 1}
   cbls150kV::Array{Array{Float32, 1}, 1}
   cbls300kV::Array{Array{Float32, 1}, 1}
end
cbls_data()=cbls_data(Array{Array{Float32, 1}, 1}(),Array{Array{Float32, 1}, 1}(),Array{Array{Float32, 1}, 1}(),Array{Array{Float32, 1}, 1}(),Array{Array{Float32, 1}, 1}(),Array{Array{Float32, 1}, 1}(),Array{Array{Float32, 1}, 1}())

#xy co-ordinate struct
#May move but an element of both transformer and cables
mutable struct xy
      x::Float64
      y::Float64
end
xy()=xy(69.69,69.69)

#gps struct
#May move but an element of both transformer and cables
mutable struct gps
      lat::Float64
      lng::Float64
end
gps()=gps(69.69,69.69)

#location and number of node struct
#May move but an element of both transformer and cables
mutable struct node
      gps::gps
      xy::xy
      num::Int32
end
node()=node(gps(),xy(),69)

#wind object
#May move but an element of both transformer and cables
mutable struct wind
      pu::Array{Float32}
      ce::Array{Float32}
      delta::Float32
      lf::Float32
end
wind()=wind(Float32[],Float32[],69.69,69.69)

#reliability structure - same for all equipment
mutable struct reliability
   fr::Float32#failure rate
   mttr::Float32#mean time to repair
   mc::Float32#mean cost of repair
end
reliability()=reliability(69.69,69.69,69.69)

###################################################################
################# Structures for Platforms ########################
###################################################################
#the structure of costs for an OSS platform
mutable struct costs_plat
   cpx::Float32#MEuros procurement and installation
   cm::Float32
   ttl::Float32
end
costs_plat()=costs_plat(0.0,0.0,0.0)

#the structure for a platform
mutable struct platform
   acdc::String#system mva
   kv::Float32#system mva
   mva::Float32#system mva
   depth::Float32#water depth
   costs::costs_plat
   wnd::String
end
platform()=platform("acdc",0.0,0.0,0.0,costs_plat(),"")

###################################################################
################# Structures for Cables ###########################
###################################################################

#the structure of costs for a cable
mutable struct costs_cbl
   perkm_cpx::Float32#MEuros per km-procurement capex only
   pf3d::Float32#MEuros per km-procurement capex only
   perkm_ttl::Float32#MEuros per km-ttl cost
   qc::Float32#MEuros compensation on OSS
   sg::Float32#MEuros switch gear on OSS
   cpx_p::Float32#MEuros procurement capex total
   cpx_i::Float32#MEuros installation capex total
   rlc::Float32#MEuros losses in cable total
   cm::Float32#MEuros cost of corrective maintenance over lifetime
   eens::Float32#MEuros cost of eens over lifetime
   ttl::Float32#MEuros total lifetime cost
   grand_ttl::Float32#includes platform and double cable if mid point compensation
end
costs_cbl()=costs_cbl(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0)

#electrical properties of a cable
mutable struct elec_cbl
   mva::Float32#per core capacity
   amp::Float32
   volt::Float32
   ohm::Float32
   farrad::Float32
   henry::Float32
   yc::Float32
   xl::Float32
   freq::Float32
end
elec_cbl()=elec_cbl(69.69,69.69,69.69,69.69,69.69,69.69,69.69,69.69,50)

#the structure used for a cable
mutable struct cable
   mva::Float32#system requirement
   length::Float32
   path::Array{node}#nodes along path in direction of upstream(nth owpp) to downstream(pcc)
   size::Float32#cross setion of core
   num::Float32#number of parallel conductors
   relia::reliability
   elec::elec_cbl
   costs::costs_cbl
   wnd::String
   reactors::Array{Int32}
   plat::platform
   mpc_ac::Bool
   mx_rng::Float32
end
cable()=cable(0.0,0.0,node[],0.0,0.0,reliability(),elec_cbl(),costs_cbl(),"",Int32[],platform(),false,0.0)

###################################################################
################# Structures for transformers #####################
###################################################################
#the structure of costs for a transformers
mutable struct costs_xfo
   cpx_p::Float32#MEuros procurement capex total
   cpx_i::Float32#MEuros installation capex total
   tlc::Float32
   cm::Float32
   eens::Float32
   ttl::Float32
end
costs_xfo()=costs_xfo(0.0,0.0,0.0,0.0,0.0,0.0)
#electrical properties of a transformer
mutable struct elec_xfo
   mva::Float32#per transformer capacity
   amp::Float32
   hv::Float32
   lv::Float32
   ohm::Float32
   farrad::Float32
   henry::Float32
   yc::Float32
   xl::Float32
   freq::Float32
end
elec_xfo()=elec_xfo(69.69,69.69,69.69,69.69,69.69,69.69,69.69,69.69,69.69,50)

#the structure for a transformer
mutable struct transformer
   mva::Float32#system mva
   num::Float32#in parallel
   eta::Float32#efficiency
   relia::reliability
   elec::elec_xfo
   costs::costs_xfo
   wnd::String
end
transformer()=transformer(0.0,0.0,0.994,reliability(),elec_xfo(),costs_xfo(),"")

###################################################################
############################# HVDC Station ########################
###################################################################
#the structure of costs for an HVDC converter
mutable struct costs_hvdc
   cpx::Float32#MEuros procurement and installation
   tlc::Float32
   eens::Float32
   cm::Float32
   ttl::Float32
end
costs_hvdc()=costs_hvdc(0.0,0.0,0.0,0.0,0.0)

#Not used other than mva but may need for testing meshing inthe future
mutable struct elec_conv
   mva::Float32#per converter capacity
   amp::Float32
   hv::Float32
   lv::Float32
   ohm::Float32
   farrad::Float32
   henry::Float32
   yc::Float32
   xl::Float32
   freq::Float32
end
elec_conv()=elec_conv(69.69,69.69,69.69,69.69,69.69,69.69,69.69,69.69,69.69,50)

#the structure for an HVDC converter
mutable struct converter
   mva::Float32#system mva
   eta::Float32#efficiency
   num::Float32#in parallel - always 1
   relia::reliability
   elec::elec_conv
   costs::costs_hvdc
   wnd::String
end
converter()=converter(0.0,0.9828,1,reliability(),elec_conv(),costs_hvdc(),"")


text:
function vertical_cut(A::FuzzyNumber; α::Real=0.2)
    B = copy(A)
    vertical_cut!(B, α=α)
    B
end

function vertical_cut!(A::FuzzyNumber; α::Real=0.2)
	_lvl = Int(floor(length(A.levels) / (1 / α)))
	for lvl = 1:_lvl
		A.grades[lvl] = A.grades[_lvl + 1]
	end
end

function vertical_cut(A::FuzzyNumber, left::Real, right::Real)
    B = copy(A)
	for lvl = 1:length(B.levels)
        _left = B.grades[lvl].left
        _right = B.grades[lvl].right
        if (!isnothing(left) && _left < left) _left = left end
        if (!isnothing(right) && _right < right) _right = right end
		B.grades[lvl] = Interval(_left, _right)
	end
	B
end

function linearly_dampen(A::FuzzyNumber; multiplier::Real=0.5)
    if isSingleton(A)
        return A
    end
    B = copy(A)
    FuzzySets.linearly_dampen!(B, multiplier=multiplier)
    B
end

function linearly_dampen!(A::FuzzyNumber; multiplier::Real=0.5)
    if isSingleton(A)
        return
    end
    x3 = support(A).right
    x2 = peak(A)
    x1 = support(A).left
    y3 = A(x3)
    y2 = A(x2)
    y1 = A(x1)
    bound_0 = (x2 - x1 == 0) ? true : false
    bound_1 = (x3 - x2 == 0) ? true : false
    left_slope = (y2 - y1) / (x2 - x1)
    right_slope = (y3 - y2) / (x3 - x2)
    left_slope /= multiplier
    right_slope /= multiplier
    left_intercept = y2 - (left_slope * x2)
    right_intercept = y2 - (right_slope * x2)
    for (lvl, α) in enumerate(A.levels)
        # x = (y - b) / m
        left = (α - left_intercept) / left_slope
        right = (α - right_intercept) / right_slope
        left = bound_0 ? 0 : left
        right = bound_1 ? 1 : right
        left = min(left, x2)
        right = max(x2, right)
		A.grades[lvl] = Interval(left, right)
	end
end

function reflection_dampen(A::FuzzyNumber)
    B = copy(A)
    reflection_dampen!(B)
    B
end

function reflection_dampen!(A::FuzzyNumber)
    if isSingleton(A)
        return
    end
    x3 = support(A).right
    x2 = peak(A)
    x1 = support(A).left

    bound_0 = (x2 - x1 == 0) ? true : false
    bound_1 = (x3 - x2 == 0) ? true : false

    is_left_narrower =  x2 - x1 < x3 - x2
    for lvl = 1:length(A.levels)
        if is_left_narrower
            left = A.grades[lvl].left
            right = x2 + (x2 - left)
        else
            right = A.grades[lvl].right
            left = x2 - (right - x2)
        end
        # left = bound_0 ? 0 : left
        # right = bound_1 ? 1 : right
        left = min(left, x2)
        right = max(x2, right)
		A.grades[lvl] = Interval(left, right)
	end
end

function midpoint_dampen(A::FuzzyNumber)
    B = copy(A)
    midpoint_dampen!(B)
    B
end

function midpoint_dampen!(A::FuzzyNumber; δ::Float64=0.5)
    if isSingleton(A)
        return
    end

    peak_mid = peak(A)
    for lvl = 1:length(A.levels)
        cut = A.grades[lvl]
        midpoint = mid(cut)
        # is_left_narrower =  midpoint - cut.left < cut.right - midpoint
        # if is_left_narrower
        #     midpoint =
        # else
        # end
        Δleft = midpoint - cut.left
        Δright = cut.right - midpoint
        left = peak_mid - (Δleft * δ)
        right = peak_mid + (Δright * δ)
        left = max(min(left, peak_mid), 0)
        right = min(max(peak_mid, right), 1)
		A.grades[lvl] = Interval(left, right)
	end
end


text:
module Hopf

using ..HarmonicBalance
using Symbolics

include("Hopf/Hopf.jl")

end


text:
using DrWatson
@quickactivate
using HierarchicalAD
using Flux, CUDA
using BSON, FileIO
using ValueHistories, DataFrames
using Plots
using ArgParse
using IterTools

using HierarchicalAD: construct_conv_classifier, train_conv_classifier!, create_input_data
using HierarchicalAD: auc_val, classifier_score

# get the data
X, y = HierarchicalAD.load_shapes2D();
datasize = size(X)

# set device id
CUDA.device!(2)

# setup paths
outpath = datadir("baseline_comparison/conv_simple")
mkpath(outpath)

# original defaults
nns = floor.(Int, 10.0 .^ range(1, 5, length=10))
nas = floor.(Int, 10.0 .^ range(1, 4, length=4))
ntst = 25000
label_vec = [
	y[!,:shape] .== 2, y[!,:shape] .!= 2,
	y[!,:scale] .< 0.75, y[!,:scale] .> 0.75,
	y[!,:posX] .> 0.5, y[!,:posY] .> 0.5,
	y[!,:normalized_orientation] .< 0.1,
	.!((y[!,:shape] .== 1) .| ((y[!,:posX] .> 0.5) .& (y[!,:posY] .> 0.5)))
	]
label_desc = [
	"ovals normal", "ovals anomalous",
	"small normal", "large normal",
	"right normal", "bottom normal",
	"rotated anomalous",
	"squares or right bottom anomalous"
	]

# ovals normal
ilabel = 1
labels = label_vec[ilabel]
desc = label_desc[ilabel]
nns = floor.(Int, 10.0 .^ range(1, 4, length=8))
nns = range(400, 3000, step=200)
nas = range(100, 1000, step = 100)

# ovals anomalous
ilabel = 2
labels = label_vec[ilabel]
desc = label_desc[ilabel]
nns = floor.(Int, 10.0 .^ range(1, 4, length=8))
nns = range(400, 3000, step=200)
nas = range(100, 1000, step = 100)

# ovals anomalous
ilabel = 8
labels = label_vec[ilabel]
desc = label_desc[ilabel]
nns = floor.(Int, 10.0 .^ range(1, 4, length=8))
nns = range(400, 3000, step=200)
nas = range(100, 1000, step = 100)

conv_classifier_params = (lr = 0.001f0, λ = 0f0, batchsize = 64, nepochs = 1000, patience = 100, verb = true)

function train_and_evaluate_classifier()
	for seed = 1:1
		aucs_per_samples = []
		scores_per_samples = []
		tst_ys_per_samples = []
		for (nn, na) in product(nns, nas)
			@info "Training classifier using $(desc), seed=$seed."
			@info "# normal training/validation = $nn, # anomalous training/validation = $na, # test = $(ntst*2)"
			try
				global tr_nX, (tr_X, tr_y), (val_X, val_y), (tst_X, tst_y) = create_input_data(X, y, labels, nn, nn, na, na,
				    ntst; seed=seed);
			catch e
				if isa(e, BoundsError)
					@info "Not enough data, skipping training..."
					push!(aucs_per_samples, NaN)
					push!(scores_per_samples, [NaN])
					push!(tst_ys_per_samples, [NaN])
					continue
				else
					rethrow(e)
				end
			end
			# construct the classifier
			classifier_conv_part = [Conv((5,5), 1=>32, relu), MaxPool((3,3))]
			odims = outdims(classifier_conv_part, datasize)
			odim = reduce(*, odims[1:3])
			classifier = gpu(Chain(classifier_conv_part..., x->reshape(x, odim, size(x,4)),  Dense(odim, 2)))

			# train it
			classifier, history, opt = train_conv_classifier!(classifier, tr_X, tr_y, val_X, val_y; conv_classifier_params...);
			scores = cpu(classifier_score(classifier, tst_X, 128))
			auc = auc_val(tst_y, scores)
			@info "Finished training, test AUC = $(auc)"
			push!(aucs_per_samples, auc)
			push!(scores_per_samples, scores)
			push!(tst_ys_per_samples, tst_y)
		end
		df = replace(desc, " " => "_")
		savef = joinpath(outpath, "classifier_output_$(df)_seed=$(seed).bson")
		outdata = Dict(
			:nns => nns, :nas => nas, :ntst => ntst,
			:conv_classifier_params => conv_classifier_params,
			:aucs => aucs_per_samples,
			:scores => scores_per_samples,
			:tst_ys => tst_ys_per_samples,
			:desc => desc
			)
		save(savef, outdata)
		@info "Saved to $savef"
	end
end

outpath = datadir("baseline_comparison/conv_simple")
mkpath(outpath)

nns = vcat(range(30, 370, step=30), range(400, 3000, step=200))
nas = vcat([10], collect(range(100, 1000, step = 200)))

for ilabel in 1:8
	labels = label_vec[ilabel]
	desc = label_desc[ilabel]
	train_and_evaluate_classifier()
end

train_and_evaluate_classifier()

nns = range(400, 3000, step=200)
nas = range(100, 1000, step = 100)

#
ilabel = 8
labels = label_vec[ilabel]
desc = label_desc[ilabel]

train_and_evaluate_classifier()

#
ilabel = 5
labels = label_vec[ilabel]
desc = label_desc[ilabel]

train_and_evaluate_classifier()

#
ilabel = 3
labels = label_vec[ilabel]
desc = label_desc[ilabel]

train_and_evaluate_classifier()

########
outpath = datadir("baseline_comparison/conv_simple_smaller_ns")
mkpath(outpath)

ilabel = 1
labels = label_vec[ilabel]
desc = label_desc[ilabel]
nns = range(10, 400, step=30)
nas = vcat([10], collect(range(100, 1000, step = 200)))
train_and_evaluate_classifier()

# the most compelte run
outpath = datadir("baseline_comparison/had_simple")
mkpath(outpath)

nns = vcat(range(10, 370, step=30), range(400, 3000, step=200))
nas = vcat([10], collect(range(100, 1000, step = 200)))

mfs = vcat(["activation=relu_channels=8-8-16-32-64_data=shapes2D_gamma=50.0_lambda=0.0_latent_dim=4_model=fvlae_model_id=20210817174044364_xdist=bernoulli.bson"],
	readdir(modeldir))
for mf in mfs
	for ilabel in 1:8
		labels = label_vec[ilabel]
		desc = label_desc[ilabel]
		autoencoder_data = load(joinpath(modeldir, mf))
		train_and_evaluate_had(mf, autoencoder_data)
	end
end

max_stars_repo_name: vitskvara/HierarchicalAD

max_issues_repo_name: vitskvara/HierarchicalAD

max_forks_repo_name: vitskvara/HierarchicalAD


text:
module NewtonsMethod
export newtonroot

using ForwardDiff
# Write your package code here.

function newtonroot(f, f_prime; x_0, tol = 1e-7, maxiter = 1000)
    err = abs(f(x_0))
    iter = 0
    x = x_0
    while err > tol && iter <= maxiter
        if f_prime(x) == 0 || f_prime(x) == Inf || f_prime(x) == -Inf
            if rand() >= 0.5
                x = x + 1e-10
            else
                x = x - 1e-10
            end
            iter = iter+1
            err = abs(f(x))
            continue
        end
        iter = iter+1
        xNew = x - f(x)/f_prime(x)
        err = abs(f(x))
        x = xNew
    end
    if iter > maxiter
        return nothing
    end
    return x
end

function newtonroot(f; x_0, tol = 1e-7, maxiter = 1000)
    f_prime = x -> ForwardDiff.derivative(f, x)
    return newtonroot(f, f_prime; x_0 = x_0, tol = tol, maxiter = maxiter)
end
end


text:
for s in enumerate((:vertextype, :facetype, :normaltype,
                       :texturecoordinatetype, :colortype))
    @eval begin
        $(s[2]){T<:HomogenousMesh}(t::Type{T}) = t.parameters[$(s[1])]
        $(s[2])(mesh::HomogenousMesh) = $(s[2])(typeof(mesh))
    end
end

hasvertices(msh) = vertextype(msh) != Void
hasfaces(msh) = facetype(msh) != Void
hasnormals(msh) = normaltype(msh) != Void
hastexturecoordinates(msh) = texturecoordinatetype(msh) != Void
hascolors(msh) = colortype(msh) != Void

vertices(msh) = msh.vertices
faces(msh) = msh.faces
normals(msh) = msh.normals
texturecoordinates(msh) = msh.texturecoordinates
colors(msh) = msh.color

# Bad, bad name! But it's a little tricky to filter out faces and verts from the attributes, after get_attribute
function attributes_noVF{T<:AbstractMesh}(m::T)
    fielddict = Dict{Symbol, Any}(map(fieldnames(T)[3:end]) do field
        field => getfield(m, field)
    end)
    return filter(fielddict) do key,val
        val != nothing && val != Void[]
    end
end
#Gets all non Void attributes from a mesh in form of a Dict fieldname => value
function attributes(m::AbstractMesh)
    filter((key,val) -> (val != nothing && val != Void[]), all_attributes(m))
end
#Gets all non Void attributes types from a mesh type fieldname => ValueType
function attributes{M <: HMesh}(m::Type{M})
    filter((key,val) -> (val != Void && val != Vector{Void}), all_attributes(M))
end

function all_attributes{M <: HMesh}(m::Type{M})
    Dict{Symbol, Any}(map(field -> (field => fieldtype(M, field)), fieldnames(M)))
end
function all_attributes{M <: HMesh}(m::M)
    Dict{Symbol, Any}(map(field -> (field => getfield(m, field)),  fieldnames(M)))
end

# Needed to not get into an stack overflow
convert{M <: AbstractMesh}(::Type{M}, mesh::AbstractGeometry) = M(mesh)
@compat (::Type{HM1}){HM1 <: AbstractMesh}(mesh::HM1) = mesh

"""
Uses decompose to get all the converted attributes from the meshtype and
creates a new mesh with the desired attributes from the converted attributs
Getindex can be defined for any arbitrary geometric type or exotic mesh type.
This way, we can make sure, that you can convert most of the meshes from one type to the other
with minimal code.
"""
@compat function (::Type{HM1}){HM1 <: AbstractMesh}(primitive::Union{AbstractMesh, GeometryPrimitive})
    result = Dict{Symbol, Any}()
    for (field, target_type) in zip(fieldnames(HM1), HM1.parameters)
        if target_type != Void
            if field == :attribute_id
                if !isa(primitive, HomogenousMesh)
                    error("Primitive $primitive doesn't hold attribute indexes")
                end
                result[field] = primitive.attribute_id
            else
                result[field] = decompose(target_type, primitive)
            end
        end
    end
    HM1(result)
end

isvoid{T}(::Type{T}) = false
isvoid(::Type{Void}) = true
isvoid{T}(::Type{Vector{T}}) = isvoid(T)
@compat function (::Type{HM1}){HM1 <: HomogenousMesh}(primitive::HomogenousMesh)
    fnames = fieldnames(HM1)
    args = ntuple(nfields(HM1)) do i
        field, target_type = fnames[i], fieldtype(HM1, i)
        soure_type = fieldtype(typeof(primitive), i)
        isa(HM1.parameters[i], TypeVar) && return getfield(primitive, field) # target is not defined
        if !isvoid(target_type) && isvoid(soure_type) # target not there yet, maybe we can decompose though (e.g. normals)
            return decompose(HM1.parameters[i], primitive)
        elseif isvoid(target_type)
            return target_type()
        else
            return convert(target_type, getfield(primitive, field))
        end
    end
    HM1(args...)
end

#Should be:
#function call{M <: HMesh, VT <: Point, FT <: Face}(::Type{M}, vertices::Vector{VT}, faces::Vector{FT})
# Haven't gotten around to implement the types correctly with abstract types in FixedSizeArrays
@compat function (::Type{M}){M <: HMesh, VT, FT <: Face}(
        vertices::Vector{Point{3, VT}}, faces::Vector{FT}
    )
    msh = PlainMesh{VT, FT}(vertices=vertices, faces=faces)
    convert(M, msh)
end
get_default(x::Union{Type, TypeVar}) = nothing
get_default{X <: Array}(x::Type{X}) = Void[]

"""
generic constructor for abstract HomogenousMesh, infering the types from the keywords (which have to match the field names)
some problems with the dispatch forced me to use this method name... need to further investigate this
"""
function homogenousmesh(attribs::Dict{Symbol, Any})
    newfields = []
    for name in fieldnames(HMesh)
        push!(newfields, get(attribs, name, get_default(fieldtype(HMesh, name))))
    end
    HomogenousMesh(newfields...)
end

"""
Creates a mesh from keyword arguments, which have to match the field types of the given concrete mesh
"""
@compat (::Type{M}){M <: HMesh}(; kw_args...) = M(Dict{Symbol, Any}(kw_args))

"""
Creates a new mesh from a dict of `fieldname => value` and converts the types to the given meshtype
"""
@compat function (::Type{M}){M <: HMesh}(attribs::Dict{Symbol, Any})
    newfields = map(zip(fieldnames(HomogenousMesh), M.parameters)) do field_target_type
        field, target_type = field_target_type
        default = fieldtype(HomogenousMesh, field) <: Vector ? Void[] : nothing
        get(attribs, field, default)
    end
    M(HomogenousMesh(newfields...))
end

"""
Creates a new mesh from an old one, with changed attributes given by the keyword arguments
"""
@compat function (::Type{M}){M <: HMesh}(mesh::AbstractMesh, attributes::Dict{Symbol, Any})
    newfields = map(fieldnames(HomogenousMesh)) do field
        get(attributes, field, getfield(mesh, field))
    end
    HomogenousMesh(newfields...)
end
"""
Creates a new mesh from an old one, with a new constant attribute (like a color)
"""
@compat function (::Type{HM}){HM <: HMesh, ConstAttrib}(mesh::AbstractMesh, constattrib::ConstAttrib)
    result = Dict{Symbol, Any}()
    for (field, target_type) in zip(fieldnames(HM), HM.parameters)
        if target_type <: ConstAttrib
            result[field] = constattrib
        elseif target_type != Void
            result[field] = decompose(target_type, mesh)
        end
    end
    HM(result)
end
function add_attribute(m::AbstractMesh, attribute)
    attribs = attributes(m) # get all attribute values as a Dict fieldname => value
    attribs[:color] = attribute # color will probably be renamed to attribute. not sure yet...
    homogenousmesh(attribs)
end

"""
Creates a new mesh from a tuple of a geometry type and a constant attribute
"""
@compat function (::Type{HM}){HM <: HMesh, ConstAttrib, P<:AbstractGeometry}(x::Tuple{P, ConstAttrib})
    any, const_attribute = x
    add_attribute(HM(any), const_attribute)
end

merge{M <: AbstractMesh}(m::Vector{M}) = merge(m...)

"""
Merges an arbitrary mesh. This function probably doesn't work for all types of meshes
parameters:
`m1` first mesh
`meshes...` other meshes
"""
function merge{M <: AbstractMesh}(m1::M, meshes::M...)
    v = copy(m1.vertices)
    f = copy(m1.faces)
    attribs = deepcopy(attributes_noVF(m1))
    for mesh in meshes
        append!(f, mesh.faces + length(v))
        append!(v, mesh.vertices)
        for (v1, v2) in zip(values(attribs), values(attributes_noVF(mesh)))
            append!(v1, v2)
        end
    end
    attribs[:vertices]  = v
    attribs[:faces]     = f
    return M(attribs)
end

"""
A mesh with one constant attribute can be merged as an attribute mesh.
Possible attributes are FSArrays
"""
function merge{_1, _2, _3, _4, ConstAttrib <: Colorant, _5, _6}(
        m1::HMesh{_1, _2, _3, _4, ConstAttrib, _5, _6},
        meshes::HMesh{_1, _2, _3, _4, ConstAttrib, _5, _6}...
    )
    vertices     = copy(m1.vertices)
    faces        = copy(m1.faces)
    attribs      = Dict{Symbol, Any}(filter((k,v) -> k != :color, attributes_noVF(m1)))
    attribs      = Dict{Symbol, Any}([(k, copy(v)) for (k,v) in attribs])
    color_attrib = RGBA{U8}[RGBA{U8}(m1.color)]
    index        = Float32[length(color_attrib)-1 for i=1:length(m1.vertices)]
    for mesh in meshes
        append!(faces, mesh.faces + length(vertices))
        append!(vertices, mesh.vertices)
        attribsb = attributes_noVF(mesh)
        for (k,v) in attribsb
            k != :color && append!(attribs[k], copy(v))
        end
        push!(color_attrib, mesh.color)
        append!(index, Float32[length(color_attrib)-1 for i=1:length(mesh.vertices)])
    end
    delete!(attribs, :color)
    attribs[:vertices]      = vertices
    attribs[:faces]         = faces
    attribs[:attributes]    = color_attrib
    attribs[:attribute_id]  = index
    return HMesh{_1, _2, _3, _4, Void, typeof(color_attrib), eltype(index)}(attribs)
end

# Fast but slightly ugly way to implement mesh multiplication
# This should probably go into FixedSizeArrays.jl, Vector{FSA} * FSA
immutable MeshMulFunctor{T}
    matrix::Mat{4,4,T}
end
@compat (m::MeshMulFunctor{T}){T}(vert) = Vec{3, T}(m.matrix*Vec{4, T}(vert..., 1))
function *{T}(m::Mat{4,4,T}, mesh::AbstractMesh)
    msh = deepcopy(mesh)
    map!(MeshMulFunctor(m), msh.vertices)
    msh
end


text:
using GeneralizedGenerated
using NestedTuples

abstract type Hypercube{k} end

struct RandHypercube{T} <: Hypercube{Inf}
    rng :: T
end

function Base.rand(ω::RandHypercube)
    rand(ω.rng)
end

# TODO: Add a test
function NestedTuples.with(m::Module, hcube_nt::NamedTuple{N,Tuple{H}}, tv::TupleVector{T,X}, ex::TypelevelExpr{E}) where {T,X,E, N, H<:Hypercube}
    n = length(tv)
    (ω,) = keys(hcube_nt)
    next!ω_expr = NestedTuples.TypelevelExpr(:($next!($ω)))
    next! = TupleVectors.next!
    NestedTuples.with(m, hcube_nt, next!ω_expr)
    result = @inbounds chainvec(NestedTuples.with(m, hcube_nt, tv[1], ex), n)
    for j in 2:n
        NestedTuples.with(m, hcube_nt, next!ω_expr)
        @inbounds result[j] = NestedTuples.with(m, hcube_nt, tv[j], ex)
    end
    return result
end

# TODO: Add a test
function NestedTuples.with(m::Module, hcube_nt::NamedTuple{N,Tuple{H}}, n :: Int, ex::TypelevelExpr{E}) where {T,X,E, N, H<:Hypercube}
    (ω,) = keys(hcube_nt)
    next!ω_expr = NestedTuples.TypelevelExpr(:($next!($ω)))
    NestedTuples.with(m, hcube_nt, next!ω_expr)
    result = @inbounds chainvec(NestedTuples.with(m, hcube_nt, ex), n)
    for j in 2:n
        NestedTuples.with(m, hcube_nt, next!ω_expr)
        @inbounds result[j] = NestedTuples.with(m, hcube_nt, ex)
    end
    return result
end

# function NestedTuples.with(ω::Hypercube, nt::NamedTuple, ex::TypelevelExpr{E}) where {E}
#     next!(ω)
#     NestedTuples.with(nt, ex)
# end

# TODO: Add interface for Hypercube

# using Colors
# using GLMakie

# function makeplot()
#     ω = SobolHypercube(100)
#     x = range(-2,2,length=100)
#     d = For(x) do xj Normal(2*xj, 1/(1 + xj^2)) end
#     y = rand(ω, d)

#     plt = scatter(x,y,markersize=2, color=colorant"rgba(0,0,0,0.1)")
#     for j in 1:100
#         next!(ω)
#         y = rand(ω, d)
#         scatter!(x, y, markersize=2, color=colorant"rgba(0,0,0,0.1)")
#     end
#     plt
# end

# makeplot()


text:
using Test
using PlanktonIndividuals

@testset "PlanktonIndividuals 3D tests:" begin
    include(dirname(pathof(PlanktonIndividuals))*"/../test/model_test_3D.jl")
    @test isapprox(TC,TCt; atol=1e2)
end #@testset "PlanktonIndividuals 3D tests:" begin

@testset "PlanktonIndividuals 2D tests:" begin
    include(dirname(pathof(PlanktonIndividuals))*"/../test/model_test_2D.jl")
    @test isapprox(TP,TPt; atol=1e1)
end #@testset "PlanktonIndividuals 2D tests:" begin

@testset "PlanktonIndividuals 1D tests:" begin
    include(dirname(pathof(PlanktonIndividuals))*"/../test/model_test_1D.jl")
    @test isapprox(TP,TPt; atol=1e1)
end #@testset "PlanktonIndividuals 1D tests:" begin

@testset "PlanktonIndividuals 0D tests:" begin
    include(dirname(pathof(PlanktonIndividuals))*"/../test/model_test_0D.jl")
    @test isapprox(TP,TPt; atol=1e1)
end #@testset "PlanktonIndividuals 0D tests:" begin


text:
using PyCall
@pyimport scipy.linalg.interpolative
println("Julia loaded")

function main(rows, cols, values)
    println("In julia")
    A = sparse(rows, cols, values)
    println(interpolative.interp_decomp(A, 2))
    println("Exiting julia")
end

max_stars_repo_name: ligmg/ligmg

max_issues_repo_name: ligmg/ligmg

max_forks_repo_name: ligmg/ligmg


text:
module Day01

using ..ReTest, ..DocStringExtensions

import ..solve, ..PuzzleInfo, ..TEST_DATA_DIRECTORY

"""
    $(FUNCTIONNAME)(values)

Count the number of positive differences between successive elements in `values`.

# Examples
```jldoctest; setup = :(using AoC2021.Day01: $(FUNCTIONNAME))
julia> $(FUNCTIONNAME)([199, 200, 208, 210, 200, 207, 240, 269, 260, 263])
7
```
"""
function countincreases(values::AbstractVector)
    count(x -> x > 0, diff(values))
end

"""
    $(FUNCTIONNAME)(values, windowsize)

Stack sliding windows of size `windowsize` of a 1D array `values` in to a 2D array.

# Examples
```jldoctest; setup = :(using AoC2021.Day01: $(FUNCTIONNAME))
julia> $(FUNCTIONNAME)([199, 200, 208, 210, 200, 207, 240, 269, 260, 263], 3)
8×3 Matrix{Int64}:
 199  200  208
 200  208  210
 208  210  200
 210  200  207
 200  207  240
 207  240  269
 240  269  260
 269  260  263
```
"""
function slidingwindows(values::AbstractVector, windowsize::Integer)
    reduce(
        vcat, @view(values[i:i+windowsize-1])' for i in 1:length(values)-windowsize+1
    )
end

"""
    $(FUNCTIONNAME)(input, type)

Read array from input stream `input` with one value of type `type` per line.
"""
function readarraydata(input::IO, type::Type)
    [parse(type, value) for value in eachline(input)]
end

abstract type Day01PuzzleInfo <: PuzzleInfo end

"""
Day 1 part 1 puzzle information.

$(FIELDS)
"""
struct Part1PuzzleInfo <: Day01PuzzleInfo
    """Path to input data file."""
    datapath::String
end

"""
Day 1 part 2 puzzle information.

$(FIELDS)
"""
struct Part2PuzzleInfo <: Day01PuzzleInfo
    """Path to input data file."""
    datapath::String
    """Size of sliding window to use."""
    windowsize::Int
end

Part2PuzzleInfo(datapath) = Part2PuzzleInfo(datapath, 3)

"""
    $(FUNCTIONNAME)(depths, info)

Apply preprocessing to depths data in `data` for puzzle part with information `info`.
"""
preprocessdepths(depths::AbstractVector, ::Day01PuzzleInfo) = depths

function preprocessdepths(depths::AbstractVector, info::Part2PuzzleInfo)
    vec(sum(slidingwindows(depths, info.windowsize), dims=2))
end

function solve(info::Day01PuzzleInfo)
    depths = open(info.datapath, "r") do input
        readarraydata(input, Int)
    end
    processeddepths = preprocessdepths(depths, info)
    countincreases(processeddepths)
end

@testset begin
    @test solve(Part1PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day01.txt"))) == 7
    @test solve(Part2PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day01.txt"), 3)) == 5
end

end # module


text:
@generated function offsetof(::Type{T}, ::Val{_field}) where {T<:CBinding.Cstruct,_field}
    offset = -1
    for field in CBinding.fields(T).parameters
        if field.parameters[1] == _field
            offset = field.parameters[2].parameters[4]
            break
        end
    end
    @assert offset >= 0 "Failed to find offset of $_field in $T"
    :($offset)
end
@generated function offsetof(::Type{T}, ::Val{_field}) where {T,_field}
    @assert isstructtype(T) "offsetof is only valid for structs"
    offset = -1
    for (idx,field) in enumerate(fieldnames(T))
        if field == _field
            offset = fieldoffset(T, idx)
            break
        end
    end
    @assert offset >= 0 "Failed to find offset of $_field in $T"
    :($offset)
end
"Returns the integer offset of `field` in `T`."
macro offsetof(T, field)
    :($offsetof($T, Val($field)))
end

@generated _fieldtype(::Type{T}, ::Val{field}) where {T,field} =
    :($(fieldtype(T, field)))
@generated _fieldtype(::Type{T}, ::Val{field}) where {T<:CBinding.Cstruct,field} =
    :($(CBinding.field(T, field).parameters[2]))
@inline _fieldtype(::Type{Cptr{T}}, ::Val{field}) where {T,field} =
    _fieldtype(T, Val(field))


text:
# Simple test UDP server
function udp_srv(port::Int)
    r = Future()

    sock = UDPSocket()
    bind(sock, ip"127.0.0.1", port)

    @async begin
        put!(r, String(recv(sock)))
        close(sock)
    end

    return r
end

@testset "Ext - Syslogs" begin
    levels = copy(Memento._log_levels)
    levels["invalid"] = 100

    # Run a local UDP server for testing
    r = udp_srv(8080)

    # Create our DefaultHandler w/ the Syslog IO type
    handler = DefaultHandler(
        Syslogs.Syslog(ip"127.0.0.1", 8080; tcp=false),
        DefaultFormatter("{level}: {msg}")
    )

    logger = Logger(
        "Syslogs",
        Dict("Syslog" => handler),
        "debug",
        levels,
        DefaultRecord,
        true
    )

    # We just want to test that our glue code works with Syslogs.jl
    @suppress Memento.info(logger, "Hello World!")
    s = fetch(r)
    @test occursin("Hello World!", s)
end


text:
"""
nk_heq()

CI for nsoli and H-equation.
"""
function nk_heq()
    n = 32;
    u0 = zeros(n);
    FS = zeros(n);
    FPS = zeros(n, 20);
    FPJ = zeros(n, n);
    FPV = zeros(n,)
    c = 0.999;
    atol = 1.e-9;
    rtol = 1.e-9;
    hdata = heqinit(u0, c);
    dout =
        nsol(heqf!, u0, FS, FPJ, heqJ!; rtol = rtol, atol = atol, pdata = hdata, sham = 1)
    kout = nsoli(
        heqf!,
        u0,
        FS,
        FPS;
        pdata = hdata,
        rtol = rtol,
        atol = atol,
        lmaxit = -1,
        eta = 0.1,
        fixedeta = false,
    );
    kout2 = nsoli(
        heqf!,
        u0,
        FS,
        FPS;
        pdata = hdata,
        rtol = rtol,
        atol = atol,
        lmaxit = 2,
        eta = 0.01,
    );
    kout3 = nsoli( heqf!, u0, FS, FPV; pdata = hdata, rtol = rtol,
      atol = atol, lmaxit = 40, eta = 0.1,
      fixedeta=true, lsolver = "bicgstab")
    ksol = kout.solution;
    dsol = dout.solution;
    ksol2 = kout2.solution;
    ksol3 = kout3.solution
    soltest = norm(ksol - dsol, Inf) + norm(ksol - ksol2, Inf) +
              norm(ksol3 - dsol,Inf)
    solpass = (soltest < 1.e-7)
    solpass || println("solpass fails")
    kfpass = (sum(kout2.stats.ikfail) == 9)
    kfpass || println("kfpass fails")
    histdiff = (dout.history - kout.history[1:8]) ./ dout.history[1]
    histpass = (norm(histdiff, Inf) < 1.e-2)
    histpass || println("histpass fails")
    histdiffb = (kout.history - kout3.history) ./ kout.history[1]
    histpassb = (norm(histdiffb, Inf) < 1.e-2)
    histpassb || println("histpassb fails")
    nkhpass = solpass && kfpass && histpass && histpassb
    return nkhpass
end


text:
#!/usr/bin/env julia

import Pkg
import Pkg.Types: VersionSpec, VersionRange, VersionBound, semver_spec
import Base: thismajor, thisminor, thispatch, nextmajor, nextminor, nextpatch

const STDLIBS = [
    "Base64"
    "CRC32c"
    "Dates"
    "DelimitedFiles"
    "Distributed"
    "FileWatching"
    "Future"
    "InteractiveUtils"
    "Libdl"
    "LibGit2"
    "LinearAlgebra"
    "Logging"
    "Markdown"
    "Mmap"
    "Pkg"
    "Printf"
    "Profile"
    "Random"
    "REPL"
    "Serialization"
    "SHA"
    "SharedArrays"
    "Sockets"
    "SparseArrays"
    "Statistics"
    "SuiteSparse"
    "Test"
    "Unicode"
    "UUIDs"
]

function uuid(name::AbstractString)
    if name == "Pkg"
        return "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
    elseif name == "Statistics"
        return "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
    else
        string(Pkg.METADATA_compatible_uuid(String(name)))
    end
end

function uses(repo::AbstractString, lib::AbstractString)
    pattern = string(raw"\b(import|using)\s+((\w|\.)+\s*,\s*)*", lib, raw"\b")
    success(`git -C $repo grep -Eq $pattern -- '*.jl'`)
end

function semver(intervals)
    spec = String[]
    for ival in intervals
        if ival.upper == v"∞"
            push!(spec, "≥ $(thispatch(ival.lower))")
        else
            lo, hi = ival.lower, ival.upper
            if lo.major < hi.major
                push!(spec, "^$(lo.major).$(lo.minor).$(lo.patch)")
                for major = lo.major+1:hi.major-1
                    push!(spec, "~$major")
                end
                for minor = 0:hi.minor-1
                    push!(spec, "~$(hi.major).$minor")
                end
                for patch = 0:hi.patch-1
                    push!(spec, "=$(hi.major).$(hi.minor).$patch")
                end
            elseif lo.minor < hi.minor
                push!(spec, "~$(lo.major).$(lo.minor).$(lo.patch)")
                for minor = lo.minor+1:hi.minor-1
                    push!(spec, "~$(hi.major).$minor")
                end
                for patch = 0:hi.patch-1
                    push!(spec, "=$(hi.major).$(hi.minor).$patch")
                end
            else
                for patch = lo.patch:hi.patch-1
                    push!(spec, "=$(hi.major).$(hi.minor).$patch")
                end
            end
        end
    end
    return join(spec, ", ")
end

if !isempty(ARGS) && ARGS[1] == "-f"
    const force = true
    popfirst!(ARGS)
else
    const force = false
end
isempty(ARGS) && (push!(ARGS, pwd()))

for arg in ARGS
    dir = abspath(expanduser(arg))
    isdir(dir) ||
        error("$arg does not appear to be a package (not a directory)")

    name = basename(dir)
    if isempty(name)
        dir = dirname(dir)
        name = basename(dir)
    end
    endswith(name, ".jl") && (name = chop(name, tail=3))

    project_file = joinpath(dir, "Project.toml")
    !force && isfile(project_file) &&
        error("$arg already has a project file")

    require_file = joinpath(dir, "REQUIRE")
    isfile(require_file) ||
        error("$arg does not appear to be a package (no REQUIRE file)")

    project = Dict(
        "name" => name,
        "uuid" => uuid(name),
        "deps" => Dict{String,String}(),
        "compat" => Dict{String,String}(),
        "extras" => Dict{String,String}(),
    )

    test_require_file = joinpath(dir, "test", "REQUIRE")

    for (file, section) in ((require_file, "deps"),
                            (test_require_file, "extras"))
        isfile(file) || continue
        reqs = Pkg.Pkg2.Reqs.read(file)
        for req in reqs
            req isa Pkg.Pkg2.Reqs.Requirement || continue
            dep = String(req.package)
            if dep != "julia"
                project[section][dep] = uuid(dep)
            end
            if req.versions != Pkg.Pkg2.Pkg2Types.VersionSet()
                project["compat"][dep] = semver(req.versions.intervals)
            end
        end

        for stdlib in STDLIBS
            if uses(dirname(file), stdlib)
                project[section][stdlib] = uuid(stdlib)
            end
        end
    end

    if !isempty(project["extras"])
        project["targets"] = Dict("test" => collect(keys(project["extras"])))
        haskey(project["deps"], "Test") && delete!(project["deps"], "Test")
    end

    println(stderr, "Generating project file for $name: $project_file")
    open(project_file, "w") do io
        Pkg.TOML.print(io, project, sorted=true)
    end
    project = Pkg.Types.read_project(project_file)
    Pkg.Types.write_project(project, project_file)
end


text:
## HiC data + loading cooler code

include("cooler.jl")

function loadhicdata(loadfuns = [load_islet_hic])

    dfs = [f() for f in loadfuns]
    ns = names(first(dfs))
    for df in dfs[2:end]
        ns = intersect(ns, names(df))
    end
    mapreduce(df -> df[!, ns], vcat, dfs)
end

function load_islet_hic(file="islet_loops_TSTFF938730.bedpe.gz", projdir=getprojectdir())
    filepath = joinpath(projdir, "data", file)
    isletloops = CSV.read(filepath, DataFrame)
    rename!(isletloops, "#chrA" => "chrA")
    isletloops[!, :Study] .= "IsletHiC"
    isletloops
end


text:
using GPUArrays, StaticArrays
using Base.Test, GPUArrays.TestSuite

function run_indexing(Typ)
    @testset "indexing" begin
        for T in (Float32, Int32#=, SVector{3, Float32}=#)
            @testset "Indexing with $T" begin
                x = rand(T, 32)
                src = Typ(x)
                GPUArrays.allowslow(true)
                for (i, xi) in enumerate(x)
                    @test src[i] == xi
                end
                GPUArrays.allowslow(false)
                @test Array(src[1:3]) == x[1:3]
                @test Array(src[3:end]) == x[3:end]
            end
        end

        for T in (Float32, Int32)
            @testset "Indexing with $T" begin
                x = zeros(T, 7)
                src = Typ(x)
                GPUArrays.allowslow(true)
                for i = 1:7
                    src[i] = i
                end
                @test Array(src) == T[1:7;]
                src[1:3] = T[77, 22, 11]
                @test Array(src[1:3]) == T[77, 22, 11]
                src[1] = T(0)
                src[2:end] = 77
                GPUArrays.allowslow(false)
                @test Array(src) == T[0, 77, 77, 77, 77, 77, 77]
            end
        end

        for T in (Float32, Int32)
            @testset "issue #42 with $T" begin
                Ac = rand(Float32, 2, 2)
                A = Typ(Ac)
                GPUArrays.allowslow(true)
                @test A[1] == Ac[1]
                @test A[end] == Ac[end]
                @test A[1, 1] == Ac[1, 1]
                GPUArrays.allowslow(false)
            end
        end
    end
end


text:
module GLExtendedFunctions

export glErrorMessage
export glPrintError
export glCheckError
export glGenOne
export glGetActiveUniform
export glGetActiveAttrib
export glGetActiveUniformsiv
export getAttributesInfo
export getUniformsInfo
export getUniformInBlockInfo
export getProgramInfo
export glSwitch
#export glDepthMask
export glShaderSource
export glGetShaderSource
export glCompile
export glGetShaderiv
export glGetProgramiv
export glGetAttachedShaders
export getValByRef
export getValsByRef
export glString
export glGetIntegerv
export glGetInteger64v
export glGetBooleanv
export glGetFloatv
export glGetDoublev
export glGetIntegeri_v
export glGetInteger64i_v
export glGetBooleani_v
export glGetFloati_v
export glGetDoublei_v
export glGet
#export glGetInternalformati64v
#export glGetInternalformativ
export setFragLocation
export glFragLocation
export glAttribLocation
export glUniformLocation
export glUniform

using MatrixMath
using ModernGL

import ModernGL.glShaderSource
import ModernGL.glGetShaderSource
import ModernGL.glGetAttachedShaders
#import ModernGL.glGetActiveUniform
#import ModernGL.glGetActiveAttrib
import ModernGL.glGetIntegerv
#import ModernGL.glGenBuffers
#import ModernGL.glGenVertexArrays
#import ModernGL.glGenTextures
#import ModernGL.glGenFramebuffers
#import ModernGL.glGetTexLevelParameteriv
#import ModernGL.glGenRenderbuffers
#import ModernGL.glDeleteTextures
#import ModernGL.glDeleteVertexArrays
#import ModernGL.glDeleteBuffers
import ModernGL.glGetShaderiv
import ModernGL.glGetProgramiv
#import ModernGL.glViewport
#import ModernGL.glScissor

#=
glErrorMessage() = (err = glGetError();	err in keys(Errors) ? Errors[err] : "Unknown OpenGL error with error code $err.")
glPrintError(message,title="") = (length(message) > 0 ? (length(title) > 0 ? error("Error ", title, ": ", message) : error("Error: ", message)) : nothing)

function glCheckError(title="")
	err = false
	for i=1:100
		msg = glErrorMessage()
		glPrintError(msg,title)
		if msg == "" break end
		if !err err = true end
	end
	err
end

glGenOne(glGenFn,args...) = (r1 = glGenFn(args...); r2 = glCheckError(string(glGenFn)); [r2,r1])
=#

""" TODO """
function glGetActiveUniform(program::GLuint, index::Integer)
    actualLength   = GLsizei[1]
    uniformSize    = GLint[1]
    typ            = GLenum[1]
    maxcharsize    = glGetProgramiv(program, GL_ACTIVE_UNIFORM_MAX_LENGTH)
    name           = Vector{GLchar}(maxcharsize)

    glGetActiveUniform(program, index, maxcharsize, actualLength, uniformSize, typ, name)

    actualLength[1] <= 0 &&  error("No active uniform at given index. Index: ", index)

    uname = unsafe_string(pointer(name), actualLength[1])
    uname = Symbol(replace(uname, r"\[\d*\]", "")) # replace array brackets. This is not really a good solution.
    (uname, typ[1], uniformSize[1])
end

""" TODO """
function glGetActiveAttrib(program::GLuint, index::Integer)
    actualLength   = GLsizei[1]
    attributeSize  = GLint[1]
    typ            = GLenum[1]
    maxcharsize    = glGetProgramiv(program, GL_ACTIVE_ATTRIBUTE_MAX_LENGTH)
    name           = Vector{GLchar}(maxcharsize)

    glGetActiveAttrib(program, index, maxcharsize, actualLength, attributeSize, typ, name)

    actualLength[1] <= 0 && error("No active uniform at given index. Index: ", index)

    uname = unsafe_string(pointer(name), actualLength[1])
    uname = Symbol(replace(uname, r"\[\d*\]", "")) # replace array brackets. This is not really a good solution.
    (uname, typ[1], attributeSize[1])
end

""" TODO """
function glGetActiveUniformsiv(program::GLuint, index::Integer, variable::GLenum)
    result = Ref{GLint}(-1)
		#glGetActiveUniformBlockiv(program, index, variable, result)
    result[]
end

""" display information for a program's attributes """
function getAttributesInfo(program::GLuint)
	# how many attribs?
	@show activeAttr = glGetProgramiv(program, GL_ACTIVE_ATTRIBUTES)
	# get location and type for each attrib
	for i=0:activeAttr-1
		@show name, typ, siz = glGetActiveAttrib(program,	i)
		@show loc = glGetAttribLocation(program, name)
	end
end

""" display info for all active uniforms in a program """
function getUniformsInfo(program::GLuint)
	# Get uniforms info (not in named blocks)
	@show activeUnif = glGetProgramiv(program, GL_ACTIVE_UNIFORMS)

	for i=0:activeUnif-1
		@show index = glGetActiveUniformsiv(program, i, GL_UNIFORM_BLOCK_INDEX)
		if (index == -1)
			@show name 		     = glGetActiveUniformName(program, i)
			@show uniType 	   	 = glGetActiveUniformsiv(program, i, GL_UNIFORM_TYPE)

			@show uniSize 	   	 = glGetActiveUniformsiv(program, i, GL_UNIFORM_SIZE)
			@show uniArrayStride = glGetActiveUniformsiv(program, i, GL_UNIFORM_ARRAY_STRIDE)

			auxSize = 0
			if (uniArrayStride > 0)
				@show auxSize = uniArrayStride * uniSize
			else
				@show auxSize = spGLSLTypeSize[uniType]
			end
		end
	end
	# Get named blocks info
	@show count = glGetProgramiv(program, GL_ACTIVE_UNIFORM_BLOCKS)

	for i=0:count-1
		# Get blocks name
		@show name 	 		 = glGetActiveUniformBlockName(program, i)
		@show dataSize 		 = glGetActiveUniformBlockiv(program, i, GL_UNIFORM_BLOCK_DATA_SIZE)

		@show index 	 		 = glGetActiveUniformBlockiv(program, i,  GL_UNIFORM_BLOCK_BINDING)
		@show binding_point 	 = glGetIntegeri_v(GL_UNIFORM_BUFFER_BINDING, index)

		@show activeUnif   	 = glGetActiveUniformBlockiv(program, i, GL_UNIFORM_BLOCK_ACTIVE_UNIFORMS)

		indices = zeros(GLuint, activeUnif)
		glGetActiveUniformBlockiv(program, i, GL_UNIFORM_BLOCK_ACTIVE_UNIFORM_INDICES, indices)
		@show indices
		for ubindex in indices
			@show name 		   = glGetActiveUniformName(program, ubindex)
			@show uniType 	   = glGetActiveUniformsiv(program, ubindex, GL_UNIFORM_TYPE)
			@show uniOffset    = glGetActiveUniformsiv(program, ubindex, GL_UNIFORM_OFFSET)
			@show uniSize 	   = glGetActiveUniformsiv(program, ubindex, GL_UNIFORM_SIZE)
			@show uniMatStride = glGetActiveUniformsiv(program, ubindex, GL_UNIFORM_MATRIX_STRIDE)
		end
	end
end

# display the values for a uniform in a named block
""" TODO """
function getUniformInBlockInfo(program::GLuint, blockName, uniName)
	@show index = glGetUniformBlockIndex(program, blockName)
	if (index == GL_INVALID_INDEX)
		println("$uniName is not a valid uniform name in block $blockName")
	end
	@show bindIndex 		= glGetActiveUniformBlockiv(program, index, GL_UNIFORM_BLOCK_BINDING)
	@show bufferIndex 	= glGetIntegeri_v(GL_UNIFORM_BUFFER_BINDING, bindIndex)
	@show uniIndex 			= glGetUniformIndices(program, uniName)

	@show uniType 			= glGetActiveUniformsiv(program, uniIndex, GL_UNIFORM_TYPE)
	@show uniOffset 		= glGetActiveUniformsiv(program, uniIndex, GL_UNIFORM_OFFSET)
	@show uniSize 			= glGetActiveUniformsiv(program, uniIndex, GL_UNIFORM_SIZE)
	@show uniArrayStride 	= glGetActiveUniformsiv(program, uniIndex, GL_UNIFORM_ARRAY_STRIDE)
	@show uniMatStride 		= glGetActiveUniformsiv(program, uniIndex, GL_UNIFORM_MATRIX_STRIDE)
end

# display program's information
""" TODO """
function getProgramInfo(program::GLuint)
	# check if name is really a program
	@show program
	# Get the shader's name
	@show shaders = glGetAttachedShaders(program)
	for shader in shaders
		@show info = GLENUM(convert(GLenum, glGetShaderiv(shader, GL_SHADER_TYPE))).name
	end
	# Get program info
	@show info = glGetProgramiv(program, GL_PROGRAM_SEPARABLE)
	@show info = glGetProgramiv(program, GL_PROGRAM_BINARY_RETRIEVABLE_HINT)
	@show info = glGetProgramiv(program, GL_LINK_STATUS)
	@show info = glGetProgramiv(program, GL_VALIDATE_STATUS)
	@show info = glGetProgramiv(program, GL_DELETE_STATUS)
	@show info = glGetProgramiv(program, GL_ACTIVE_ATTRIBUTES)
	@show info = glGetProgramiv(program, GL_ACTIVE_UNIFORMS)
	@show info = glGetProgramiv(program, GL_ACTIVE_UNIFORM_BLOCKS)
	@show info = glGetProgramiv(program, GL_ACTIVE_ATOMIC_COUNTER_BUFFERS)
	@show info = glGetProgramiv(program, GL_TRANSFORM_FEEDBACK_BUFFER_MODE)
	@show info = glGetProgramiv(program, GL_TRANSFORM_FEEDBACK_VARYINGS)
end

""" TODO """
glSwitch(id::Symbol, on::Bool) = (option = Options[id]; on ? glEnable(option) : glDisable(option))
#glDepthMask(on)

""" TODO """
glShaderSource(shaderID::GLuint, source::String) = (shadercode=Vector{UInt8}(string(source,"\x00")); glShaderSource(shaderID, 1, Ptr{UInt8}[pointer(shadercode)], Ref{GLint}(length(shadercode))))

""" TODO """
function glGetShaderSource(shaderID::GLuint)
	len = Ref(GLint(0))
	glGetShaderiv(shaderID, GL_SHADER_SOURCE_LENGTH, len)
	len=len.x

	if len <= 0 return "" end
	println(len)

	sourceLen = len
	source = zeros(GLchar,len)
	len = zeros(GLsizei,1)
	glGetShaderSource(shaderID, sourceLen, pointer(len), pointer(source)) #(shader::GLuint, bufSize::GLsizei, length::Ptr{GLsizei}, source::Ptr{GLchar})::Void

	source=convert(String, source[1:sourceLen-1])
	println(source)
	chomp(readline())
end

""" TODO """
function glCompile(shaderID::GLuint, source::String)
	glShaderSource(shaderID, source)
	glCompileShader(shaderID)
end

""" TODO """
glGetShaderiv(shaderID::GLuint, variable::GLenum) = (result = Ref{GLint}(-1); glGetShaderiv(shaderID, variable, result); result[])

""" TODO """
glGetProgramiv(programID::GLuint, variable::GLenum) = (result = Ref{GLint}(-1); glGetProgramiv(programID, variable, result); result[])

""" TODO """
function glGetAttachedShaders(program::GLuint)
    shader_count   = glGetProgramiv(program, GL_ATTACHED_SHADERS)
    length_written = GLsizei[0]
    shaders        = zeros(GLuint, shader_count)

    glGetAttachedShaders(program, shader_count, length_written, shaders)
    shaders[1:first(length_written)]
end

""" TODO """
function getValByRef(f::Function, arg)
	ref = Ref(arg)
	f(ref)
	ref.x
end

""" TODO """
function getValsByRef(f::Function, args...)
	#refs=Array{Base.RefValue{Any},1}() # problem: different types
	#for arg in args push!(refs, Ref(arg)) end
	#f(refs...)
	#if length(refs) == 1 return refs[1] end
	#refs
	[]
end

""" TODO """
glString(name::GLenum) = unsafe_string(glGetString(name))

""" TODO """
glString(name::GLenum, index::GLuint) = unsafe_string(glGetStringi(name,index))

""" TODO """
glGetIntegerv(name::GLenum) = getValByRef((r)->glGetIntegerv(name,r), GLint(-1))

""" TODO """
glGetInteger64v(name::GLenum) = getValByRef((r)->glGetInteger64v(name,r), GLint64(-1))

""" TODO """
glGetBooleanv(name::GLenum) = getValByRef((r)->glGetBooleanv(name,r), GLboolean(false))

""" TODO """
glGetFloatv(name::GLenum) = getValByRef((r)->glGetFloatv(name,r), GLfloat(-1))

""" TODO """
glGetDoublev(name::GLenum) = getValByRef((r)->glGetDoublev(name,r), GLdouble(-1))

""" TODO """
glGetIntegeri_v(name::GLenum, index::GLuint) = getValByRef((r)->glGetIntegeri_v(name,r), GLint(-1))

""" TODO """
glGetInteger64i_v(name::GLenum, index::GLuint) = getValByRef((r)->glGetInteger64i_v(name,r), GLint64(-1))

""" TODO """
glGetBooleani_v(name::GLenum, index::GLuint) = getValByRef((r)->glGetBooleani_v(name,r), GLboolean(false))

""" TODO """
glGetFloati_v(name::GLenum, index::GLuint) = getValByRef((r)->glGetFloati_v(name,r), GLfloat(-1))

""" TODO """
glGetDoublei_v(name::GLenum, index::GLuint) = getValByRef((r)->glGetDoublei_v(name,r), GLdouble(-1))

""" TODO """
function glGet(typ::DataType, name::GLenum)

	fn = nothing
	ti = typeof(index)

	if typ == GLboolean fn = glGetBooleanv
	elseif typ == GLint fn = glGetIntegerv
	elseif typ == GLint64 fn = glGetInteger64v
	elseif typ == GLfloat fn = glGetFloatv
	elseif typ == GLdouble fn = glGetDoublev
	end

	if fn != nothing return fn(name) end #getValByRef((r)->(fn)(name,r), typ(-1))	end
	nothing
end

""" TODO """
function glGet(typ::DataType, name::GLenum, index::GLuint)

	fn = nothing
	ti = typeof(index)

	if typ == GLboolean fn = glGetBooleani_v
	elseif typ == GLint fn = glGetIntegeri_v
	elseif typ == GLint64 fn = glGetInteger64i_v
	elseif typ == GLfloat fn = glGetFloati_v
	elseif typ == GLdouble fn = glGetDoublei_v
	end

	if fn != nothing return fn(name,index) end #getValByRef((r)->(fn)(name,index,r), typ(-1))	end
	nothing
end

#glGetInternalformati64v
#glGetInternalformativ
""" TODO """
setFragLocation(id::GLuint, name::String, colorNumber=GLuint(0)) = glFragLocation(id, name, colorNumber)

""" TODO """
glFragLocation(id::GLuint, name::String, colorNumber::GLuint) = glBindFragDataLocation(id, colorNumber, pointer(string(name,"\x00")))

""" TODO """
glAttribLocation(id::GLuint, name::String) = glGetAttribLocation(id, pointer(string(name,"\x00")))

""" TODO """
glUniformLocation(id::GLuint, name::String) = glGetUniformLocation(id, pointer(string(name,"\x00")))

""" TODO """
glUniform(location::GLint, value::Bool) = glUniform1ui(location, GLuint(value?1:0))
glUniform(location::GLint, value::UInt32) = glUniform1ui(location, GLuint(value))
glUniform(location::GLint, value::UInt64) = glUniform1ui(location, GLuint(value))
glUniform(location::GLint, value::Int32) = glUniform1i(location, GLint(value))
glUniform(location::GLint, value::Int64) = glUniform1i(location, GLint(value))
glUniform(location::GLint, value::Float32) = glUniform1f(location, GLfloat(value))
glUniform(location::GLint, value::Float64) = glUniform1d(location, GLdouble(value))
glUniform(location::GLint, value::Array{Float32,1}) = glUniform1fv(location, length(value), pointer(value))
glUniform(location::GLint, value::Array{Float64,1}) = glUniform1dv(location, length(value), pointer(value))
#glUniform(location::GLint, value::Array{Vec{2, Float32},1}) = glUniform2fv(location, length(value), pointer(value))
#glUniform(location::GLint, value::Array{Vec{3, Float32},1}) = glUniform3fv(location, length(value), pointer(value))
#glUniform(location::GLint, value::Array{Vec{4, Float32},1}) = glUniform4fv(location, length(value), pointer(value))
glUniform(location::GLint, value::Vec2f) = glUniform2fv(location, 1, pointer(convert(Array, value)))
glUniform(location::GLint, value::Vec3f) = glUniform3fv(location, 1, pointer(convert(Array, value)))
glUniform(location::GLint, value::Vec4f) = glUniform4fv(location, 1, pointer(convert(Array, value)))
glUniform(location::GLint, value::Mat2x2f, transpose=false) = glUniformMatrix2fv(location, 1, transpose, pointer(convert(Array, value)))
glUniform(location::GLint, value::Mat3x3f, transpose=false) = glUniformMatrix3fv(location, 1, transpose, pointer(convert(Array, value)))
glUniform(location::GLint, value::Mat4x4f, transpose=false) = glUniformMatrix4fv(location, 1, transpose, pointer(convert(Array, value)))

#glGetSamplerParameterIiv
#glGetPointerv
#glGetQueryObjectui64v
#glGetProgramiv
#glGetBufferPointerv
#glGetTexParameterIuiv
#glGetShaderiv
#glGetRenderbufferParameteriv
#glGetUniformuiv
#glGetProgramPipelineiv
#glGetVertexAttribfv
#glGetVertexAttribLdv
#glGetVertexAttribiv
#glGetSamplerParameteriv
#glGetActiveUniformBlockiv
#glGetVertexAttribdv
#glGetQueryObjectuiv
#glGetProgramResourceiv
#glGetUniformfv
#glGetUniformdv
#glGetProgramInterfaceiv
#glGetVertexAttribIuiv
#glGetFramebufferParameteriv
#glGetActiveSubroutineUniformiv
#glGetVertexAttribPointerv
#glGetBufferParameteriv
#glGetUniformiv
#glGetQueryObjecti64v
#glGetTexLevelParameteriv
#glGetBufferParameteri64v
#glGetQueryObjectiv
#glGetActiveUniformsiv
#glGetTexParameterfv
#glGetTexLevelParameterfv
#glGetVertexAttribIiv
#glGetFramebufferAttachmentParameteriv
#glGetActiveAtomicCounterBufferiv
#glGetSynciv
#glGetSamplerParameterfv
#glGetQueryiv
#glGetTexParameterIiv
#glGetUniformSubroutineuiv
#glGetQueryIndexediv
#glGetProgramStageiv
#glGetSamplerParameterIuiv
#glGetTexParameteriv
#glGetMultisamplefv

end #GLExtendedFunctions

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
using RecipesBase
using ColorTypes

@recipe function plot(grid::AbstractCartesianMesh{1}, w, var=1)
    seriestype := :steppre
    if var isa Symbol
        title --> var
    end

    x = [cell_center(grid, j)[1] + dx(grid)[1]/2 for j in all_cells(grid)]
    data = [wi[var] for wi in w]
    x, data
end

@recipe function plot(grid::AbstractCartesianMesh{2}, w, var=1)
    seriestype := :heatmap
    seriescolor --> :viridis
    aspect_ratio --> 1
    if var isa Symbol
        title --> var
    end

    x = LinRange(grid.x_min[1], grid.x_max[1], grid.nb_cells[1])
    y = LinRange(grid.x_min[2], grid.x_max[2], grid.nb_cells[2])
    field = permutedims(map(wi -> wi[var], w), (2, 1))
    x, y, field
end

function get_base_colors(N)
    if N == 3
        return (RGB(0.643, 0.216, 0.255), RGB(0.666, 0.655, 0.224), RGB(0.176, 0.275, 0.443))
    elseif N == 4
        return (RGB(0.812, 0.408, 0.435), RGB(0.392, 0.69, 0.345), RGB(0.831, 0.757, 0.416), RGB(0.373, 0.322, 0.576))
    elseif N == 5
        return (RGB(0.812, 0.408, 0.435), RGB(0.392, 0.69, 0.345), RGB(0.831, 0.757, 0.416), RGB(0.373, 0.322, 0.576), RGB(0.2, 0.2, 0.2))
    end
end

@recipe function plot(grid::AbstractCartesianMesh{2}, w, vars::NTuple{N, Int};
                      base_colors=nothing) where N
    seriestype := :image
    yflip --> false

    if isnothing(base_colors)
        base_colors = get_base_colors(N)
    end
    base_colors = XYZ.(base_colors)

    x = LinRange(grid.x_min[1], grid.x_max[1], grid.nb_cells[1])
    y = LinRange(grid.x_min[2], grid.x_max[2], grid.nb_cells[2])
    pixels = map(wi -> sum(wi[vars[k]]*base_colors[k] for k in 1:N), w)
    field = permutedims(pixels, (2, 1))
    x, y, field
end


text:
"""
x = admm(A, B; ...)

Solves non-negative least-squares problem by the Alternating Direction Method
of Multipliers (ADMM).

Optional arguments:
    ρ: penalty parameter (set heuristically by default)
    ε: tolerance for stopping (small number times sqrt[m*n] by default)

References:
    S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein (2011). Distributed
    Optimization and Statistical Learning via the Alternating Direction Method
    of Multipliers. Foundations and Trends in Machine Learning.

    http://stanford.edu/~eryu/nnlsqr.html
"""
function admm(A,
              B::AbstractMatrix;
              ρ=max(1.0, norm(A)^2/size(A,2)),
              ε=sqrt(size(A,2)*size(B,2))*1e-15,
              max_iter=1000*size(A,2),
              kwargs...)

    # Dimensions
    m,k = size(A)
    n = size(B,2)

    # Cache matrices
    AtB = A'*B

    # Cache cholesky factorization
    L = cholesky(A'*A + ρ*I)

    # Initialize variables
    Z,U = zeros(k,n),zeros(k,n)

    # Matrix storing the solutions
    X = L \ (AtB+ρ*(Z-U))

    # Solve
    err = Inf
    for i = 1:max_iter
        Z = max(0,X+U)
        U = U+X-Z
        X = L \ (AtB+ρ*(Z-U))
        #vec norm deprecated in 0.7
        println(size(norm(X-Z)))
        norm(X-Z) < ε && break
    end

    return max(0, X)
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Functions to compute the satellite ground trace.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export ground_trace

"""
    function ground_trace(orbp::OrbitPropagator{N}, eop_data::Union{Nothing, EOPData_IAU1980, EOPData_IAU2000A} = nothing; ECI = TEME(), ECEF = PEF(), span = 1.0) where N

Compute the ground trace of the object with orbit defined by `orbp`.

By default, it considers that the orbit elements on the propagator are
represented in the True Equator, Mean Equinox (TEME) reference frame and the
ground trace will be computed in the Pseudo-Earth Fixed (PEF) reference frame.
Hence, no EOP data is needed. However, this can be changed by the keywords
presented as follows.

# Keywords

* `eop_data`: EOP data that will be used to convert the ECI reference frame to
              the ECEF reference frame. If `nothing`, then it will not be used
              (see `rECItoECEF`). (**Default** = `nothing`)
* `ECI`: ECI frame in which the orbit elements in `orbp` are represented.
         (**Default** = `TEME()`)
* `ECEF`: ECEF frame that will be used to compute the ground trace.
          (**Default** = `PEF()`)
* `span`: Defines for how much time the ground trace will be computed. The unit
          is the orbit period. (**Default** = 1.0)
* `dt`: Time interval between two samples [s]. (**Default** = 10.0)

# Returns

A vector of tuples with the pairs `(latitude,longitude)` of the ground trace.

"""
function ground_trace(orbp::OrbitPropagator{N};
                      eop_data::Union{Nothing, EOPData_IAU1980, EOPData_IAU2000A} = nothing,
                      ECI = TEME(), ECEF = PEF(), span = 1.0, dt = 10.0) where N

    # Get the orbit period.
    T = period(orbp.orb, :J2)

    # Copy orbit structure so that it is not modified by `propagate`.
    orbp_c = deepcopy(orbp)

    # Compute the points represented in the inertial reference frame.
    (o, r_i, ~) = propagate!(orbp_c, 0.0:dt:T*span)

    # Get the epochs in Julian Day of each instant.
    JD = map(x->x.t, o)

    # Convert from the ECI to the ECEF frame.
    if eop_data == nothing
        r_e = map( (t,v_i)->rECItoECEF(ECI, ECEF, t)*v_i, JD, r_i )
    else
        r_e = map( (t,v_i)->rECItoECEF(ECI, ECEF, t, eop_data)*v_i, JD, r_i )
    end

    # Convert to Geodetic coordinates.
    geod = ECEFtoGeodetic.(r_e)
    return map(x->(x[1],x[2]), geod)
end


text:
module TestMatrix

# TODO: Test `diag_dat` argument to `upload!` and `replace_coefficients!`

import ..repl_output
using AMGX, Defer, Test, JSON, CUDA, SparseArrays
using AMGX: Config, Resources, AMGXMatrix, dDDI, dFFI

@scope @testset "Matrix" begin
    c = @! Config("")
    r = @! Resources(c)
    m = @! AMGXMatrix(r, dDDI)
    @test occursin("dDDI", repl_output(m))

    @testset "upload" begin
        AMGX.upload!(m,
            Cint[0, 1, 3],
            Cint[1, 0, 1],
            [1.0, 2.0, 3.0]
        )
        @test nnz(m) == 3
        @test AMGX.matrix_get_size(m) == (2, (1,1))
        @test size(m) == (2, 2)
        @test occursin("of size 2×2 with 3 stored entries", repl_output(m))

        @testset "replace coefficients" begin
            # TODO: Should test this does something
            AMGX.replace_coefficients!(m, [2.0, 3.0, 4.0])
            @test_throws ArgumentError AMGX.replace_coefficients!(m, [2.0, 3.0, 4.0, 5.0])
            @test_throws ArgumentError AMGX.replace_coefficients!(m, Float32[2.0, 3.0, 4.0])
        end

        # Wrong element type
        m2 = @! AMGXMatrix(r, dFFI)
        @test_throws ArgumentError AMGX.upload!(m2,
            Cint[0, 1, 3],
            Cint[1, 0, 1],
            [1.0, 2.0, 3.0]
        )
    end

    @testset "blocked upload" begin
        # blocked
        blocks = [[1.0 2.0; 3.0 4.0],
                  [5.0 6.0; 7.0 8.0],
                  [9.0 10.0; 11.0 12.0]]
        blocks_flatten = collect(Iterators.flatten(blocks))
        AMGX.upload!(m,
            Cint[0, 1, 3],
            Cint[1, 0, 1],
            blocks_flatten;
            block_dims = (2,2)
        )
        @test AMGX.matrix_get_size(m) == (2, (2,2))
        @test size(m) == (4, 4)
        @test nnz(m) == 12
        @test occursin("of size 2⋅2×2⋅2 with 3 stored block entries", repl_output(m))

        @testset "replace coefficients" begin
            AMGX.replace_coefficients!(m, ones(Float64, length(blocks_flatten)))
            @test_throws ArgumentError AMGX.replace_coefficients!(m, ones(Float64, length(blocks_flatten)-1))
            @test_throws ArgumentError AMGX.replace_coefficients!(m, Float32.(blocks_flatten))
        end
    end

    @testset "upload rectangular" begin
        # rectangular
        AMGX.upload!(m,
            Cint[0, 1, 3],
            Cint[1, 0, 2],
            [1.0, 2.0, 3.0]
        )
        @test nnz(m) == 3
    end

    @testset "wrong element type" begin
        @test_throws ArgumentError AMGX.upload!(m,
            Cint[0, 1, 3],
            Cint[1, 0, 1],
            Float32[1.0, 2.0, 3.0]
        )
    end

    @scope @testset "upload device" begin
        AMGX.upload!(m,
            CuArray(Cint[0, 1, 3]),
            CuArray(Cint[1, 0, 1]),
            CuArray([1.0, 2.0, 3.0])
        )
        @test nnz(m) == 3
    end

    @testset "upload CUDA sparse matrix" begin
        c = CUDA.CUSPARSE.CuSparseMatrixCSR(sprand(Float64, 10, 10, 0.5))
        AMGX.upload!(m, c)
        @test nnz(m) == nnz(c)
        @test size(m) == size(c)
    end
end

end # module


text:
using ReachabilityAnalysis, Plots

@taylorize function lorenz!(dx, x, params, t)
    local σ = 10.0
    local β = 8.0 / 3.0
    local ρ = 28.0
    dx[1] = σ * (x[2] - x[1])
    dx[2] = x[1] * (ρ - x[3]) - x[2]
    dx[3] = x[1] * x[2] - β * x[3]
    return dx
end

X0 = Hyperrectangle(low=[0.9, 0.0, 0.0], high=[1.1, 0.0, 0.0])
prob = @ivp(x' = lorenz!(x), dim=3, x(0) ∈ X0);

alg = TMJets(abstol=1e-15, orderT=10, orderQ=2, maxsteps=50_000);

sol = solve(prob, T=10.0, alg=alg);

solz = overapproximate(sol, Zonotope);

plot(solz, vars=(0, 1), xlab="t", ylab="x")

plot(solz(0.0 .. 1.5), vars=(0, 1), xlab="t", ylab="x", lw=0.0)
plot!(x -> 20.0, c=:red, xlims=(0.0, 1.5), lab="")

ρ([1.0, 0.0, 0.0], solz(0.0 .. 1.5))

plot(solz, vars=(0, 2), xlab="t", ylab="y")

-ρ([0.0, -1.0, 0.0], solz)

ρ([0.0, 1.0, 0.0], solz)

plot(solz, vars=(0, 3), xlab="t", ylab="z")

plot(solz, vars=(1, 3))


text:
using JITrench

f(x) = 2x^2 + 5x + 10

x = Variable(2.0, name="x")
#=
name: x
values: 2.0
creator: User-Defined (nothing)
=#

y = f(x)
#=
name: nothing
values: 28.0
creator: JITrench.Add
=#

backward!(y)

@show x.grad
# x.grad = Variable(13.0)