text:
Base.eltype(::TurbulenceClosure{T}) where T = T

function Base.convert(::TurbulenceClosure{T2}, closure::TurbulenceClosure{T1}) where {T1, T2}
    paramdict = Dict((p, convert(T2, getproperty(closure, p))) for p in propertynames(closure))
    Closure = typeof(closure).name.wrapper
    return Closure(T2; paramdict...)
end

tracer_diffusivities(tracers, κ::Number) = with_tracers(tracers, (), (tracers, init) -> κ)

function tracer_diffusivities(tracers, κ::NamedTuple)

    all(name ∈ propertynames(κ) for name in tracers) ||
        throw(ArgumentError("Tracer diffusivities or diffusivity parameters must either be a constants
                            or a `NamedTuple` with a value for every tracer!"))

    return κ
end

convert_diffusivity(T, κ::Number) = convert(T, κ)
convert_diffusivity(T, κ::NamedTuple) = convert(NamedTuple{propertynames(κ), NTuple{length(κ), T}}, κ)

@inline geo_mean_Δᶠ(i, j, k, grid::RegularCartesianGrid{T}) where T = (grid.Δx * grid.Δy * grid.Δz)^T(1/3)

ext: jl


text:
using GeometryTypes, GLVisualize, GLAbstraction, ImageMagick
using FileIO, ColorTypes, Reactive

if !isdefined(:runtests)
    window = glscreen()
end
const record_interactive = true

type Mario{T}
    x 			::T
    y 			::T
    vx 			::T
    vy 			::T
    direction 	::Symbol
end

gravity(dt, mario) = (mario.vy = (mario.y > 0.0 ? mario.vy - (dt/4.0) : 0.0); mario)

function physics(dt, mario)
    mario.x = mario.x + dt * mario.vx
    mario.y	= max(0.0, mario.y + dt * mario.vy)
    mario
end

function walk(keys, mario)
    mario.vx = keys[1]
    mario.direction = keys[1] < 0.0 ? :left : keys[1] > 0.0 ? :right : mario.direction
    mario
end

function jump(keys, mario)
    if keys[2] > 0.0 && mario.vy == 0.0
    	mario.vy = 6.0
    end
	mario
end

function update(dt, keys, mario)
    mario = gravity(dt, mario)
    mario = jump(keys, 	mario)
    mario = walk(keys, 	mario)
    mario = physics(dt, mario)
    mario
end

mario2model(mario) = translationmatrix(Vec3f0(mario.x, mario.y, 0f0))*scalematrix(Vec3f0(5f0))

const mario_images = Dict()

function play(x::Vector)
	const_lift(getindex, x, loop(1:length(x)))
end

function read_sequence(path)
    if isdir(path)
        return map(load, sort(map(x->joinpath(path, x), readdir(path))))
    else
        return fill(load(path), 1)
    end
end

for verb in ["jump", "walk", "stand"], dir in ["left", "right"]
    pic = dir
    if verb != "walk" # not a sequemce
        pic *= ".png"
    end
    path = assetpath("mario", verb, pic)
    sequence = read_sequence(path)
	gif = map(img->map(RGBA{U8}, img), sequence)
	mario_images[verb*dir] = play(gif)
end
function mario2image(mario, images=mario_images)
	verb = mario.y > 0.0 ? "jump" : mario.vx != 0.0 ? "walk" : "stand"
	mario_images[verb*string(mario.direction)].value # is a signal of pictures itself (animation), so .value samples the current image
end
function arrows2vec(direction)
	direction == :up 	&& return Vec2f0( 0.0,  1.0)
	direction == :down 	&& return Vec2f0( 0.0, -1.0)
	direction == :right && return Vec2f0( 3.0,  0.0)
	direction == :left 	&& return Vec2f0(-3.0,  0.0)
	Vec2f0(0.0)
end

# Put everything together
arrows 			= sampleon(bounce(1:10), window.inputs[:arrow_navigation])
keys 			= const_lift(arrows2vec, arrows)
mario_signal 	= const_lift(update, 8.0, keys, Mario(0.0, 0.0, 0.0, 0.0, :right))
image_stream 	= const_lift(mario2image, mario_signal)
modelmatrix 	= const_lift(mario2model, mario_signal)

mario = visualize(image_stream, model=modelmatrix)

_view(mario, window, camera=:fixed_pixel)

if !isdefined(:runtests)
    renderloop(window)
end


text:
@testset "AbsAlgAss" begin

  @testset "Decomposition" begin
    Fp = GF(3)
    G = small_group(8, 4)
    FpG = group_algebra(Fp, G)

    dec = decompose(FpG)
    @test length(dec) == 5
    dim1 = 0
    dim4 = 0
    for (A, toFpG) in dec
      if dim(A) == 1
        dim1 += 1
      elseif dim(A) == 4
        dim4 += 1
      end
    end
    @test dim1 == 4
    @test dim4 == 1

    # Check whether the maps "work"
    A, toFpG = dec[1]
    e = toFpG(one(A))
    @test e^2 == e
    ee = toFpG\one(FpG)
    @test ee^2 == ee

    # And now the same for AlgAss
    FpG = AlgAss(FpG)[1]
    dec = decompose(FpG)
    @test length(dec) == 5
    dim1 = 0
    dim4 = 0
    for (A, toFpG) in dec
      if dim(A) == 1
        dim1 += 1
      elseif dim(A) == 4
        dim4 += 1
      end
    end
    @test dim1 == 4
    @test dim4 == 1

    # Check whether the maps "work"
    A, toFpG = dec[1]
    e = toFpG(one(A))
    @test e^2 == e
    ee = toFpG\one(FpG)
    @test ee^2 == ee

    # And now for AlgMat
    A = matrix_algebra(Fp, 2)

    dec = decompose(A)
    @test length(dec) == 1

    B, BtoA = dec[1]
    @test dim(B) == 4

    e = BtoA(one(B))
    @test e^2 == e
    ee = BtoA\one(A)
    @test ee^2 == ee

    A = matrix_algebra(Fp, [ matrix(Fp, [ 1 1; 0 1 ]) ]) # not semisimple!
    @test_throws AssertionError decompose(A)

    Qx, x = FlintQQ["x"]
    A = AlgAss((x^2 + 1)*(x^2 + 3))
    dec = Hecke.as_number_fields(A)

    @test length(dec) == 2
    @test degree(dec[1][1]) == 2
    @test degree(dec[2][1]) == 2

    K, AtoK = dec[1]
    e = AtoK(one(A))
    @test e^2 == e
    ee = AtoK\one(K)
    @test ee^2 == ee
  end

  @testset "Generators" begin
    Qx, x = FlintQQ["x"]
    A = AlgAss((x^2 + 1)*(x^2 + 3))
    g, full_basis, v = gens(A, Val{true})

    @test length(full_basis) == dim(A)

    M = zero_matrix(FlintQQ, dim(A), dim(A))
    for i = 1:dim(A)
      Hecke.elem_to_mat_row!(M, i, full_basis[i])
    end
    @test rank(M) == dim(A)

    for i = 1:dim(A)
      b = full_basis[i]
      a = one(A)
      for (j, k) in v[i]
        a *= g[j]^k
      end
      @test b == a
    end
  end

  @testset "Radical" begin
    Qx, x = FlintQQ["x"]
    # f = x^2 + 1
    # g = x^3 + 3x^2 + 5x - 5
    f2g3 = x^13 + 9x^12 + 44x^11 + 120x^10 + 205x^9 + 153x^8 + 32x^7 - 168x^6 - 5x^5 - 485x^4 + 500x^3 - 400x^2 + 375x - 125 # = f^2*g^3
    A = AlgAss(f2g3)
    fg = A(fmpq[-5, 5, -2, 6, 3, 1, 0, 0, 0, 0, 0, 0, 0]) # = f*g
    J = radical(A)
    I = ideal(A, fg)
    @test I == J

    f = x^2 + 1
    K, a = number_field(f, "a")
    Ky, y = K["y"]
    # g = y^3 - 3y^2 - 3y + 2
    # h = y^2 + 5y + 5
    g2h3 = y^12 + 9y^11 + 3y^10 - 198y^9 - 603y^8 + 423y^7 + 4829y^6 + 8430y^5 + 4335y^4 - 2675y^3 - 3075y^2 + 500 # = g^2*h^3
    A = AlgAss(g2h3)
    gh = A(map(K, [10, -5, -28, -13, 2, 1, 0, 0, 0, 0, 0, 0])) # = g*h
    J = radical(A)
    I = ideal(A, gh)
    @test I == J

    G = small_group(8, 4)
    F2 = GF(2)
    A = group_algebra(F2, G)
    I = radical(A)
    @test nrows(basis_matrix(I, copy = false)) == 7
    A = AlgAss(A)[1]
    I = radical(A)
    @test nrows(basis_matrix(I, copy = false)) == 7

    F3 = GF(3)
    A = group_algebra(F3, G)
    I = radical(A)
    @test nrows(basis_matrix(I, copy = false)) == 0

    F4 = GF(2, 2)
    A = group_algebra(F4, G)
    I = radical(A)
    @test nrows(basis_matrix(I, copy = false)) == 7
    A = AlgAss(A)[1]
    I = radical(A)
    @test nrows(basis_matrix(I, copy = false)) == 7

    A = group_algebra(FlintQQ, G)
    I = radical(A)
    @test nrows(basis_matrix(I, copy = false)) == 0

    for K in [ F2, F4, FlintQQ ]
      A = matrix_algebra(K, [ matrix(K, 2, 2, [ 1, 0, 0, 0 ]), matrix(K, 2, 2, [ 0, 1, 0, 0 ]), matrix(K, 2, 2, [ 0, 0, 0, 1]) ]) # i. e. upper triangular matrices
      I = radical(A)
      @test nrows(basis_matrix(I, copy = false)) == 1
    end

  end

  @testset "rand" begin
    Fp = GF(3)
    G = small_group(8, 4)
    FpG = group_algebra(Fp, G)
    A = AlgAss(FpG)[1]
    @assert A isa Hecke.AbsAlgAss

    E = elem_type(A)
    @test rand(A) isa E
    @test rand(rng, A) isa E
    @test rand(A, 2, 3) isa Matrix{E}

    Random.seed!(rng, rand_seed)
    a = rand(rng, A)
    Random.seed!(rng, rand_seed)
    @test a == rand(rng, A)
  end
end


text:
import ForwardDiff

using JuLIP: vecs, mat, JVec

export ADPotential

# Implementation of a ForwardDiff site potential
# ================================================

"""
`abstract FDPotential <: SitePotential`

To implement a site potential that uses ForwardDiff.jl for AD, create
a concrete type and overload `JuLIP.Potentials.ad_evaluate`.

Example
```julia
mutable struct P1 <: FDPotential
end
@pot P1
JuLIP.Potentials.ad_evaluate(pot::P1, R::Matrix{T}) where {T<:Real} =
               sum( exp(-norm(R[:,i])) for i = 1:size(R,2) )
```

Notes:

* For an `FDPotential`, the argument to `ad_evaluate` will be
`R::Matrix` rather than `R::Vector{<:JVec}`; this is an unfortunate current limitation of
`ForwardDiff`. Hopefully it can be fixed.
* TODO: allow arguments `(r, R)` then use chain-rule to combine them.
"""
struct ADPotential{TV, T, FR} <: SitePotential
   V::TV
   rcut::T
   gradfun::FR
end

@pot ADPotential

ADPotential(V, rcut) = ADPotential(V, rcut, ForwardDiff.gradient)

cutoff(V::ADPotential) = V.rcut

evaluate!(tmp, V::ADPotential, R::AbstractVector{<: JVec}) = V.V(R)

function evaluate_d!(dEs, tmp, V::ADPotential, R::AbstractVector{<: JVec})
   dV = V.gradfun( S -> V.V(vecs(S)), collect(mat(R)[:]) ) |> vecs
   copyto!(dEs, dV)
   return dEs
end

####### Prototype AD Hessian Implementation
# using JuLIP, JuLIP.Potentials
# using JuLIP.Potentials: evaluate_d
# JJ = JuLIP
# JPP = JuLIP.Preconditioners
#
# function ad_grad(V::SitePotential, S)
#    R = reshape(S, 3, length(S) ÷ 3) |> vecs
#    r = [norm(u) for u in R]
#    dV = evaluate_d(V, r, R)
#    return mat(dV)[:]
# end
#
# ad_hess(V::SitePotential, S::Vector{Float64}) = ForwardDiff.jacobian( S_ -> ad_grad(V, S_), S )
#
# ad_hess(V, R::AbstractVector{JVecF}) = ad_hess(V, mat(R)[:])
#
# atind2lininds(i::Integer) = (i-1) * 3 + [1,2,3]
#
# localinds(j::Integer) =  3 * (j-1) + [1,2,3]
#
# function ad_hessian(V, at::AbstractAtoms)
#    I = Int[]; J = Int[]; Z = Float64[]
#    for (i0, neigs, r, R, _) in sites(at, cutoff(V))
#       ii = atind2lininds(i0)
#       jj = [atind2lininds(j_) for j_ in neigs]
#       # compute positive version of hessian of V(R)
#       hV = ad_hess(V, R)
#
#       nneigs = length(neigs)
#       for j1 = 1:nneigs, j2 = 1:nneigs
#          LJ1, LJ2 = localinds(j1), localinds(j2)       # local indices
#          GJ1, GJ2 = jj[j1], jj[j2]                     # global indices
#          H = view(hV, LJ1, LJ2)
#          for a = 1:3, b = 1:3
#             append!(I, [ GJ1[a],    ii[a],  GJ1[a],  ii[a] ])
#             append!(J, [ GJ2[b],   GJ2[b],   ii[b],  ii[b] ])
#             append!(Z, [ H[a,b], - H[a,b], -H[a,b], H[a,b] ])
#          end
#       end
#    end
#    N = 3*length(at)
#    return sparse(I, J, Z, N, N)
# end
#
# ad_hessian(at::AbstractAtoms, x::Vector) =
#     ad_hessian(calculator(at), set_dofs!(at, x))
#

# # ========= STUFF FOR FF Preconditioner =============
#
# """
# `ADP`: this computes a basic FF preconditioner by AD-ing the site
# energy derivatives. It is slow and most of the time not as effective as
# the hand-coded ones.
# """
# type ADP{VT <: SitePotential} <: SitePotential
#  V::VT
# end
# cutoff(V::ADP) = cutoff(V.V)
#
# function _ad_grad(V::SitePotential, S)
#    R = reshape(S, 3, length(S) ÷ 3) |> vecs
#    r = [norm(u) for u in R]
#    dV = evaluate_d(V, r, R)
#    return mat(dV)[:]
# end
#
# _ad_hess(V::SitePotential, S) = ForwardDiff.jacobian( S_ -> _ad_grad(V, S_), S )
#
# function precon(Vad::ADP, r, R)
#    V = Vad.V
#    hV = _ad_hess(V, mat(R)[:])
#    # positive factorisation
#    hV = 0.5 * (hV + hV')
#    L = cholfact(Positive, hV)[:L]
#    H = L * L'
#    return 0.9 * H + 0.1 * maximum(diag(H)) * eye(size(L,1))
# end
#
# hinds(j) =  3 * (j-1) + [1,2,3]


text:
"""
    normalize_attribute_name(name)

For `String` names, this simply verifies that they pass the attribute
name production, but are otherwise untouched.

For `Symbol` names, this converts `snake_case` Symbol objects to their
`kebab-case` equivalent. So that keywords, such as `for` could be used,
we strip leading underscores.
"""
function normalize_attribute_name(name::Symbol)
    name = String(name)
    if '_' in name
       while length(name) > 0 && name[1] == '_'
           name = name[2:end]
       end
       name = replace(name, "_" => "-")
    end
    return normalize_attribute_name(name)
end

function normalize_attribute_name(name::AbstractString)
    # Attribute names are unquoted and do not have & escaping;
    # the &, % and \ characters don't seem to be prevented by the
    # specification, but they likely signal a programming error.
    for invalid in "/>='<&%\\\"\t\n\f\r\x20\x00"
        if invalid in name
            throw(DomainError(name, "Invalid character ('$invalid') " *
               "found within an attribute name."))
        end
    end
    if isempty(name)
        throw("Attribute name must not be empty.")
    end
    return name
end

"""
    rewrite_inside_tag(expr)

Attempt to speed up serialization of inside_tag by exploring the
expression tree at macro expansion time.
"""
function rewrite_inside_tag(expr)::Vector{Expr}
    if Meta.isexpr(expr, :tuple)
        return [rewrite_attribute(pair) for pair in expr.args]
    end
    if Meta.isexpr(expr, :call) && expr.args[1] == :Dict
        return [rewrite_attribute(pair) for pair in expr.args[2:end]]
    end
    return [rewrite_attribute(expr)]
end

function rewrite_attribute(expr)::Expr
    if expr isa QuoteNode && expr.value isa Symbol
        (name, value) = (expr, "")
    elseif Meta.isexpr(expr, :(=), 2)
        (name, value) = expr.args
    elseif Meta.isexpr(expr, :string, 1) && typeof(expr.args[1]) == String
        (name, value) = (expr.args[1], "")
    elseif Meta.isexpr(expr, :call, 3) && expr.args[1] == :(=>)
        (_, name, value) = expr.args
    else
        return :(inside_tag($(esc(expr))))
    end
    if name isa QuoteNode && name.value isa Symbol
        name = name.value
    end
    if name isa Expr || name isa QuoteNode
        return :(inside_tag($(esc(expr))))
    end
    attribute = normalize_attribute_name(name)
    return :(attribute_pair($attribute, $(esc(value))))
end


text:
using SimString
using Documenter

DocMeta.setdocmeta!(SimString, :DocTestSetup, :(using SimString); recursive=true)

makedocs(;
    modules=[SimString],
    authors="Bernard Brenyah",
    repo="https://github.com/PyDataBlog/SimString.jl/blob/{commit}{path}#{line}",
    sitename="SimString.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://PyDataBlog.github.io/SimString.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/PyDataBlog/SimString.jl",
    devbranch="main",
)


text:
module SSLTestOverloads
using SimpleSDMLayers
using Test

M = rand(Bool, (5, 10))
S = SimpleSDMPredictor(M, 0.0, 1.0, 0.0, 1.0)

for i in eachindex(M)
   @test S[i] == M[i]
end
for i in eachindex(S)
   @test S[i] == M[i]
end

@test typeof(S[0.2, 0.6]) == eltype(M)
@test S[0.2, 0.6] == S[Point(0.2, 0.6)]
@test isnothing(S[1.2, 0.3])
@test isnothing(S[1.2, 1.3])
@test isnothing(S[0.2, 1.3])

C = (left=0.2, bottom=0.5)
@test typeof(clip(S; C...)) == typeof(S)

Y = SimpleSDMResponse(zeros(Float64, (5,5)), 0.0, 1.0, 0.0, 1.0)
Y[0.1,0.1] = 0.2
@test Y[Point(0.1,0.1)] == 0.2

Z = convert(SimpleSDMPredictor, Y)
Y[0.1,0.1] = 4.0
@test Z[0.1,0.1] != Y[0.1,0.1]
@test Z[0.1,0.1] == 0.2

Z = convert(SimpleSDMPredictor, Y)
V = collect(Z)
@test typeof(V) == Vector{Float64}

# typed similar
c2 = similar(Y, Bool)
@test eltype(c2) == Bool

@test eltype(convert(Int64, c2)) == Int64
@test eltype(convert(Float32, c2)) == Float32

# replacement
s1 = SimpleSDMResponse(collect(reshape(1:9, 3, 3)))
@test length(s1) == 9
replace!(s1, 1 => 2, 3 => 2, 9 => nothing)
@test length(s1) == 8
@test s1.grid[1,1] == 2
@test s1.grid[3,1] == 2
@test s1.grid[1,3] == 7
@test isnothing(s1.grid[3,3])

s1 = SimpleSDMPredictor(collect(reshape(1:9, 3, 3)))
@test length(s1) == 9
s2 = replace(s1, 1 => 2, 3 => 2, 9 => nothing)
@test length(s2) == 8
@test s2.grid[1,1] == 2
@test s2.grid[3,1] == 2
@test s2.grid[1,3] == 7

# ==, isequal, hash
l1, l2 = SimpleSDMPredictor(WorldClim, BioClim, 1:2; left = 0.0, right = 10.0, bottom = 0.0, top = 10.0)
l3 = copy(l1)
l4 = similar(l1)
replace!(l4, nothing => NaN)
l5 = SimpleSDMPredictor(replace(l1.grid, nothing => missing), l1)

@test l1 == l1
@test l1 === l1
@test l2 != l1
@test l3 == l1
@test l3 !== l1

@test l4 != l1
@test l4 != l4
@test !isequal(l4, l1)
@test isequal(l4, l4)

@test ismissing(l5 == l1)
@test ismissing(l5 == l5)
@test !isequal(l5, l1)
@test isequal(l5, l5)

@test hash(l1) == hash(l1)
@test hash(l2) != hash(l1)
@test hash(l3) == hash(l1)
@test hash(l4) != hash(l1)
@test hash(l4) == hash(l4)
@test hash(l5) != hash(l1)
@test hash(l5) == hash(l5)

# getindex(layer1, layer2)
l1, l2 = SimpleSDMPredictor(WorldClim, BioClim, 1:2; left = 0.0, right = 10.0, bottom = 0.0, top = 10.0)
l3 = SimpleSDMPredictor(WorldClim, BioClim, 1; left = 5.0, right = 10.0, bottom = 5.0, top = 10.0)
@test stride(l1) == stride(l3)

l4 = clip(l1, l2)
@test l4 == l1

l5 = clip(l1, l3)
@test l5 == l3
@test_throws ArgumentError clip(l3, l1)

end


text:
function sympifyallexpressions(expressions::Array{String,1},leads_sym::Array{SymPy.Sym,1},vars::Array{String,1})
    expressions_sym = Array{SymPy.Sym}(length(expressions))

    lags = createsymversions(vars,"_lag") # Create a list of all lagged variable - these are _potential_ states
    writtenlags = createsymversions(vars,"(-1)") # Preallocating this saves time, by not having to repeatedly call Sym
    writtenleads = createsymversions(vars,"(+1)") # As above

    for (eqno,exp) in enumerate(expressions)
        Expr = sympify_catch(exp,eqno)
        expressions_sym[eqno] = Expr
    end

    subsargs = createpairedtuplelist([writtenlags writtenleads],[lags leads_sym])
    if !isempty(subsargs)
        expressions_sym = SymPy.subs(expressions_sym,subsargs...)
    end

    isstate = BitArray(length(lags))
    fill!(isstate,false)

    for expression in expressions_sym
        isstate = isstate .| expressioncontains(expression,lags)
    end

    # Now create the states array and indices array
    numstates = countnz(isstate)
    states_sym = Array{SymPy.Sym}(numstates)
    states_indices = Array{Int}(numstates)
    counter = 0
    for (i,state) in enumerate(lags)
        if isstate[i]
            counter += 1
            states_sym[counter] = state
            states_indices[counter] = i
        end
    end

    return expressions_sym,isstate,states_sym,states_indices
end

function rearrangeequations(equations)
    out = similar(equations)
    for (i,eq) in enumerate(equations)
        numequals = length(matchall(r"=", eq))
        if numequals > 1
            error("You have more than 1 equals sign in equation $(i).")
        elseif numequals == 1
            index = search(eq,'=')
            out[i] = string(eq[1:(index-1)]," - (", eq[(index+1):end] , ")")
        else
            out[i] = eq
        end
    end
    return out
end


text:
module TestAqua

using Aqua
using Reagents

test() = Aqua.test_all(Reagents; unbound_args = false)

end  # module


text:
"""
    VectorTuple{T}

A collection DataType for storing a `Tuple` of singular elements of type `T`
and/or `AbstractArray{<:T}`s in a convenient vector form that utilizes
linear indexing. Moreover `VectorTuple`s can be modified using standard vector
operations such as `empty!`, `push!`, and `deleteat!`. `VectorTuple`s should
be defined from an original tuple via `VectorTuple(tuple)` or by listing the
tuple elements `VectorTuple(items...)`. Note this is still an experimental
type and is primarily intended to store infinite parameter reference tuples and
point variable support value tuples. Some of the notable capabilities are
exemplified below.

**Example**
```julia-repl
julia> tuple = (3, [-2, 4], ones(2, 2))
(3, [-2, 4], [1.0 1.0; 1.0 1.0])

julia> vt = VectorTuple(tuple) # make by listing items (notice everything is a vector)
([3.0], [3.0, 4.0], [1.0, 1.0, 1.0, 1.0])

julia> vt[2] # linear indexing
-2.0

julia> vt[2, 2] # tuple indexing (note the second index is treated linearly)
4.0

julia> vt[6:end] # linear slicing
2-element Array{Float64,1}:
 1.0
 1.0

julia> vt[2:3, :] # tuple slicing
2-element Array{Array{Float64,1},1}:
 [-2.0, 4.0]
 [1.0, 1.0, 1.0, 1.0]

julia> tuple2 = Tuple(vt) # rebuild original Tuple with original indices
([3.0], [-2.0, 4.0], [1.0 1.0; 1.0 1.0])

julia> push!(vt, [42., 42]) # add new tuple element
([3.0], [-2.0, 4.0], [1.0, 1.0, 1.0, 1.0], [42.0, 42.0])

julia> deleteat!(vt, 4) # delete an element via linear indexing
([3.0], [-2.0, 4.0], [1.0, 1.0, 1.0])

julia> Tuple(vt) # The 3rd element becomes a SparseAxisArray because of deletion
([3.0], [-2.0, 4.0],   [1, 2]  =  1.0
  [2, 2]  =  1.0
  [2, 1]  =  1.0)

julia> deleteat!(vt, 3, tuple_index = true) # delete a whole tuple element
([3.0], [-2.0, 4.0])
```
"""
struct VectorTuple{T}
    values::Vector{T}
    ranges::Vector{UnitRange{Int}}
    indices::Vector{Any}
end

## Parse the proper indices for each container type
# Array
function _get_indices(
    arr::Array{T, N}
    )::CartesianIndices{N, NTuple{N, Base.OneTo{Int}}} where {T, N}
    return CartesianIndices(arr)
end

# DenseAxisArray
function _get_indices(arr::JuMPC.DenseAxisArray{T, N, Ax}
    )::Tuple{CartesianIndices{N, NTuple{N, Base.OneTo{Int}}}, Ax} where {T, N, Ax}
    return CartesianIndices(arr), axes(arr)
end

# SparseAxisArray
function _get_indices(arr::JuMPC.SparseAxisArray{T, N, K})::Vector{K} where {T, N, K}
    return sort!(collect(InfiniteOpt._keys(arr)))
end

# Fallback
_get_indices(arr)::Nothing = nothing

## Order arrays by index if necessary (to take care of SparseAxisArrays)
# SparseAxisArray
function _make_ordered(arr::JuMPC.SparseAxisArray{T},
                       indices::Vector)::Vector{T} where {T}
    return [arr[i] for i in indices]
end

# other
_make_ordered(arr, indices) = arr

# internal function for constructor
function _make_vt_components(tuple::Tuple)
    # form the range vector
    ranges = Vector{UnitRange{Int}}(undef, length(tuple))
    for i in eachindex(ranges)
        element_length = length(tuple[i])
        if i === 1
            ranges[i] = UnitRange(1, element_length)
        else
            ranges[i] = UnitRange(ranges[i-1].stop + 1, ranges[i-1].stop + element_length)
        end
    end
    # get the indices and reorder the tuple if needed
    indices = Any[_get_indices(tuple[i]) for i in eachindex(tuple)]
    if any(k isa JuMPC.SparseAxisArray for k in tuple)
        tuple = Tuple(_make_ordered(tuple[i], indices[i]) for i in eachindex(tuple))
    end
    # expand the elements into a single vector
    vals = [j for i in tuple for j in i]
    return vals, ranges, indices
end

# Constructor from a Tuple with type
function VectorTuple{T}(tuple::Tuple)::VectorTuple{T} where {T}
    if isempty(tuple)
        return VectorTuple{T}()
    end
    vals, ranges, indices = _make_vt_components(tuple)
    return VectorTuple{T}(vals, ranges, indices)
end

# Constructor from a Tuple without type
function VectorTuple(tuple::Tuple)::VectorTuple
    if isempty(tuple)
        return VectorTuple()
    end
    vals, ranges, indices = _make_vt_components(tuple)
    return VectorTuple(vals, ranges, indices)
end

# Constructor from various arguments (like splatting the tuple) with type
function VectorTuple{T}(items...)::VectorTuple{T} where {T}
    return VectorTuple{T}(items)
end

# Constructor from various arguments (like splatting the tuple) without type
function VectorTuple(items...)::VectorTuple
    return VectorTuple(items)
end

# Default Constructor with type
function VectorTuple{T}()::VectorTuple{T} where {T}
    return VectorTuple(T[], UnitRange{Int}[], Any[])
end

# Default Constructor without type
function VectorTuple()::VectorTuple
    return VectorTuple(Any[], UnitRange{Int}[], Any[])
end

# Extend Base.:(==)
function Base.:(==)(vt1::VectorTuple, vt2::VectorTuple)::Bool
    return all(getproperty(vt1, f) == getproperty(vt2, f) for f in fieldnames(VectorTuple))
end

# Extend Base.size
function Base.size(vt::VectorTuple, d::Int)::Int
    d > 1 && throw(ArgumentError("invalid VectorTuple dimension $d"))
    return length(vt.ranges)
end

# Determine if 2 VectorTuples have the same structure
function same_structure(vt1::VectorTuple, vt2::VectorTuple;
                        use_indices = true)::Bool
    same_size = vt1.ranges == vt2.ranges
    return same_size && (!use_indices || (vt1.indices == vt2.indices))
end

# Extend Base.copy
function Base.copy(vt::VectorTuple{T})::VectorTuple{T} where {T}
    return VectorTuple(copy(vt.values), copy(vt.ranges), copy(vt.indices))
end

## Extend first and last indices
# Base.firstindex (linear)
Base.firstindex(vt::VectorTuple)::Int = firstindex(vt.values)

# Base.firstindex (tuple)
Base.firstindex(vt::VectorTuple, d::Int)::Int = 1

# Base.lastindex (linear)
Base.lastindex(vt::VectorTuple)::Int = lastindex(vt.values)

# Base.lastindex (tuple)
function Base.lastindex(vt::VectorTuple, d::Int)::Int
    return d == 1 ? size(vt, 1) : throw(ArgumentError("lastindex(::VectorTuple, " *
                                                      "d::Int > 1) not defined."))
end

## Extend Base.getindex
# linear index along vector
Base.getindex(vt::VectorTuple, i) = getindex(vt.values, i)

# Tuple indexing (single position)
function Base.getindex(vt::VectorTuple, i::Int, j)
    return getindex(vt.values[vt.ranges[i]], j)
end

# Tuple indexing (with a colon)
function Base.getindex(vt::VectorTuple, i::Colon, j)
    return [getindex(vt, i, j) for i = 1:size(vt, 1)]
end

# Tuple indexing (with an iterable) --> assumed backup
function Base.getindex(vt::VectorTuple, is, j)
    return [getindex(vt, i, j) for i in is]
end

## Extend Base.setindex!
# linear index along vector
function Base.setindex!(vt::VectorTuple, v, i)
    return setindex!(vt.values, v, i)
end

# Tuple indexing (with integer i and colon j)
function Base.setindex!(vt::VectorTuple, v, i::Int, j::Colon)
    return setindex!(vt.values, v, vt.ranges[i])
end

# Tuple indexing (with integer i and any non colon j)
function Base.setindex!(vt::VectorTuple, v, i::Int, j)
    linear_index = (vt.ranges[i].start - 1) .+ j
    if any(linear_index .> vt.ranges[i].stop)
        throw(BoundsError(vt, [i, j]))
    end
    return setindex!(vt.values, v, linear_index)
end

# Tuple indexing (with colon i and integer j)
function Base.setindex!(vt::VectorTuple, v, i::Colon, j::Int)
    linear_indexes = [vt.ranges[k].start - 1 + j for k in eachindex(vt.ranges)]
    return setindex!(vt.values, v, linear_indexes)
end

# Tuple indexing (with colon i and noninteger j)
function Base.setindex!(vt::VectorTuple, v, i::Colon, j)
    return setindex!(vt, v, 1:size(vt, 1), j)
end

# Tuple indexing (with array i and any j)
function Base.setindex!(vt::VectorTuple, v, i, j)
    @assert length(v) == length(i)
    return [setindex!(vt, v[k], i[k], j) for k in 1:length(i)]
end

# Extend simple Base redirect methods with VectorTuple arguments
for op = (:length, :isempty, :eachindex, :keys)
    @eval Base.$op(vt::VectorTuple) = $op(vt.values)
end

# Extend simple Base 2 argument methods
for op = (:findfirst, :findall)
    @eval Base.$op(f::Function, vt::VectorTuple) = $op(f, vt.values)
end

# Extend Base.in
Base.in(i, vt::VectorTuple)::Bool = in(i, vt.values)

# Extend Base.iterate
Base.iterate(vt::VectorTuple) = iterate(vt.values)
Base.iterate(vt::VectorTuple, i) = iterate(vt.values, i)

# Extend Base.empty!
function Base.empty!(vt::VectorTuple{T})::VectorTuple{T} where {T}
    empty!(vt.values)
    empty!(vt.ranges)
    empty!(vt.indices)
    return vt
end

## Extend Base.push!
# Single item
function Base.push!(vt::VectorTuple{T}, item::T)::VectorTuple{T} where {T}
    push!(vt.values, item)
    push!(vt.ranges, UnitRange(vt.ranges[end].stop + 1, vt.ranges[end].stop + 1))
    push!(vt.indices, nothing)
    return vt
end

# AbstractArray
function Base.push!(vt::VectorTuple{T},
                    items::AbstractArray{<:T})::VectorTuple{T} where {T}
    indices = _get_indices(items)
    push!(vt.values, _make_ordered(items, indices)...)
    push!(vt.ranges, UnitRange(vt.ranges[end].stop + 1, length(vt.values)))
    push!(vt.indices, indices)
    return vt
end

## Update the indices if an element is deleted
# SparseAxisArray
function _update_indices(inds::Vector{<:Tuple}, i)
    new_inds = deleteat!(inds, i)
    return length(new_inds) >= 2 ? new_inds : nothing
end

# DenseAxisArray
function _update_indices(inds::Tuple, i)
    new_inds = reduce(vcat, Base.Iterators.product(inds[2]...))
    return _update_indices(new_inds, i) # call the SparseAxisArray version
end

# Array
function _update_indices(inds::CartesianIndices, i)
    new_inds = [Tuple(inds[k]) for k in 1:length(inds)]
    return _update_indices(new_inds, i) # call the SparseAxisArray version
end

# Extend Base.deleteat!
function Base.deleteat!(vt::VectorTuple{T}, i::Int;
                        tuple_index::Bool = false)::VectorTuple{T} where {T}
    if tuple_index
        len = length(vt.ranges[i])
        start = vt.ranges[i].start
        deleteat!(vt.values, vt.ranges[i])
        deleteat!(vt.ranges, i)
        deleteat!(vt.indices, i)
        for j in eachindex(vt.ranges)
            if vt.ranges[j].start > start
                vt.ranges[j] = vt.ranges[j] .- len
            end
        end
    else
        group = findfirst(x->in(i, x), vt.ranges)
        deleteat!(vt.values, i)
        if length(vt.ranges[group]) == 1
            deleteat!(vt.ranges, group)
            deleteat!(vt.indices, group)
        else
            vt.indices[group] = _update_indices(vt.indices[group],
                                                i - vt.ranges[group].start + 1)
        end
        for j in eachindex(vt.ranges)
            if vt.ranges[j].start > i
                vt.ranges[j] = vt.ranges[j] .- 1
            elseif vt.ranges[j].stop >= i
                vt.ranges[j] = UnitRange(vt.ranges[j].start, vt.ranges[j].stop - 1)
            end
        end
    end
    return vt
end

# Extend Base.deleteat! (boolean list)
function Base.deleteat!(vt::VectorTuple{T}, inds::AbstractVector{<:Bool};
                        tuple_index::Bool = false)::VectorTuple{T} where {T}
    if tuple_index
        deleteat!(vt.values, [j for i in vt.ranges[inds] for j in i])
        deleteat!(vt.indices, inds)
        prev_sum = 0
        for i in eachindex(vt.ranges)
            if inds[i]
                prev_sum += length(vt.ranges[i])
            else
                vt.ranges[i] = vt.ranges[i] .- prev_sum
            end
        end
        deleteat!(vt.ranges, inds)
    else
        deleteat!(vt.values, inds)
        prev_sum = 0
        delete_inds = []
        for i in eachindex(vt.ranges)
            delete_sum = sum(inds[vt.ranges[i]])
            if length(vt.ranges[i]) == delete_sum
                push!(delete_inds, i)
            else
                group_inds = vt.ranges[i][inds[vt.ranges[i]]] .+ (1 - vt.ranges[i].start)
                if length(group_inds) != 0
                    vt.indices[i] = _update_indices(vt.indices[i], group_inds)
                end
                vt.ranges[i] = UnitRange(vt.ranges[i].start - prev_sum,
                                         vt.ranges[i].stop - prev_sum - delete_sum)
            end
            prev_sum += delete_sum
        end
        deleteat!(vt.ranges, delete_inds)
        deleteat!(vt.indices, delete_inds)
    end
    return vt
end

# TODO add deleteat! for sequential indices

## Make the original array for a particular tuple element
# Array
function _make_array(values::Vector{T},
                     indices::CartesianIndices{N})::Array{T, N} where {T, N}
    return [values[i] for i in LinearIndices(indices)]
end

# DenseAxisArray
function _make_array(values::Vector, indices::Tuple)
    array = _make_array(values, indices[1])
    return JuMPC.DenseAxisArray(array, indices[2]...)
end

# SparseAxisArray
function _make_array(values::Vector, indices::Vector)
    data = Dict(indices[i] => values[i] for i in eachindex(indices))
    return JuMPC.SparseAxisArray(data)
end

# Other
function _make_array(values::Vector, indices::Nothing)
    return length(values) == 1 ? first(values) : values
end

# Define Tuple construction from a VectorTuple
function Base.Tuple(vt::VectorTuple; use_indices = true)::Tuple
    if use_indices
        return Tuple(_make_array(vt.values[vt.ranges[i]], vt.indices[i])
                     for i in eachindex(vt.ranges))
    else
        return Tuple(length(vt.ranges[i]) == 1 ? vt.values[vt.ranges[i].start] : vt.values[vt.ranges[i]]
                     for i in eachindex(vt.ranges))
    end
end

# Define Tuple using construction from a values vector in combination with a VT
function Base.Tuple(values::Vector, vt::VectorTuple; use_indices = true)::Tuple
    @assert length(values) == length(vt)
    if use_indices
        return Tuple(_make_array(values[vt.ranges[i]], vt.indices[i])
                     for i in eachindex(vt.ranges))
    else
        return Tuple(length(vt.ranges[i]) == 1 ? values[vt.ranges[i].start] : values[vt.ranges[i]]
                     for i in eachindex(vt.ranges))
    end
end

## Enable pretty printing
# Extend Base.string
function Base.string(vt::VectorTuple)::String
    return string(Tuple(vt, use_indices = false))
end

# Show VectorTuples in REPLMode
function Base.show(io::IO, vt::VectorTuple)
    print(io, string(vt))
    return
end

# Show VectorTuples in IJuliaMode
function Base.show(io::IO, ::MIME"text/latex", vt::VectorTuple)
    print(io, string(vt))
end


text:
module SymbolServer

export SymbolServerProcess
export getstore

using Serialization, Pkg, SHA

include("clientprocess/from_static_lint.jl")

mutable struct SymbolServerProcess
    process::Base.Process
    context::Union{Nothing,Pkg.Types.Context}
    depot::Dict{String,ModuleStore}
    process_stderr::Union{IOBuffer,Nothing}

    function SymbolServerProcess(;environment=nothing, depot=nothing)
        jl_cmd = joinpath(Sys.BINDIR, Base.julia_exename())
        client_process_script = joinpath(@__DIR__, "clientprocess", "clientprocess_main.jl")

        env_to_use = copy(ENV)
        if depot!==nothing
            if depot==""
                delete!(env_to_use, "JULIA_DEPOT_PATH")
            else
                env_to_use["JULIA_DEPOT_PATH"] = depot
            end
        end

        stderr_for_client_process = VERSION < v"1.1.0" ? nothing : IOBuffer()

        p = if environment===nothing
            open(pipeline(Cmd(`$jl_cmd $client_process_script`, env=env_to_use), stderr=stderr_for_client_process), read=true, write=true)
        else
            open(pipeline(Cmd(`$jl_cmd --project=$environment $client_process_script`, dir=environment, env=env_to_use), stderr=stderr_for_client_process), read=true, write=true)
        end
        ssp = new(p, nothing, Dict(), stderr_for_client_process)
        get_context(ssp)
        return ssp
    end
end

function request(server::SymbolServerProcess, message::Symbol, payload)
    serialize(server.process, (message, payload))
    ret_val = try
        deserialize(server.process)
    catch err
        # Only Julia 1.1 and newer support capturing stderr into an IOBuffer
        if server.process_stderr!==nothing
            stderr_from_client_process = String(take!(server.process_stderr))

            complete_error_message = string(sprint(showerror, err), "\n\nstderr from client process:\n\n", stderr_from_client_process)

            error(complete_error_message)
        else
            complete_error_message = string(sprint(showerror, err), "\n\nCouldn't capture stderr from client process on julia 1.0.\n\n")

            error(complete_error_message)
        end
    end

    !(ret_val isa Tuple{Symbol,<:Any}) && error("Invalid response:\n", ret_val)
    return ret_val
end

function get_context(server::SymbolServerProcess)
    status, payload = request(server, :get_context, nothing)
    if status == :success
        server.context = payload
        return
    else
        error(payload)
    end
end

function load_store_from_disc(file::String)
    io = open(file)
    store = deserialize(io)
    close(io)
    typeof(store) != ModuleStore && @info "WARNING: Type mismatch in loaded store"
    return store
end

function shouldreload(pkgid, pkg, c)
    path = pkg_path(pkgid, c)
    isempty(path) && return false
    !endswith(pkg.ver, "+") && return false
    !isgitrepo(path) && return false
    return getgithash(path) != pkg.sha
end

function safe_load_store(pkg::PackageID, server::SymbolServerProcess, allowfail = true)
    storedir = abspath(joinpath(@__DIR__, "..", "store"))
    try
        server.depot[pkg.name] = load_store_from_disc(joinpath(storedir, "$(pkg.uuid).jstore"))
        !(server.depot[pkg.name] isa ModuleStore) && error("Type mismatch")

        if shouldreload(pkg, server.depot[pkg.name], server.context)
            parents = Base.UUID(pkg.uuid) in keys(server.context.env.manifest) ? [pkg] : find_parent(server.context, pkg.uuid)
            isempty(parents) && return
            loaded_pkgs = load_package(server, first(parents))
            for pkg1 in loaded_pkgs
                if haskey(server.context.env.manifest, Base.UUID(pkg1[1]))
                    safe_load_store(PackageID(pkg1[2], pkg1[1]), server, false)
                end
            end
        end
    catch e
        !allowfail && return
        parents = Base.UUID(pkg.uuid) in keys(server.context.env.manifest) ? [pkg] : find_parent(server.context, pkg.uuid)
        isempty(parents) && return
        loaded_pkgs = load_package(server, first(parents))
        for pkg1 in loaded_pkgs
            if haskey(server.context.env.manifest, Base.UUID(pkg1[1]))
                safe_load_store(PackageID(pkg1[2], pkg1[1]), server, false)
            end
        end
    end
end# loaded_pkgs = load_package(server, pkg)

# Public API

function getstore(server::SymbolServerProcess)
    storedir = abspath(joinpath(@__DIR__, "..", "store"))
    try
        server.depot["Base"] = load_store_from_disc(joinpath(storedir, "Base.jstore"))
        server.depot["Core"] = load_store_from_disc(joinpath(storedir, "Core.jstore"))
    catch e
        load_core(server)
        try
            server.depot["Base"] = load_store_from_disc(joinpath(storedir, "Base.jstore"))
            server.depot["Core"] = load_store_from_disc(joinpath(storedir, "Core.jstore"))
        catch e
            error("Couldn't load core stores")
        end
    end
    for pkg in get_manifest(server.context)
        pkg.name in keys(server.depot) && continue
        safe_load_store(pkg, server)
    end
    return server.depot
end

function Base.kill(server::SymbolServerProcess)
    serialize(server.process, (:close, nothing))
end

function load_core(server::SymbolServerProcess)
    status, payload = request(server, :load_core, nothing)
    if status == :success
        return payload
    else
        error(payload)
    end
end

function load_package(server::SymbolServerProcess, pkg::PackageID)
    status, payload = request(server, :load_package, (pkg.name, pkg.uuid))
    if status == :success
        return payload
    else
        error(payload)
    end
end

function load_all(server::SymbolServerProcess)
    status, payload = request(server, :load_all, nothing)
    if status == :success
        return payload
    else
        error(payload)
    end
end

function clear_disc_store()
    storedir = abspath(joinpath(@__DIR__, "..", "store"))
    for f in readdir(storedir)
        if endswith(f, ".jstore")
            rm(joinpath(storedir, f))
        end
    end
end

end # module


text:
################################################################################
###                    Parallel siman algorithm                           ######
################################################################################

#Beta: for running ATA with multiple cores do:
#add julia.exe to the windows path and run Julia with the command prompt cmd:
#julia -p <numberOfCores>

using Distributed  #this is not needed if Julia has been run with <numberOfCores> >1
@everywhere using ATA
@everywhere cd("where your input files are")

#Input files needed: settings_ata.jl, bank.csv, Categoricalconstraints.csv, overlap_matrix.csv, BSpar.jld2 (only for Chance-Contrained (CC))
#settings_ata.jl : overall features of the tests, such as length, expected score, item use, ecc...
#Categoricalconstraints.csv : categorical constraints
#overlap_matrix.csv : maximum overlap allowed between test forms, it is a TxT matrix. If the assembly is non parallel (i.e. there is more than one group of parallel tests, n_groups>1)
#it is a n_groups x n_groups matrix.
#BSpar.jld2 : it is a nPar way array (Array{Float64,nPar}) where nPar is the number of IRT parameters. Each sub array is a n_items x R matrix (Matrix{Float64}(.,n_items,R)).

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

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

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

# 4. Add objective function (Optional)
add_obj_fun!(ata_model);
print_infos(ata_model)

#Assembly settings

# Set the solver, "siman" for simulated annealing, "jump" for MILP solver.
solver = "siman";

# SIMAN (Suggested for Large scale ATA):

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

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

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

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

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

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

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

#! Termination criteria:

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

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

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

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

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

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

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

#to stop all the processes do:
#ctrl+C
#interrupt()
#addprocs(numberOfProcs)
#@everywhere cd("where your input files are")
#include("LocalTest.jl")
#Julia will compile the package again.


text:
module TunedSmoother

include("candidate_smoothers.jl")
export Jacobi!,GS!,SOR!,pseudo!,state,precond!,tunedJ!

include("create_synthetic.jl")
export create_synthetic,static,dipole,sphere

include("create_waterlily.jl")
export create_waterlily,circle,TGV,donut,wing,shark

include("metrics.jl")
export itcount,avecount,loss,fit,mg!

end # module


text:
# ==============================================================================
# Xavier Gandibleux - November 2021
#   Implemented in Julia 1.6

# ==============================================================================
using PyPlot

# ------------------------------------------------------------------------------
# kung Algorithm for extracting SN from a static set S of points
function kungAlgorithm(S1, S2)
    S = []
    for i=1:length(S1)
        push!(S, (S1[i] , S2[i]) )
    end
    sort!(S, by = x -> x[1])
    SN=[] ; push!(SN, S[1]) ; minS2 = S[1][2]
    for i=2:length(S1)
        if S[i][2] < minS2
            push!(SN, S[i]) ; minS2 = S[i][2]
        end
    end
    return SN
end

# ------------------------------------------------------------------------------
# compute the corner points of a set (S1,S2) and its lower envelop
function computeCornerPointsLowerEnvelop(S1, S2)
    # when all objectives have to be maximized
    Env1=[]; Env2=[]
    for i in 1:length(S1)-1
        push!(Env1, S1[i]); push!(Env2, S2[i])
        push!(Env1, S1[i]); push!(Env2, S2[i+1])
    end
    push!(Env1, S1[end]);push!(Env2, S2[end])
    return Env1,Env2
end

# ------------------------------------------------------------------------------
# display different results
function displayGraphics(fname,YN)
    PlotOrthonormedAxis = true  # Axis orthonormed or not
    DisplayYN   = true          # Non-dominated points corresponding to efficient solutions
    DisplayUBS  = false         # Points belonging to the Upper Bound Set
    DisplayLBS  = false         # Points belonging to the Lower Bound Set
    DisplayInt  = false         # Points corresponding to integer solutions
    DisplayProj = false         # Points corresponding to projected solutions
    DisplayFea  = false         # Points corresponding to feasible solutions
    DisplayPer  = false         # Points corresponding to perturbated solutions

    YN_1=[];YN_2=[]
    for i in 1:length(YN)
        push!(YN_1, YN[i][1])
        push!(YN_2, YN[i][2])
    end

    # --------------------------------------------------------------------------
    # Setup
    figure("Project MOMH 2021",figsize=(6.5,5))
    if PlotOrthonormedAxis
        vmin = 0.99 * min(minimum(YN_1),minimum(YN_2))
        vmax = 1.01 * max(maximum(YN_1),maximum(YN_2))
        xlim(vmin,vmax)
        ylim(vmin,vmax)
    end
    xlabel(L"z^1(x)")
    ylabel(L"z^2(x)")
    PyPlot.title("Bi-01BKP | $fname")

    # --------------------------------------------------------------------------
    # Display Non-Dominated points
    if DisplayYN
        # display only the points corresponding to non-dominated points
        scatter(YN_1, YN_2, color="black", marker="+", label = L"y \in Y_N")
        # display segments joining adjacent non-dominated points
        plot(YN_1, YN_2, color="black", linewidth=0.75, marker="+", markersize=1.0, linestyle=":")
        # display segments joining non-dominated points and their corners points
        Env1,Env2 = computeCornerPointsLowerEnvelop(YN_1, YN_2)
        plot(Env1, Env2, color="black", linewidth=0.75, marker="+", markersize=1.0, linestyle=":")
    end

    # --------------------------------------------------------------------------
    # Display a Upper bound set (primal, by default)
    if DisplayUBS
        plot(xU, yU, color="green", linewidth=0.75, marker="+", markersize=1.0, linestyle=":")
        scatter(xU, yU, color="green", label = L"y \in U", s = 150,alpha = 0.3)
        scatter(xU, yU, color="green", marker="o", label = L"y \in U")
    end

    # --------------------------------------------------------------------------
    # Display a Lower bound set (dual, by excess)
    if DisplayLBS
        plot(xL, yL, color="blue", linewidth=0.75, marker="+", markersize=1.0, linestyle=":")
        scatter(xL,yL, color="blue", marker="x", label = L"y \in L")
    end

    # --------------------------------------------------------------------------
    # Display integer points (feasible and non-feasible in GravityMachine)
    if DisplayInt
        scatter(XInt,YInt, color="orange", marker="s", label = L"y"*" rounded")
    end

    # --------------------------------------------------------------------------
    # Display projected points (points Δ(x,x̃) in GravityMachine)
    if DisplayProj
        scatter(XProj,YProj, color="red", marker="x", label = L"y"*" projected")
    end

    # --------------------------------------------------------------------------
    # Display feasible points
    if DisplayFea
        scatter(XFeas,YFeas, color="green", marker="o", label = L"y \in F")
    end

    # --------------------------------------------------------------------------
    # Display perturbed points (after a cycle in GravityMachine)
    if DisplayPer
        scatter(XPert,YPert, color="magenta", marker="s", label ="pertub")
    end

    # --------------------------------------------------------------------------
    legend(bbox_to_anchor=[1,1], loc=0, borderaxespad=0, fontsize = "x-small")
end

# ==============================================================================
# Example on how to use it :

#fname="example"
#YN = [[902.0, 1193.0], [962.0, 1187.0], [967.0, 1186.0], [994.0, 1163.0],
#[1080.0, 1155.0], [1108.0, 1131.0], [1117.0, 1122.0], [1119.0, 1118.0],
#[1142.0, 1083.0], [1148.0, 1056.0], [1159.0, 1044.0], [1167.0, 1005.0],
#[1180.0, 1004.0], [1181.0, 1003.0], [1182.0, 979.0], [1198.0, 977.0],
#[1219.0, 969.0], [1220.0, 818.0]]

#displayGraphics(fname,YN)

max_stars_repo_name: xgandibleux/momhMKPstu

max_issues_repo_name: xgandibleux/momhMKPstu

max_forks_repo_name: xgandibleux/momhMKPstu


text:
# Kernel module
module Kernel

export GaussianKernel!, GaussianKernel, GaussianKernelPreImage, GaussianKernelDenoise, LinearKernel, PolynomialKernel, ANOVAKernel, ProjectionError, IPKernel, InverseIPKernel, AEPFeature, AEPKernel, AEPCovariance, InverseAEPKernel, LinearMean, LogMap, ExpMap, FrechetMean, PGAFeature, PGAKernel, PGACovariance, InversePGAKernel, SphericalFeature, SphericalFeatureInverse, SphericalKernel, SphericalKernelLowMem, SphericalCosKernel, SphericalCovariance, SphericalCovarianceLowMem

# Kernel definitions
include("Kernels.jl")

# Kerne preimage definitions
include("PreImage.jl")

# Kernels specific to orientation data
include("KernelOrientation.jl")

# Computes projection error. Uses a TSVD approximation if necessary.
# @param Kxx - K(x,x) the feature space inner product of x with itself.
# @param Kx -  K(:,x) the feature space inner produce of x with all atoms of K.
# @param K - K(:,:) the kernel matrix / gram matrix in feature space of all atoms - this may be a rank-k approximation.
# @param S,U - Eigen decomposition of K - these may be rank-k approximations
# @param condK - condition number of K
# @param epsilon - the TSVD cut off.
function ProjectionError(Kxx, Kx, K, S, U, condK, epsilon)
	proj_err = 1e100
	if condK > 1e4 # remember that ||x-x'||/||x|| <= cond(A) * ||Ax'-b||/||b||, where x is the exact solution to Ax=b
		# attempt at computing this more carefully:
		Utkj = U' * Kx
		S[ find( S .< epsilon ) ] = 0
		S[ find( S .> epsilon ) ] = 1 ./ S[ find( S .> 1e-30 ) ]
		S = sqrt(S)
		Utkj = diagm(S) * Utkj
		proj_err = (Kxx - Utkj'*Utkj)[1,1]
	else
		proj_err = (Kxx - (Kx' * (K \ Kx)))[1,1]
	end

	proj_err = sqrt( 2 * abs(proj_err) )
	return proj_err
end

end # module


text:
function roundTT(mpt::MPT, ϵ::Float64)

    N           = order(mpt);
    cores       = Vector{Array{Float64,3}}(undef,N);
    middlesizes = size(mpt);
    frobnorm    = norm(mpt);
    err2        = 0;

    # Right-to-left orthogonalization
        G = mpt[N];
        for k = N:-1:2
            G        = unfold(G,[1],"left");
            F        = qr(G');
            Q        = Matrix(F.Q);
            R        = Matrix(F.R);
            G        = reshape(Q', ( size(Q,2), middlesizes[k][1], Int(length(Q)/(size(Q,2)*middlesizes[k][1])) ) );
            cores[k] = G;
            G        = nmodeproduct(R,mpt[k-1],3);
        end
        cores[1] = G;
        # Compression of the orthogonalized representation
        δ = ϵ / sqrt(N-1) * frobnorm;
        for k = 1:N-1
            F     = svd(unfold(cores[k],[3],"right"),full=false);
            rcurr = length(F.S);
            sv2   = cumsum(reverse(F.S).^2);
            tr    = Int(findfirst(sv2 .> δ^2))-1;
            if tr > 0
                rcurr = length(F.S) - tr;
                err2 += sv2[tr];
            end
            Utr = F.U[:,1:rcurr];
            Str = F.S[1:rcurr];
            Vtr = F.Vt[1:rcurr,:];

            cores[k]   = reshape( Utr,( Int(length(Utr)/(size(cores[k],2)*size(Vtr,1))), size(cores[k],2), size(Vtr,1) ) );
            cores[k+1] = nmodeproduct(Matrix((Vtr'*Diagonal(Str))'),cores[k+1],1);
        end
        return MPT(cores);
    end

    function roundTT(ttm::MPT{4}, ϵ::Float64)
        tt = mpo2mps(ttm);
        tt = roundTT(tt, ϵ)
        return mps2mpo(tt,size(ttm))
    end

max_stars_repo_name: clarazen/BigMat

max_issues_repo_name: clarazen/BigMat

max_forks_repo_name: clarazen/BigMat


text:
if basename(pwd()) == "aoc"
    cd("2019/22")
end

function shuffle(filename::AbstractString, len)
    cards = [0:len-1...]
    for line in eachline(filename)
        if startswith(line, "deal with increment")
            n = parse(Int, line[21:end])
            cards[[((i - 1) * n) % len + 1 for i in eachindex(cards)]] .= cards
        elseif startswith(line, "deal into new stack")
            reverse!(cards)
        elseif startswith(line, "cut")
            cards = circshift(cards, -parse(Int, line[5:end]))
        end
    end
    cards
end

@assert shuffle("example1.txt", 10) == [0, 3, 6, 9, 2, 5, 8, 1, 4, 7]
@assert shuffle("example2.txt", 10) == [3, 0, 7, 4, 1, 8, 5, 2, 9, 6]
@assert shuffle("example3.txt", 10) == [6, 3, 0, 7, 4, 1, 8, 5, 2, 9]
@assert shuffle("example4.txt", 10) == [9, 2, 5, 8, 1, 4, 7, 0, 3, 6]

function part1(filename::AbstractString, len)
    findfirst(==(2019), shuffle(filename, len)) - 1
end

part1("input.txt", 10007)

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
using TableIOInterface
using Test

@testset "TableIOInterface.jl" begin

    fnames = ["test.csv", "test.jdf", "test.parquet", "test.arrow", "test.xlsx", "test.json",
                "test.zip", "test.sqlite", "test.hdf", "test.jld2"]

    @testset "File type ID" begin

        expected_file_types = [TableIOInterface.CSVFormat, TableIOInterface.JDFFormat, TableIOInterface.ParquetFormat,
        TableIOInterface.ArrowFormat, TableIOInterface.ExcelFormat, TableIOInterface.JSONFormat,
        TableIOInterface.ZippedFormat, TableIOInterface.SQLiteFormat, TableIOInterface.HDF5Format,
        TableIOInterface.JLD2Format]

        for (i, fname) ∈ enumerate(fnames)
            file_type = get_file_type(fname)
            @test file_type == expected_file_types[i]()
        end
    end

    @testset "Sample code generation" begin
        expected_sample_code = [
            """df_test_csv = let
                import TableIO
                import DataFrames
                DataFrames.DataFrame(TableIO.read_table(joinpath(split(@__FILE__, '#')[1] * ".assets", "test.csv")); copycols=false)
            end""",

            """df_test_jdf = let
                import TableIO
                import DataFrames
                DataFrames.DataFrame(TableIO.read_table(joinpath(split(@__FILE__, '#')[1] * ".assets", "test.jdf")); copycols=false)
            end""",

            """df_test_parquet = let
                import TableIO
                import DataFrames
                DataFrames.DataFrame(TableIO.read_table(joinpath(split(@__FILE__, '#')[1] * ".assets", "test.parquet")); copycols=false)
            end""",

            """df_test_arrow = let
                import TableIO
                import DataFrames
                DataFrames.DataFrame(TableIO.read_table(joinpath(split(@__FILE__, '#')[1] * ".assets", "test.arrow")); copycols=false)
            end""",

            """df_test_xlsx = let
                import TableIO
                import DataFrames
                # table_list = TableIO.list_tables(joinpath(split(@__FILE__, '#')[1] * ".assets", "test.xlsx")) # uncomment to get a list of all tables in this file
                DataFrames.DataFrame(TableIO.read_table(joinpath(split(@__FILE__, '#')[1] * ".assets", "test.xlsx"),
                    # "tablename", # define the table name here, if not defined the alphabetically first table is loaded
                    ); copycols=false)
            end""",

            """df_test_json = let
                import TableIO
                import DataFrames
                DataFrames.DataFrame(TableIO.read_table(joinpath(split(@__FILE__, '#')[1] * ".assets", "test.json")); copycols=true)
            end""",

            """df_test_zip = let
                import TableIO
                import DataFrames
                # table_list = TableIO.list_tables(joinpath(split(@__FILE__, '#')[1] * ".assets", "test.zip")) # uncomment to get a list of all tables in this file
                DataFrames.DataFrame(TableIO.read_table(joinpath(split(@__FILE__, '#')[1] * ".assets", "test.zip"),
                    # "tablename", # define the table name here, if not defined the alphabetically first table is loaded
                    ); copycols=false)
            end""",

            """df_test_sqlite = let
                import TableIO
                import DataFrames
                # table_list = TableIO.list_tables(joinpath(split(@__FILE__, '#')[1] * ".assets", "test.sqlite")) # uncomment to get a list of all tables in this file
                DataFrames.DataFrame(TableIO.read_table(joinpath(split(@__FILE__, '#')[1] * ".assets", "test.sqlite"),
                    # "tablename", # define the table name here, if not defined the alphabetically first table is loaded
                    ); copycols=false)
            end""",

            """df_test_hdf = let
                import TableIO
                import DataFrames
                # table_list = TableIO.list_tables(joinpath(split(@__FILE__, '#')[1] * ".assets", "test.hdf")) # uncomment to get a list of all tables in this file
                DataFrames.DataFrame(TableIO.read_table(joinpath(split(@__FILE__, '#')[1] * ".assets", "test.hdf"),
                    # "tablename", # define the table name here, if not defined the alphabetically first table is loaded
                    ); copycols=false)
            end""",

            """df_test_jld2 = let
                import TableIO
                import DataFrames
                # table_list = TableIO.list_tables(joinpath(split(@__FILE__, '#')[1] * ".assets", "test.jld2")) # uncomment to get a list of all tables in this file
                DataFrames.DataFrame(TableIO.read_table(joinpath(split(@__FILE__, '#')[1] * ".assets", "test.jld2"),
                    # "tablename", # define the table name here, if not defined the alphabetically first table is loaded
                    ); copycols=false)
            end"""

        ]

        for (i, fname) ∈ enumerate(fnames)
            sample_code = get_example_code("my_notebook.jl.assets", fname)
            @test sample_code == expected_sample_code[i]
            # println(sample_code)
        end

        @test is_extension_supported("csv") && !is_extension_supported("foo") && is_extension_supported("xlsx")

    end

end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

module GLWindow

using ModernGL
using GLAbstraction
using GLFW
using Reactive
using GeometryTypes
using ColorTypes
using FixedPointNumbers
using FixedSizeArrays
using Compat
using FileIO

import GLFW.Window
import GLFW.Monitor
import GLAbstraction.render
import GeometryTypes.widths

include("types.jl")

include("core.jl")
include("events.jl")
export pressed, dragged, clicked

include("callbacks.jl")
include("render.jl")
include("screen.jl")

export createwindow
export swapbuffers
export pollevents
export Screen
export UnicodeInput
export KeyPressed
export MouseClicked
export MouseMoved
export EnteredWindow
export WindowResized
export MouseDragged
export Scrolled
export Window
export leftclickdown
export Screen
export primarymonitorresolution
export renderloop
export render_frame
export screenshot
export screenbuffer
export zeroposition
export create_glcontext
export renderlist
export destroy!
export robj_from_camera

end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

importall Base

element_has_nodes(::Type{Val{:C3D4}}) = 4
element_has_type( ::Type{Val{:C3D4}}) = :Tet4

element_has_nodes(::Type{Val{:C3D8}}) = 8
element_has_type( ::Type{Val{:C3D8}}) = :Hex8

element_has_nodes(::Type{Val{:C3D10}}) = 10
element_has_type(::Type{Val{:C3D10}}) = :Tet10

element_has_nodes(::Type{Val{:C3D10}}) = 10
element_has_nodes(::Type{Val{:C3D20}}) = 20

element_has_nodes(::Type{Val{:C3D20E}}) = 20

element_has_nodes(::Type{Val{:S3}}) = 3
element_has_type( ::Type{Val{:S3}}) = :Tri3

element_has_nodes(::Type{Val{:STRI65}}) = 6
element_has_type(::Type{Val{:STRI65}}) = :Tri6

element_has_nodes(::Type{Val{:CPS4}}) = 4
element_has_type(::Type{Val{:CPS4}}) = :Quad4

"""
Checks for if line is a comment line or just empty
"""
function empty_or_comment_line{T<:AbstractString}(line::T)
    startswith(line, "**") || (length(line) == 0)
end

"""
Function for parsing nodes from the file
"""
function parse_section(model, lines, key::Symbol, idx_start,
    idx_end, ::Type{Val{:NODE}})

    info("Parsing nodes")
    ids = Integer[]
    definition = lines[idx_start]
    for line in lines[idx_start + 1: idx_end]
        if !(empty_or_comment_line(line))
            m = matchall(r"[-0-9.]+", line)
            node_id = parse(Int, m[1])
            coords = float(m[2:end])
            model["nodes"][node_id] = coords
        end
    end
    has_set_def = match(r"NSET=([\w\_\-]+)", definition)
    if has_set_def != nothing
      set_name = has_set_def[1]
      info("Creating nset $set_name")
      model["nsets"][set_name] = ids
  end
end

"""
Custon regex to find match from string. Index used if there are multiple matches
"""
function regex_match(regex_str, line, idx)
    return match(regex_str, line).captures[idx]
    println(eltype_sym)
end

"""
Simple iterator for comsuming element list. Depending
on the used element, connectivity nodes might be listed
in multiple lines, which is why iterator is used to handle
this problem.
"""
function consumeList(arr, start, stop)
    function _it()
        for i=start:stop
            produce(arr[i])
        end
    end
    Task(_it)
end

"""
Parse elements from input.
"""
function parse_section(model, lines, key, idx_start, idx_end, ::Type{Val{:ELEMENT}})
    #definition = uppercase(lines[idx_start])
    definition = lines[idx_start]
    element_type = regex_match(r"TYPE=([\w\-\_]+)", definition, 1)
    eltype_sym = Symbol(element_type)
    eltype_nodes = element_has_nodes(Val{eltype_sym})
    element_type = element_has_type(Val{eltype_sym})
    info("Parsing elements. Type: $(element_type)")
    list_iterator = consumeList(lines, idx_start+1, idx_end)
    ids = Integer[]
    line = consume(list_iterator)
    while line != nothing
        arr_num_as_str = matchall(r"[0-9]+", line)
        numbers = map(x-> parse(Int, x), arr_num_as_str)
        if !(empty_or_comment_line(line))
            id = numbers[1]
            push!(ids, id)
            connectivity = numbers[2:end]
            while length(connectivity) != eltype_nodes
                @assert length(connectivity) < eltype_nodes
                line = consume(list_iterator)
                arr_num_as_str = matchall(r"[0-9]+", line)
                numbers = map(x-> parse(Int, x), arr_num_as_str)
                push!(connectivity, numbers...)
            end
            model["elements"][id] = Dict(("type"=>element_type),
                                ("connectivity"=>connectivity))
        end
        line = consume(list_iterator)
    end
    has_set_def = match(r"ELSET=([\w\_\-]+)", definition)
    if has_set_def != nothing
      set_name = has_set_def[1]
      info("Creating elset $set_name")
      model["elsets"][set_name] = ids
  end
end

"""
Parsing Node- and ElementSets.
"""
function parse_section(model, lines, key, idx_start, idx_end, ::Union{Type{Val{:NSET}}, Type{Val{:ELSET}}})
    set_regex_string = Dict(:NSET  => r"NSET=([\w\-\_]+)",
                            :ELSET => r"ELSET=([\w\-\_]+)" )
    #definition = uppercase(lines[idx_start])
    # FIXME: do not uppercase set names
    definition = lines[idx_start]
    regex_string = set_regex_string[key]
    set_name = regex_match(regex_string, definition, 1)
    info("Creating $(lowercase(string(key))) $set_name")
    data = Integer[]
    if endswith(strip(definition), "GENERATE")
        line = lines[idx_start + 1]
        m = matchall(r"[0-9]+", line)
        first_id, last_id, step = map(x-> parse(Int, x), m)
        set_ids = collect(first_id:step:last_id)
        push!(data, set_ids...)
    else
        for line in lines[idx_start + 1: idx_end]
            if !(empty_or_comment_line(line))
                m = matchall(r"[0-9]+", line)
                set_ids = map(s -> parse(Int, s), m)
                push!(data, set_ids...)
            end
        end
    end
    selected_set = key == :NSET ? "nsets" : "elsets"
    model[selected_set][set_name] = data
end

"""
Parse SURFACE keyword
"""
function parse_section(model, lines, key, idx_start, idx_end, ::Type{Val{:SURFACE}})
    info("Parsing surface")
    #definition = uppercase(lines[idx_start])
    definition = lines[idx_start]
    has_set_def = match(r"TYPE=([\w\_\-]+),.*NAME=([\w\_\-]+)", definition)
    has_set_def != nothing || return
    set_type = Symbol(has_set_def[1])
    set_name = Symbol(has_set_def[2])
    data = Vector{Tuple{Int64, Symbol}}()
    for line in lines[idx_start + 1: idx_end]
        empty_or_comment_line(line) && continue
        m = match(r"(?P<element_id>\d+),.+(?P<element_side>S\d+).*", line)
        element_id = parse(Int, m[:element_id])
        element_side = Symbol(m[:element_side])
        push!(data, (element_id, element_side))
    end
    model["surface_types"][set_name] = set_type
    model["surfaces"][set_name] = data
    return
end

"""
Find lines, which contain keywords, for example "*NODE"
"""
function find_keywords(lines)
    indexes = Integer[]
    for (idx, line) in enumerate(lines)
        if startswith(line, "*") && !startswith(line, "**")
            push!(indexes, idx)
        end
    end
    push!(indexes, length(lines) + 1)
    return indexes
end

"""
Main function for parsing Abaqus input file.
"""
function parse_abaqus(fid::IOStream)
    lines = readlines(fid)
    keyword_indexes = find_keywords(lines)
    idx_start = keyword_indexes[1]
    keyword_sym::Symbol = :none
    parser::Function = x->()
    model = Dict{AbstractString, Any}()
    model["nodes"] = Dict{Int64, Vector{Float64}}()
    model["nsets"] = Dict{AbstractString, Vector{Int64}}()
    model["elsets"] = Dict{AbstractString, Vector{Int64}}()
    model["elements"] = Dict{Integer, Any}()
    model["surfaces"] = Dict{Symbol, Vector{Tuple{Int64, Symbol}}}()
    model["surface_types"] = Dict{Symbol, Symbol}()
    for idx_end in keyword_indexes[2:end]
        keyword_line = uppercase(lines[idx_start])
        keyword = regex_match(r"\s*([\w ]+)", keyword_line, 1)
        k_sym = Symbol(keyword)
        args = Tuple{Dict, Vector{Int}, Symbol, Int, Int, Type{Val{k_sym}}}
        if method_exists(parse_section, args)
            parse_section(model, lines, k_sym, idx_start, idx_end-1, Val{k_sym})
#       else
#           warn("Unknown section: $(keyword)")
        end
        idx_start = idx_end
    end
    return model
end

function abaqus_read_mesh(fn)
    model = open(parse_abaqus, fn)
    mesh = Mesh()
    mesh.nodes = model["nodes"]
    for (nset_name, node_ids) in model["nsets"]
        mesh.node_sets[Symbol(nset_name)] = Set(node_ids)
    end
    for (elid, eldata) in model["elements"]
        eltype = eldata["type"]
        elcon = eldata["connectivity"]
        mesh.elements[elid] = elcon
        mesh.element_types[elid] = eltype
    end
    for (elset_name, element_ids) in model["elsets"]
        mesh.element_sets[Symbol(elset_name)] = Set(element_ids)
    end
    mesh.surfaces = model["surfaces"]
    mesh.surface_types = model["surface_types"]
    return mesh
end


text:
extractpoint(a::Point) = (a.data...,)
function pointview(::Type{Point}, A::AbstractArray)
    @assert size(A,1) == 3 "First dimension must be of size 3."
    mappedarray(Point, extractpoint,
                A[expand2indices(A,(1,))...],
                A[expand2indices(A,(2,))...],
                A[expand2indices(A,(3,))...])
end


text:
using Test

# Test that the macro call `m` throws an error exception during pre-compilation
macro test_macro_throws(errortype, m)
    # See https://discourse.julialang.org/t/test-throws-with-macros-after-pr-23533/5878
    :(@test_throws $errortype try @eval $m catch err; throw(err.error) end)
end

using JuMP
const MOIT = MOI.Test

function bridged_mock(mock_optimize!::Function...;
                      model = JuMP._MOIModel{Float64}())
    mock = MOI.Utilities.MockOptimizer(model)
    bridged = MOI.Bridges.full_bridge_optimizer(mock, Float64)
    MOI.Utilities.set_mock_optimize!(mock, mock_optimize!...)
    return bridged
end


text:
"""
A `Problem` stores the high-level description of a motion planning problem.

This is the main structure of the package.
The idea is for users to start with an empty problem for a given robot.
Then, in order to specify the requirements for the motion they want to compute, they can add different types of constraints, one by one.
See [Creating a Problem](@ref) and [Defining Constraints](@ref) for detailed examples.

Eventually, the problem becomes a complex description with all the requirements.
Once the description is finalized, it can be optimized.
This happens in one of the `TORA.solve_with_*` functions, which does the heavy-lifting of converting a `Problem` (the high-level description) into an [NLP](https://en.wikipedia.org/wiki/Nonlinear_programming) problem that can be solved using off-the-shelf solvers, such as [Ipopt](https://github.com/coin-or/Ipopt) and [Knitro](https://www.artelys.com/solvers/knitro/).

See also: [`solve_with_ipopt`](@ref), [`solve_with_knitro`](@ref)
"""
mutable struct Problem
    num_knots::Int  # Total number of knots
    dt::Float64  # Time step between two knots

    fixed_q::Dict{Int64,Array{Float64}}  # Fixed joint positions
    fixed_v::Dict{Int64,Array{Float64}}  # Fixed joint velocities
    fixed_τ::Dict{Int64,Array{Float64}}  # Fixed joint torques

    ee_pos::Dict{Int64,Point{3,Float64}}  # End-effector target positions

    # Jacobian data (e.g., for sparse Jacobian of dynamics with AD)
    jacdata_fwd_dyn::JacobianData
    jacdata_inv_dyn::JacobianData
    jacdata_ee_position::JacobianData

    @doc """
        Problem(robot, num_knots, dt)

    Create an empty `Problem`.

    # Arguments
    - `robot::Robot`: the robot associated with this `Problem`.
    - `num_knots::Int`: the number of knots used to represent the motion.
    - `dt::Float64`: the time step duration between each two consecutive knots.

    See also: [`show_problem_info`](@ref)
    """
    function Problem(robot::Robot, num_knots::Int, dt::Float64)
        @assert num_knots > 0
        @assert dt > 0

        input, output = rand(2 * (robot.n_q + robot.n_v) + robot.n_τ), rand(robot.n_q + robot.n_v)
        jacdata_fwd_dyn = JacobianData((out, x) -> forward_dynamics_defects!(out, robot, x, dt), output, input)

        input, output = rand(2 * (robot.n_q + robot.n_v) + robot.n_τ), rand(robot.n_q + robot.n_τ)
        jacdata_inv_dyn = JacobianData((out, x) -> inverse_dynamics_defects!(out, robot, x, dt), output, input)

        input, output = rand(robot.n_q), rand(3)
        jacdata_ee_position = JacobianData((out, x) -> ee_position!(out, robot, x), output, input)

        new(
            num_knots, dt,
            Dict{Int64,Array{Float64}}(),
            Dict{Int64,Array{Float64}}(),
            Dict{Int64,Array{Float64}}(),
            Dict{Int64,Point{3,Float64}}(),
            jacdata_fwd_dyn,
            jacdata_inv_dyn,
            jacdata_ee_position
        )
    end
end

"""
    fix_joint_positions!(problem, robot, knot, q)

Fix the joint positions of the `robot` at a specific `knot`.

See also: [`fix_joint_velocities!`](@ref), [`fix_joint_torques!`](@ref)
"""
function fix_joint_positions!(problem::Problem, robot::Robot, knot, q)
    @assert 1 <= knot <= problem.num_knots
    @assert all(robot.q_lo .<= q .<= robot.q_hi)
    problem.fixed_q[knot] = q
    return problem
end

"""
    fix_joint_velocities!(problem, robot, knot, v)

Fix the joint velocities of the `robot` at a specific `knot`.

See also: [`fix_joint_positions!`](@ref), [`fix_joint_torques!`](@ref)
"""
function fix_joint_velocities!(problem::Problem, robot::Robot, knot, v)
    @assert 1 <= knot <= problem.num_knots
    @assert all(robot.v_lo .<= v .<= robot.v_hi)
    problem.fixed_v[knot] = v
    return problem
end

"""
    fix_joint_torques!(problem, robot, knot, τ)

Fix the joint torques of the `robot` at a specific `knot`.

See also: [`fix_joint_positions!`](@ref), [`fix_joint_velocities!`](@ref)
"""
function fix_joint_torques!(problem::Problem, robot::Robot, knot, τ)
    @assert 1 <= knot <= problem.num_knots - 1
    @assert all(robot.τ_lo .<= τ .<= robot.τ_hi)
    problem.fixed_τ[knot] = τ
    return problem
end

"""
    constrain_ee_position!(problem, knot, position)

Fix the end-efector positions of the `robot` at a specific `knot`.

See also: [`fix_joint_positions!`](@ref), [`fix_joint_velocities!`](@ref), [`fix_joint_torques!`](@ref)
"""
function constrain_ee_position!(problem::Problem, knot, position)
    @assert 1 <= knot <= problem.num_knots
    @assert length(position) == 3
    problem.ee_pos[knot] = position
    return problem
end

"""
    show_problem_info(problem)

Output a summary of the problem, including the number of knots with constraints.
"""
function show_problem_info(problem::Problem)
    t = (problem.num_knots - 1) * problem.dt

    println("Motion duration ...................................... $(t) seconds")
    println("Number of knots with constrained joint positions ..... $(length(problem.fixed_q))")
    println("Number of knots with constrained joint velocities .... $(length(problem.fixed_v))")
    println("Number of knots with constrained joint torques ....... $(length(problem.fixed_τ))")
    println("Number of knots with constrained ee position ......... $(length(problem.ee_pos))")
end

export
    Problem,
    fix_joint_positions!,
    fix_joint_velocities!,
    fix_joint_torques!,
    constrain_ee_position!,
    show_problem_info


text:
"""
    vel_copy!(vel::NamedTuple, u, v, w, g::AbstractGrid{TX, TY, TZ}) where {TX, TY, TZ}
Copy external velocities into `PlanktonModel`
"""
function vel_copy!(vel::NamedTuple, u, v, w, g::AbstractGrid{TX, TY, TZ}) where {TX, TY, TZ}
    copy_interior_u!(vel.u.data, u, g, TX())
    copy_interior_v!(vel.v.data, v, g, TY())
    copy_interior_w!(vel.w.data, w, g, TZ())

    fill_halo_vel!(vel, g)
end

@inline function copy_interior!(c, t, g::AbstractGrid)
    copyto!(view(c, g.Hx+1:g.Hx+g.Nx, g.Hy+1:g.Hy+g.Ny, g.Hz+1:g.Hz+g.Nz), t)
end
@inline function copy_interior_u!(c, t, g::AbstractGrid, ::Periodic)
    copyto!(view(c, g.Hx+1:g.Hx+g.Nx, g.Hy+1:g.Hy+g.Ny, g.Hz+1:g.Hz+g.Nz), t)
end
@inline function copy_interior_v!(c, t, g::AbstractGrid, ::Periodic)
    copyto!(view(c, g.Hx+1:g.Hx+g.Nx, g.Hy+1:g.Hy+g.Ny, g.Hz+1:g.Hz+g.Nz), t)
end
@inline function copy_interior_w!(c, t, g::AbstractGrid, ::Periodic)
    copyto!(view(c, g.Hx+1:g.Hx+g.Nx, g.Hy+1:g.Hy+g.Ny, g.Hz+1:g.Hz+g.Nz), t)
end

@inline function copy_interior_u!(c, t, g::AbstractGrid, ::Bounded)
    copyto!(view(c, g.Hx+1:g.Hx+g.Nx+1, g.Hy+1:g.Hy+g.Ny, g.Hz+1:g.Hz+g.Nz), t)
end
@inline function copy_interior_v!(c, t, g::AbstractGrid, ::Bounded)
    copyto!(view(c, g.Hx+1:g.Hx+g.Nx, g.Hy+1:g.Hy+g.Ny+1, g.Hz+1:g.Hz+g.Nz), t)
end
@inline function copy_interior_w!(c, t, g::AbstractGrid, ::Bounded)
    copyto!(view(c, g.Hx+1:g.Hx+g.Nx, g.Hy+1:g.Hy+g.Ny, g.Hz+1:g.Hz+g.Nz+1), t)
end

function validate_temp(sim::PlanktonSimulation, g::AbstractGrid{TX, TY, TZ}) where {TX, TY, TZ}
    temp_size = (g.Nx, g.Ny, g.Nz)
    validation = true

    if size(sim.input.temp)[1:3] ≠ temp_size
        validation = false
        throw(ArgumentError("Dimension mismatch: the size of temperature must be $(temp_size)."))
    end

    @assert sim.input.ΔT_temp % sim.ΔT == 0.0

    if size(sim.input.temp)[4] < floor(Int,sim.iterations*sim.ΔT/sim.input.ΔT_temp)
        throw(ArgumentError("Temperature provided cannot cover the time period of (iterations+1)*ΔT ($(sim.iterations*sim.ΔT) seconds) with ΔT_temp = $(sim.input.ΔT_temp)."))
    end

    return validation
end

function validate_PARF(sim::PlanktonSimulation, g::AbstractGrid{TX, TY, TZ}) where {TX, TY, TZ}
    PARF_size = (g.Nx, g.Ny)
    validation = true

    if size(sim.input.PARF)[1:2] ≠ PARF_size
        validation = false
        throw(ArgumentError("Dimension mismatch: the size of PARF must be $(PARF_size)."))
    end

    @assert sim.input.ΔT_PAR % sim.ΔT == 0.0

    if size(sim.input.PARF)[3] < floor(Int,sim.iterations*sim.ΔT/sim.input.ΔT_PAR)
        throw(ArgumentError("Surface PAR provided cannot cover the time period of (iterations+1)*ΔT ($(sim.iterations*sim.ΔT) seconds) with ΔT_PAR = $(sim.input.ΔT_PAR)."))
    end

    return validation
end

function validate_velocity(sim::PlanktonSimulation, g::AbstractGrid{TX, TY, TZ}) where {TX, TY, TZ}
    if sim.input.vels ≠ (;)
        u_size = (g.Nx, g.Ny, g.Nz)
        v_size = (g.Nx, g.Ny, g.Nz)
        w_size = (g.Nx, g.Ny, g.Nz)
        validation = true

        if isa(TX(), Bounded)
            u_size = (g.Nx+1, g.Ny, g.Nz)
        end

        if isa(TY(), Bounded)
            v_size = (g.Nx, g.Ny+1, g.Nz)
        end

        if isa(TZ(), Bounded)
            w_size = (g.Nx, g.Ny, g.Nz+1)
        end

        if size(sim.input.vels.u)[1:3] ≠ u_size
            validation = false
            throw(ArgumentError("Dimension mismatch: the size of u must be $(u_size), for $(TX) topology."))
        end

        if size(sim.input.vels.v)[1:3] ≠ v_size
            validation = false
            throw(ArgumentError("Dimension mismatch: the size of v must be $(v_size), for $(TY) topology."))
        end

        if size(sim.input.vels.w)[1:3] ≠ w_size
            validation = false
            throw(ArgumentError("Dimension mismatch: the size of w must be $(w_size), for $(TZ) topology."))
        end

        @assert sim.input.ΔT_vel % sim.ΔT == 0.0

        if size(sim.input.vels.u)[4] < floor(Int,sim.iterations*sim.ΔT/sim.input.ΔT_vel)
            throw(ArgumentError("Velocities provided cannot cover the time period of (iterations+1)*ΔT ($(sim.iterations*sim.ΔT) seconds) with ΔT_vel = $(sim.input.ΔT_vel)."))
        end
        return validation
    else
        return true
    end
end

function set_vels_fields!(sim::PlanktonSimulation, uv, vv, wv)
    sim.input.vels = (u = uv, v = vv, w = wv)
    validate_velocity(sim, sim.model.grid)
    return nothing
end
function set_PARF_fields!(sim::PlanktonSimulation, PARF)
    sim.input.PARF = PARF
    validate_PARF(sim, sim.model.grid)
    return nothing
end
function set_temp_fields!(sim::PlanktonSimulation, temp)
    sim.input.temp = temp
    validate_temp(sim, sim.model.grid)
    return nothing
end


text:
using Nextion
using Test

@testset "pbar" begin

    nexSerial = NexSerial("/dev/ttyUSB0")

    nexProgressBar = NexProgressBar(nexSerial, Name("j0"), cid=ComponentID(1))

    send(nexSerial, "page pg_pbar")

    sleep(1)

    nexProgressBar.value = 30

    @test_throws ErrorException nexProgressBar.value = 105  # should raise error because value must be in 0-100
    @test_throws ErrorException nexProgressBar.value = -1  # should raise error because value must be in 0-100

    sleep(1)

    close(nexSerial)

    @test true

end


text:
module test_lordeckcodes_unittest1

using Test
using LoRDeckCodes

# EncodeDecodeRecommendedDecks
function read_remove_bom(path::String)
    bytes = read(path)
    bytes[1:3] == [0xef, 0xbb, 0xbf] && deleteat!(bytes, 1:3)
    String(bytes)
end

file_string = replace(read_remove_bom(normpath(@__DIR__, "DeckCodesTestData.txt")), "\r" => "")
deck_blocks = split(file_string, "\n\n")
for deck_block in deck_blocks
    isempty(deck_block) && continue
    deckcode_and_card_list = split(deck_block, "\n")
    deckcode = first(deckcode_and_card_list)
    card_list = filter(x -> !isempty(x),deckcode_and_card_list[2:end])
    deck = Deck(CardCodeAndCount.(card_list))
    @test DeckEncoder.encode_deck(deck) == deckcode
end

# SmallDeck
deck = Deck([CardCodeAndCount("01DE002", 1)])
@test Deck(DeckEncoder.encode_deck(deck)) == deck

# LargeDeck
deck = Deck([CardCodeAndCount("01DE002", 3),
             CardCodeAndCount("01DE003", 3),
             CardCodeAndCount("01DE004", 3),
             CardCodeAndCount("01DE005", 3),
             CardCodeAndCount("01DE006", 3),
             CardCodeAndCount("01DE007", 3),
             CardCodeAndCount("01DE008", 3),
             CardCodeAndCount("01DE009", 3),
             CardCodeAndCount("01DE010", 3),
             CardCodeAndCount("01DE011", 3),
             CardCodeAndCount("01DE012", 3),
             CardCodeAndCount("01DE013", 3),
             CardCodeAndCount("01DE014", 3),
             CardCodeAndCount("01DE015", 3),
             CardCodeAndCount("01DE016", 3),
             CardCodeAndCount("01DE017", 3),
             CardCodeAndCount("01DE018", 3),
             CardCodeAndCount("01DE019", 3),
             CardCodeAndCount("01DE020", 3),
             CardCodeAndCount("01DE021", 3)])
deckcode = DeckEncoder.encode_deck(deck)
@test Deck(DeckEncoder.encode_deck(deck)) == deck

# DeckWithCountsMoreThan3Small
deck = Deck([CardCodeAndCount("01DE002", 4)])
@test Deck(DeckEncoder.encode_deck(deck)) == deck

# DeckWithCountsMoreThan3Large
deck = Deck([CardCodeAndCount("01DE002", 3),
             CardCodeAndCount("01DE003", 3),
             CardCodeAndCount("01DE004", 3),
             CardCodeAndCount("01DE005", 3),
             CardCodeAndCount("01DE006", 4),
             CardCodeAndCount("01DE007", 5),
             CardCodeAndCount("01DE008", 6),
             CardCodeAndCount("01DE009", 7),
             CardCodeAndCount("01DE010", 8),
             CardCodeAndCount("01DE011", 9),
             CardCodeAndCount("01DE012", 3),
             CardCodeAndCount("01DE013", 3),
             CardCodeAndCount("01DE014", 3),
             CardCodeAndCount("01DE015", 3),
             CardCodeAndCount("01DE016", 3),
             CardCodeAndCount("01DE017", 3),
             CardCodeAndCount("01DE018", 3),
             CardCodeAndCount("01DE019", 3),
             CardCodeAndCount("01DE020", 3),
             CardCodeAndCount("01DE021", 3)])
deckcode = DeckEncoder.encode_deck(deck)
@test Deck(deckcode) == Deck(deck.cards)

# SingleCardCodeAndCount40Times
deck = Deck([CardCodeAndCount("01DE002", 40)])
@test Deck(DeckEncoder.encode_deck(deck)) == deck

# WorstCaseLength
deck = Deck([CardCodeAndCount("01DE002", 4),
             CardCodeAndCount("01DE003", 4),
             CardCodeAndCount("01DE004", 4),
             CardCodeAndCount("01DE005", 4),
             CardCodeAndCount("01DE006", 4),
             CardCodeAndCount("01DE007", 5),
             CardCodeAndCount("01DE008", 6),
             CardCodeAndCount("01DE009", 7),
             CardCodeAndCount("01DE010", 8),
             CardCodeAndCount("01DE011", 9),
             CardCodeAndCount("01DE012", 4),
             CardCodeAndCount("01DE013", 4),
             CardCodeAndCount("01DE014", 4),
             CardCodeAndCount("01DE015", 4),
             CardCodeAndCount("01DE016", 4),
             CardCodeAndCount("01DE017", 4),
             CardCodeAndCount("01DE018", 4),
             CardCodeAndCount("01DE019", 4),
             CardCodeAndCount("01DE020", 4),
             CardCodeAndCount("01DE021", 4)])
@test Deck(DeckEncoder.encode_deck(deck)) == deck

# OrderShouldNotMatter1
deck1 = Deck([CardCodeAndCount("01DE002", 1), CardCodeAndCount("01DE003", 2), CardCodeAndCount("02DE003", 3)])
deck2 = Deck([CardCodeAndCount("01DE003", 2), CardCodeAndCount("02DE003", 3), CardCodeAndCount("01DE002", 1)])
@test DeckEncoder.encode_deck(deck1) == DeckEncoder.encode_deck(deck2)
deck3 = Deck([CardCodeAndCount("01DE002", 4), CardCodeAndCount("01DE003", 2), CardCodeAndCount("02DE003", 3)])
deck4 = Deck([CardCodeAndCount("01DE003", 2), CardCodeAndCount("02DE003", 3), CardCodeAndCount("01DE002", 4)])
@test DeckEncoder.encode_deck(deck3) == DeckEncoder.encode_deck(deck4)

# OrderShouldNotMatter2
deck1 = Deck([CardCodeAndCount("01DE002", 4), CardCodeAndCount("01DE003", 2), CardCodeAndCount("02DE003", 3), CardCodeAndCount("01DE004", 5)])
deck2 = Deck([CardCodeAndCount("01DE004", 5), CardCodeAndCount("01DE003", 2), CardCodeAndCount("02DE003", 3), CardCodeAndCount("01DE002", 4)])
@test DeckEncoder.encode_deck(deck1) == DeckEncoder.encode_deck(deck2)

# BadCardCodeAndCountCodes
@test  DeckEncoder.isvalid(CardCodeAndCount("01DE002", 1))
@test !DeckEncoder.isvalid(CardCodeAndCount("01DE02", 1))
@test !DeckEncoder.isvalid(CardCodeAndCount("01XX002", 1))
using LoRDeckCodes: ArgumentException
@test_throws ArgumentException DeckEncoder.encode_deck(Deck([CardCodeAndCount("01DE02", 1)]))

# BadCount
@test !DeckEncoder.isvalid(CardCodeAndCount("01DE002", 0))
@test !DeckEncoder.isvalid(CardCodeAndCount("01DE002", -1))

# GarbageDecoding
badEncodingNotBase32 = "I'm no card code!"
badEncoding32 = "ABCDEFG"
badEncodingEmpty = ""
@test_throws ArgumentException Deck(badEncodingNotBase32)
@test_throws ArgumentException Deck(badEncoding32)
@test_throws ArgumentException Deck(badEncodingEmpty)

@test  DeckEncoder.isvalid(deckcode="CEAAECABAIDASDASDISC2OIIAECBGGY4FAWTINZ3AICACAQXDUPCWBABAQGSOKRM")
@test !DeckEncoder.isvalid(deckcode=badEncodingNotBase32)
@test !DeckEncoder.isvalid(deckcode=badEncoding32)
@test !DeckEncoder.isvalid(deckcode=badEncodingEmpty)
@test  DeckEncoder.isvalid(cardcode="01DE002")
@test !DeckEncoder.isvalid(cardcode="01DE02")
@test !DeckEncoder.isvalid(cardcode="01XX002")

end # module test_lordeckcodes_unittest1


text:
# FUNCTIONS

function is_prime_6k1(n)
    if n <= 3
        return n > 1
    elseif n % 2 == 0 || n % 3 == 0
        return false
    end

    i = 5
    while i * i <= n
        if n % i == 0 || n % (i + 2) == 0
            return false
        end
        i += 6
    end

    return true
end

# SCRIPT

function main()
    i = 1
    n = 1
    while i != 10001
        n += 2
        i += is_prime_6k1(n)
    end

    println(n)
end

main()

max_stars_repo_name: shanedrabing/project-euler

max_issues_repo_name: shanedrabing/project-euler

max_forks_repo_name: shanedrabing/project-euler


text:
using EvidentialDeepLearning, Documenter, Literate

ENV["GKS_ENCODING"] = "utf-8"

DocMeta.setdocmeta!(EvidentialDeepLearning, :DocTestSetup, :(using EvidentialDeepLearning); recursive=true)

inputdir = joinpath(@__DIR__, "..", "examples")
outputdir = joinpath(@__DIR__, "src", "examples")
mkpath(outputdir)
mkpath(joinpath(outputdir, "examples"))
for filename in readdir(inputdir)
    inpath = joinpath(inputdir, filename)
    cp(inpath, joinpath(outputdir, "examples", filename), force = true)
    Literate.markdown(inpath, outputdir; documenter = true)
end
cp(joinpath(@__DIR__, "..", "README.md"), joinpath(@__DIR__, "src", "index.md"), force=true)

makedocs(;
    modules=[EvidentialDeepLearning],
    authors="Jan Weidner <jw3126@gmail.com> and contributors",
    repo="https://github.com/jw3126/EvidentialDeepLearning.jl/blob/{commit}{path}#{line}",
    sitename="EvidentialDeepLearning.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jw3126.github.io/EvidentialDeepLearning.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Tutorials" => ["Basic Regression" => "examples/basic_regression.md",],
        # "Explanation" => [],
        "Reference" => ["Docstrings" => "docstrings.md"],
        # "How-to guides" => [],
    ],
    strict=false,
)

deploydocs(;
    repo="github.com/jw3126/EvidentialDeepLearning.jl",
    push_preview=true,
)


text:
const pi = 3.14159265358979323846264338327   #Pi
const grav = 9.8                               #Gravitational acceleration (m / s^2)
const cp = 1004.0                             #Specific heat of dry air at constant pressure
const cv = 717.0                              #Specific heat of dry air at constant volume
const rd = 287.0                              #Dry air constant for equation of state (P=rho*rd*T)
const p0 = 1.e5                              #Standard pressure at the surface in Pascals
const C0 = 27.5629410929725921310572974482   #Constant to translate potential temperature into pressure (P=C0*(rho*theta)**gamma)
const gamma = 1.40027894002789400278940027894   #gamma=cp/Rd
#Define domain and stability-related constants
const xlen = 2.e4    #Length of the domain in the x-direction (meters)
const zlen = 1.e4    #Length of the domain in the z-direction (meters)
const hv_beta = 0.25     #How strong to diffuse the solution: hv_beta \in [0:1]
const cfl = 1.50    #"Courant, Friedrichs, Lewy" number (for numerical stability)
const max_speed = 450        #Assumed maximum wave speed during the simulation (speed of sound + speed of wind) (meter / sec)
const hs = 2          #"Halo" size: number of cells beyond the MPI tasks's domain needed for a full "stencil" of information for reconstruction
const sten_size = 4          #Size of the stencil used for interpolation

#Parameters for indexing and flags
const NUM_VARS = 4           #Number of fluid state variables
const ID_DENS = 1           #index for density ("rho")
const ID_UMOM = 2           #index for momentum in the x-direction ("rho * u")
const ID_WMOM = 3           #index for momentum in the z-direction ("rho * w")
const ID_RHOT = 4           #index for density * potential temperature ("rho * theta")
const DIR_X = 1              #Integer constant to express that this operation is in the x-direction
const DIR_Z = 2              #Integer constant to express that this operation is in the z-direction


text:
vl2 = DummyDisplayable{MIME"application/vnd.vegalite.v2+json"}("""
{
  "data": {
    "values": [
      {"a": "A","b": 28}, {"a": "B","b": 55}, {"a": "C","b": 43}
    ]
  },
  "mark": "bar",
  "encoding": {
    "x": {"field": "a", "type": "ordinal"},
    "y": {"field": "b", "type": "quantitative"}
  }
}
""")

vl3 = DummyDisplayable{MIME"application/vnd.vegalite.v3+json"}("""
{
  "data": {
    "values": [
      {"a": "A","b": 28}, {"a": "B","b": 55}, {"a": "C","b": 43}
    ]
  },
  "mark": "bar",
  "encoding": {
    "x": {"field": "a", "type": "ordinal"},
    "y": {"field": "b", "type": "quantitative"}
  }
}
""")

vg3 = DummyDisplayable{MIME"application/vnd.vega.v3+json"}("""
{
  "width": 400,
  "height": 200,
  "padding": 5,

  "data": [
    {
      "name": "table",
      "values": [
        {"category": "A", "amount": 28},
        {"category": "B", "amount": 55},
        {"category": "C", "amount": 43},
        {"category": "D", "amount": 91},
        {"category": "E", "amount": 81},
        {"category": "F", "amount": 53},
        {"category": "G", "amount": 19},
        {"category": "H", "amount": 87}
      ]
    }
  ],

  "signals": [
    {
      "name": "tooltip",
      "value": {},
      "on": [
        {"events": "rect:mouseover", "update": "datum"},
        {"events": "rect:mouseout",  "update": "{}"}
      ]
    }
  ],

  "scales": [
    {
      "name": "xscale",
      "type": "band",
      "domain": {"data": "table", "field": "category"},
      "range": "width",
      "padding": 0.05,
      "round": true
    },
    {
      "name": "yscale",
      "domain": {"data": "table", "field": "amount"},
      "nice": true,
      "range": "height"
    }
  ],

  "axes": [
    { "orient": "bottom", "scale": "xscale" },
    { "orient": "left", "scale": "yscale" }
  ],

  "marks": [
    {
      "type": "rect",
      "from": {"data":"table"},
      "encode": {
        "enter": {
          "x": {"scale": "xscale", "field": "category"},
          "width": {"scale": "xscale", "band": 1},
          "y": {"scale": "yscale", "field": "amount"},
          "y2": {"scale": "yscale", "value": 0}
        },
        "update": {
          "fill": {"value": "steelblue"}
        },
        "hover": {
          "fill": {"value": "red"}
        }
      }
    },
    {
      "type": "text",
      "encode": {
        "enter": {
          "align": {"value": "center"},
          "baseline": {"value": "bottom"},
          "fill": {"value": "#333"}
        },
        "update": {
          "x": {"scale": "xscale", "signal": "tooltip.category", "band": 0.5},
          "y": {"scale": "yscale", "signal": "tooltip.amount", "offset": -2},
          "text": {"signal": "tooltip.amount"},
          "fillOpacity": [
            {"test": "datum === tooltip", "value": 0},
            {"value": 1}
          ]
        }
      }
    }
  ]
}
""")

vg4 = DummyDisplayable{MIME"application/vnd.vega.v4+json"}("""
{
  "width": 400,
  "height": 200,
  "padding": 5,

  "data": [
    {
      "name": "table",
      "values": [
        {"category": "A", "amount": 28},
        {"category": "B", "amount": 55},
        {"category": "C", "amount": 43},
        {"category": "D", "amount": 91},
        {"category": "E", "amount": 81},
        {"category": "F", "amount": 53},
        {"category": "G", "amount": 19},
        {"category": "H", "amount": 87}
      ]
    }
  ],

  "signals": [
    {
      "name": "tooltip",
      "value": {},
      "on": [
        {"events": "rect:mouseover", "update": "datum"},
        {"events": "rect:mouseout",  "update": "{}"}
      ]
    }
  ],

  "scales": [
    {
      "name": "xscale",
      "type": "band",
      "domain": {"data": "table", "field": "category"},
      "range": "width",
      "padding": 0.05,
      "round": true
    },
    {
      "name": "yscale",
      "domain": {"data": "table", "field": "amount"},
      "nice": true,
      "range": "height"
    }
  ],

  "axes": [
    { "orient": "bottom", "scale": "xscale" },
    { "orient": "left", "scale": "yscale" }
  ],

  "marks": [
    {
      "type": "rect",
      "from": {"data":"table"},
      "encode": {
        "enter": {
          "x": {"scale": "xscale", "field": "category"},
          "width": {"scale": "xscale", "band": 1},
          "y": {"scale": "yscale", "field": "amount"},
          "y2": {"scale": "yscale", "value": 0}
        },
        "update": {
          "fill": {"value": "steelblue"}
        },
        "hover": {
          "fill": {"value": "red"}
        }
      }
    },
    {
      "type": "text",
      "encode": {
        "enter": {
          "align": {"value": "center"},
          "baseline": {"value": "bottom"},
          "fill": {"value": "#333"}
        },
        "update": {
          "x": {"scale": "xscale", "signal": "tooltip.category", "band": 0.5},
          "y": {"scale": "yscale", "signal": "tooltip.amount", "offset": -2},
          "text": {"signal": "tooltip.amount"},
          "fillOpacity": [
            {"test": "datum === tooltip", "value": 0},
            {"value": 1}
          ]
        }
      }
    }
  ]
}
""")

vg5 = DummyDisplayable{MIME"application/vnd.vega.v5+json"}("""
{
  "width": 400,
  "height": 200,
  "padding": 5,

  "data": [
    {
      "name": "table",
      "values": [
        {"category": "A", "amount": 28},
        {"category": "B", "amount": 55},
        {"category": "C", "amount": 43},
        {"category": "D", "amount": 91},
        {"category": "E", "amount": 81},
        {"category": "F", "amount": 53},
        {"category": "G", "amount": 19},
        {"category": "H", "amount": 87}
      ]
    }
  ],

  "signals": [
    {
      "name": "tooltip",
      "value": {},
      "on": [
        {"events": "rect:mouseover", "update": "datum"},
        {"events": "rect:mouseout",  "update": "{}"}
      ]
    }
  ],

  "scales": [
    {
      "name": "xscale",
      "type": "band",
      "domain": {"data": "table", "field": "category"},
      "range": "width",
      "padding": 0.05,
      "round": true
    },
    {
      "name": "yscale",
      "domain": {"data": "table", "field": "amount"},
      "nice": true,
      "range": "height"
    }
  ],

  "axes": [
    { "orient": "bottom", "scale": "xscale" },
    { "orient": "left", "scale": "yscale" }
  ],

  "marks": [
    {
      "type": "rect",
      "from": {"data":"table"},
      "encode": {
        "enter": {
          "x": {"scale": "xscale", "field": "category"},
          "width": {"scale": "xscale", "band": 1},
          "y": {"scale": "yscale", "field": "amount"},
          "y2": {"scale": "yscale", "value": 0}
        },
        "update": {
          "fill": {"value": "steelblue"}
        },
        "hover": {
          "fill": {"value": "red"}
        }
      }
    },
    {
      "type": "text",
      "encode": {
        "enter": {
          "align": {"value": "center"},
          "baseline": {"value": "bottom"},
          "fill": {"value": "#333"}
        },
        "update": {
          "x": {"scale": "xscale", "signal": "tooltip.category", "band": 0.5},
          "y": {"scale": "yscale", "signal": "tooltip.amount", "offset": -2},
          "text": {"signal": "tooltip.amount"},
          "fillOpacity": [
            {"test": "datum === tooltip", "value": 0},
            {"value": 1}
          ]
        }
      }
    }
  ]
}
""")

# based on https://github.com/mathiasbynens/small
dummy_png_bytes = UInt8[
    0x89, 0x50, 0x4e, 0x47, 0x0d, 0x0a, 0x1a, 0x0a, 0x00, 0x00, 0x00, 0x0d,
    0x49, 0x48, 0x44, 0x52, 0x00, 0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x01,
    0x08, 0x06, 0x00, 0x00, 0x00, 0x1f, 0x15, 0xc4, 0x89, 0x00, 0x00, 0x00,
    0x0a, 0x49, 0x44, 0x41, 0x54, 0x78, 0x9c, 0x63, 0x00, 0x01, 0x00, 0x00,
    0x05, 0x00, 0x01,
]

struct DummyVegaLitePNGHybrid end
Base.show(io::IO, ::MIME"application/vnd.vegalite.v3+json", ::DummyVegaLitePNGHybrid) =
    print(io, vl3)
Base.show(io::IO, ::MIME"image/png", ::DummyVegaLitePNGHybrid) =
    write(io, dummy_png_bytes)

vl3png = DummyVegaLitePNGHybrid()