text:
"""
    half_measurement_index(lx, m, n)

Returns the index corresponding to the half_measurement of lattice
size `lx`. That means the index in the vector When `n > m` for all `m` and
`n`.

"""
# more readable version : div(lx*(lx-1),2) - div((lx-m+1)*(lx-m),2) + (n-m)
half_measurement_index(lx::Int, m::Int, n::Int) =
    1 <= m < n <= lx ? div((2*lx-m)*(m-1), 2) + (n-m) : error("1 <= $m < $n <= $lx")

"""
    randisometry(d1, d2)

Returns a uniformly chosen random isometry matrix of size `d1` × `d2`
"""
randisometry(T,rng::AbstractRNG, d1, d2) =
    d1 >= d2 ? Matrix(qr!(randn(rng, T, d1, d2)).Q) : Matrix(lq!(randn(rng, T, d1, d2)).Q)
randisometry(rng, d1, d2) = randisometry(Float64, rng, d1, d2)

"""
    svdtrunc(A [;maxdim, tol])

Performs the singular value decomposition (SVD) and truncated the
singular values (SV) accoding to the tolerance `tol` or `maxdim`
specified, whichever results in a lower number of SVs.

Return value is the tuple of truncated `U`, `S`, `Vt`. Where `U`, and
`Vt` are left and right isometries and `S` is a Diagonal matrix
containing the SVs.

"""
function svdtrunc(A::Matrix{T};
                  maxdim::Int=size(A, 2),
                  tol::Float64=1.e-12) where{T<:Number}
    fact = svd(A, full=false)
    n = min(maxdim, sum(fact.S .> fact.S[1]*tol))
    n == 0 && error("The largest singular value is smaller than tol! $m < $tol")

    fact.U[:, 1:n], Diagonal(fact.S[1:n]), fact.Vt[1:n, :]
end

"""
    picklargests(Ss, maxdim)

for a vector of descending sorted vectors `Ss`, pick the `maxdim`
largest numbers and discard the rest.

"""
##TODO: make this better and find the best version!
function picklargests(Ss::Vector{Vector{Float64}},
                      maxdim::Int=200, tol::Float64=1.e-14)
    #vs = Vector{Vector{Float64}}(Float64[], length(Ss))
    vs = [Float64[] for i=1:length(Ss)]

    pointers = ones(Int, length(Ss))
    c = [Ss[i][pointers[i]] for i in eachindex(Ss)]

    for n=1:min(maxdim, sum([length(s) for s in Ss]))

        m, index = findmax(c)
        if (m < tol)
            n==1 && error("The largest singular value is smaller than tol! $m < $tol")
            break
        end
        push!(vs[index], Ss[index][pointers[index]])
        if pointers[index] == length(Ss[index])

            c[index] = 0.0
        else
            pointers[index] += 1
            c[index] = Ss[index][pointers[index]]
        end
    end
    return vs
end

"""
    truncatedsizes(S [; maxdim, tol])

for a vector of float vectors `S` of singular values (or basically
it is just assumed that each sector is sorted descendingly), returns
the vector of new sizes when the singular values are to be truncated
to the `maxdim` largest numbers or when `tol` is reached whichever is
lower.

"""
function truncatedsizes(S::Vector{Vector{Float64}};
                        maxdim::Int=200,
                        tol::Float64=1.e-14)

    pointers = ones(Int, length(S))

    c = [S[i][pointers[i]] for i in eachindex(S)]
    for n=1:min(maxdim, sum([length(s) for s in S]))

        m, index = findmax(c)
        if (m < tol)
            n==1 && error("The largest singular value is smaller than tol! $m < $tol")
            break
        end
        if pointers[index] >= size(S[index], 1)
            pointers[index] += 1
            c[index] = 0.0
        else
            pointers[index] += 1
            c[index] = S[index][pointers[index]]
        end
    end
    pointers .- 1
end

######NOTE
###TODO: this needs to be tested!
function svdtrunc(A::AbstractSymTensor;
                  maxdim::Int=maximum(size(A)),
                  tol::Float64=1.e-12)
    T = eltype(A)
    S = vtype(A)

    n_sects = length(A.blocks)
    sects_u  = Vector{Sector{S, 2}}(undef, n_sects)
    sects_s  = Vector{Sector{S, 2}}(undef, n_sects)
    sects_v = Vector{Sector{S, 2}}(undef, n_sects)

    blks_u  = Vector{Matrix{T}}(undef, n_sects)
    blks_s  = Vector{Vector{Float64}}(undef, n_sects)
    blks_v = Vector{Matrix{T}}(undef, n_sects)

    semits = collect(onlysemitokens(A.blocks))
    i=1
    fact = svd(zeros(T, 1,1), full=false)
    for (sect, blk) in A.blocks
        c1, c2 = sect
        sects_u[i] = Sector{S}(c1, c1)
        sects_s[i] = Sector{S}(c1, c2)
        sects_v[i] = Sector{S}(c2, c2)

        try
            fact = svd(blk, full=false)
        catch e
            @warn "error $e"
            fact = svd(blk, full=false, alg=LinearAlgebra.QRIteration())
        end
        blks_u[i] = fact.U
        blks_s[i] = fact.S
        blks_v[i] = fact.Vt
        i += 1
    end

    ns = truncatedsizes(blks_s, maxdim=maxdim, tol=tol)
    indices = findall(x->x>0, ns)

    sects_u  = sects_u[indices]
    sects_s  = sects_s[indices]
    sects_v = sects_v[indices]

    blks_u = [blks_u[index][:, 1:ns[index]] for index in indices]
    blks_s = [blks_s[index][1:ns[index]] for index in indices]
    blks_v = [blks_v[index][1:ns[index], :] for index in indices]

    middledims = ns[indices]

    n_sectors =length(middledims)
    Vl = VectorSpace{S}([sects_s[i][1]=>middledims[i] for i=1:n_sectors],
                        isdual(A.space[1]))
    Vr = VectorSpace{S}([sects_s[i][2]=>middledims[i] for i=1:n_sectors],
                        isdual(A.space[2]))

    u  = SymMatrix{S,T}(zero(S), (A.space[1], dual(Vl)),
                        SortedDict(zip(sects_u,blks_u)))
    s  = SymDiagonal{S,Float64}(A.charge, (Vl, Vr),
                                SortedDict([sects_s[i]=>Diagonal(blks_s[i]) for i in 1:n_sectors]))
    v = SymMatrix{S,T}(zero(S), (dual(Vr), A.space[2]),
                       SortedDict(zip(sects_v,blks_v)))

    u, s, v
end

"""
    entorpy(P [; alpha])

Calculate the entropy of a vector of numbers `P`. The numbers are
assumed to form a probability distribution which means they add up to
one.

 If `alpha=1` (the default) it calculates the usual Shannon
(Von-Neumann) entropy and if `alpha > 1` calculates the Renyi entorpy.

"""
function entropy(P::Vector{T};
                 alpha::Int=1) where {T<:Number}

    alpha < 1 && error("Renyi entropy only defined for positive alpha")
    any(P .< 0.0) && error("Negative value in P vector ", P)
    sum(P) ≈ 1 || println("P vector doesn't add up to one : $(sum(P))")
    Pnz = filter(x->x>0.0, P)
    if alpha == 1
        return -sum(Pnz .* log.(Pnz))
    else
        return log(sum(Pnz.^alpha))/(1-alpha)
    end
end

function _realwithcheck(a::Number, tol::Float64=1.e-12)
    imag(a) > tol && error("Not real! Im = $(imag(a))")
    real(a)
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
# This file (`LibWasm.jl` not `prologue.jl`) is automatically generated using
# Clang.jl and should not be edited manually. Take a look at the `gen/` folder
# if something is to be changed.

using Pkg.Artifacts
using Pkg.BinaryPlatforms

tripletnolibc(platform) = replace(triplet(platform), "-gnu" => "")
wasmtime_folder_name(platform) =
    "wasmtime-v$release_version-$(tripletnolibc(platform))-c-api"

function get_libwasmtime_location()
    artifact_info = artifact_meta("libwasmtime", joinpath(@__DIR__, "..", "Artifacts.toml"))
    artifact_info === nothing && return nothing

    parent_path = artifact_path(Base.SHA1(artifact_info["git-tree-sha1"]))
    child_folder = readdir(parent_path)[1]
    return joinpath(
        parent_path,
        child_folder,
        "lib/libwasmtime"
    )
end

const libwasmtime_env_key = "LIBWASMTIME_LOCATION"
const libwasmtime = haskey(ENV, libwasmtime_env_key) ?
    ENV[libwasmtime_env_key] : get_libwasmtime_location()

ext: jl


text:
module contactForceLaw_ballWithBall

using Modia3D
using Modia3D.StaticArrays
import Modia3D.ModiaMath

vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)    # material of SolidFileMesh
vmat11 = Modia3D.Material(color="Red" , transparency=0.5)
vmat2 = deepcopy(vmat1)                                           # material of convex decomposition of SolidFileMesh
vmat2.transparency = 0.7

cmat = "Steel"

@assembly BallWithBall begin
  world = Modia3D.Object3D(visualizeFrame=true)

  sphereMoving     = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidSphere(0.5) , "Steel", vmat11; contactMaterial = cmat); r=[0.0, 0.0, 0.0],  fixed=false ) # , R=ModiaMath.rot2(-pi/2) )

  sphere     = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidSphere(1.0) , "Steel", vmat1; contactMaterial = cmat); r=[-3.0, 0.0, 0.0], fixed=true) #, R=ModiaMath.rot2(-pi/3) )
end

gravField = Modia3D.UniformGravityField(g=9.81, n=[-1,0,0])
threeD = BallWithBall(sceneOptions=Modia3D.SceneOptions(gravityField=gravField,visualizeFrames=true, defaultFrameLength=0.7,nz_max = 100, enableContactDetection=true, visualizeContactPoints=true, visualizeSupportPoints=true))

model = Modia3D.SimulationModel( threeD )
result = ModiaMath.simulate!(model; stopTime=10.0, tolerance=1e-8,interval=0.001, log=false)

println("... success of contactForceLaw_ballWithBall.jl!")
end


text:
# println("\nmain.jl: DEBUT DES INCLUDES à $(ms()).")
@ms include("args.jl");
@ms include("utils/log.jl")
@ms include("utils/console.jl")
@ms include("plane.jl")
@ms include("instance.jl")
@ms include("processing/instance_read_alp.jl")
@ms include("utils/array.jl")
@ms include("processing/file.jl")
@ms include("solution.jl")
@ms include("processing/solution_read.jl")

@ms using Dates
@ms using Statistics

@ms using Random: randperm
@ms using Distributions: NegativeBinomial
@ms using StatsBase: sample, ProbabilityWeights
@ms include("solvers/neighbour_operators.jl")

@ms include("utils/PlotUtils.jl")
@ms using .PlotUtils: plot_save_costs

# Les autres fichiers sont inclue seulement si nécessaire
# Dans les solveurs spécifiques ou dans le fichier usings.jl

# TODO :
# - faire les macros @p (pour imprimer) et @s (pour to_string générique)
#
# using DataStructures # risque de pollution du scope
# import DataStructures: OrderedDict

function main()
    @ms args = Args.parse_commandline(ARGS)
    ln1("main() BEGIN") # ln1 n'est utilisable que si level est connu
    lg1() && Args.show_args()

    # Le fichier suivant déclare différents solveurs PL/PLNE (par des instructions
    # `using ...`) selon les options de la ligne de commande.
    # Les instructions `using ...` **doivent** être dans des fichiers séparés et
    # non pas dans des fonctions
    #
    @time include("$APPDIR/src/usings.jl")

    # date1= now() # en secondes entières
    time1= time() # secondes précision microsecondes
    action = args[:action]
    if args[:loop]
        if action == :descent
            include("$APPDIR/src/loop_descent.jl")
        elseif action == :annealing
            include("$APPDIR/src/loop_annealing.jl")
        elseif action == :greedy
            include("$APPDIR/src/loop_greedy.jl")
        end
    else
        if action == :validate
            include("$APPDIR/src/main_validate.jl")
        elseif action == :timing
            include("$APPDIR/src/main_timing.jl")
        elseif action == :carlo
            include("$APPDIR/src/main_carlo.jl")
        elseif action == :explore
            include("$APPDIR/src/main_explore.jl")
        elseif action == :descent
            include("$APPDIR/src/main_descent.jl")
        elseif action == :greedy
            include("$APPDIR/src/main_greedy.jl")
        elseif action == :annealing
            include("$APPDIR/src/main_annealing.jl")
        elseif action == :dmip
            include("$APPDIR/src/main_dmip.jl")
        elseif action == :stats
            include("$APPDIR/src/main_stats.jl")
        elseif action == :test
            include("$APPDIR/src/main_test.jl")
        elseif action == :none
            println("Aucune action indiquée")
            println(Args.get_syntaxe())
            exit(1)
        else
            println("Erreur : action $action non implémentée (dans main.jl)")
            println(Args.get_syntaxe())
            exit(1)
        end
    end

    # if costs !== nothing
    #     inst_name, _ = splitext(basename(args[:infile]))
    #     save_path = "$APPDIR/_tmp/figures/$(inst_name)_$(action)_$(args[:itermax])"
    #     plot_save_costs(costs, save_path, plot=args[:plot], save=args[:plot])
    # else
    #     println("No costs recorded.")
    # end

    # heure de fin du traitement
    time2= time()
    sec = round((time() - time1), digits=3) # on veut limiter la précision à la ms
    ln1("Durée totale du main 1000*(time2-time1) : $(sec)s")
    ln1("main() END")
end


text:
"""
    Tensor{k,label}

Abstract base for any tensor of rank `k`.
"""
abstract type Tensor{k,label} end

(::Type{T})(k::Int) where {T<:Tensor} = T{k}()

LinearAlgebra.rank(::Tensor{k}) where k = k

"""
    system(::Tensor)

A general tensor acts on the full system, i.e. all coordinates.
"""
system(::Type{<:Tensor}) = FullSystem()

couples(a::SpinOrbital{<:Orbital}, ::Type{T}, b::SpinOrbital{<:Orbital}) where {T<:Tensor} =
    isequal(other_quantum_numbers(system(T), a, b)...)

couples(a::SpinOrbital{<:RelativisticOrbital}, ::Type{T}, b::SpinOrbital{<:RelativisticOrbital}) where {T<:Tensor} =
    isequal(first.(other_quantum_numbers(system(T), a, b))...)

# This is only true for SO(3)
components(::Tensor{k}) where k = -k:k

function Base.show(io::IO, ::Tensor{k,label}) where {k,label}
    write(io,to_boldface(label)*"̂")
    write(io,"⁽",to_superscript(k),"⁾")
end

# * Tensor components

"""
    TensorComponent(tensor, q)

Represents the `q`th component of a `tensor`; `abs(q) ≤ rank(tensor)`.
"""
struct TensorComponent{T<:Tensor}
    tensor::T
    q::Int
    function TensorComponent(tensor::T, q::Int) where {T<:Tensor}
        k = rank(tensor)
        q ∈ components(tensor) ||
            throw(ArgumentError("Tensor component $(q) not possible for tensor $(tensor) of rank $k"))
        new{T}(tensor, q)
    end
end

function Base.show(io::IO, Tq::TensorComponent)
    show(io, Tq.tensor)
    write(io, to_subscript(Tq.q))
end

Base.parent(Tᵏq::TensorComponent) = Tᵏq.tensor
component(Tᵏq::TensorComponent) = Tᵏq.q

system(::TensorComponent{T}) where T = system(T)

Base.hash(Tq::TensorComponent, h::UInt) = hash(Tq.tensor, hash(Tq.q, h))

# * Tensor products

@doc raw"""
    TensorProduct{K}(T, U)

A tensor of rank `K` formed from the product of the tensors `T` and
`U`, according to

```math
\begin{equation}
\tag{V3.1.20}
\tensor{X}^{(K)}_Q \equiv
\{\tensor{T}^{(k_1)}\tensor{U}^{(k_2)}\}^{(K)}_Q \defd
\tensor{T}^{(k_1)}_{q_1}
\tensor{U}^{(k_2)}_{q_2}
C_{k_1q_1k_2q_2}^{KQ}
\end{equation}
```
"""
struct TensorProduct{K,A<:Tensor,B<:Tensor} <: Tensor{K,'X'}
    T::A
    U::B
end
TensorProduct(K::Int, T::A, U::B) where {A<:Tensor, B<:Tensor} =
    TensorProduct{K,A,B}(T, U)

function Base.show(io::IO, X::TensorProduct{K}) where K
    write(io,"{")
    show(io, X.T)
    show(io, X.U)
    write(io,"}⁽",to_superscript(K),"⁾")
end

system(::TensorProduct{<:Any,A,B}) where {A,B} = (system(A), system(B))

# ** Scalar product

@doc raw"""
    TensorScalarProduct(T, U)

A tensor of rank 0 (and thus implicitly 0 projection) formed from the
product of the tensors `T` and `U` (which have to have the same rank),
according to

```math
\begin{equation}
\tag{V3.1.30,35}
(\tensor{T}^{(k)} \cdot
\tensor{U}^{(k)}) \defd
(-)^k\angroot{k}
\{\tensor{T}^{(k)}
\tensor{U}^{(k)}\}^{(0)}_0 \equiv
(-)^q
\tensor{T}^{(k)}_{q}
\tensor{U}^{(k)}_{-q}
\end{equation}
```
"""
struct TensorScalarProduct{A<:Tensor,B<:Tensor} <: Tensor{0,'X'}
    T::A
    U::B
    function TensorScalarProduct(T::A, U::B) where {A<:Tensor,B<:Tensor}
        rank(T) == rank(U) ||
            throw(ArgumentError("Cannot form scalar product of tensors of differing ranks"))
        new{A,B}(T,U)
    end
end

function Base.show(io::IO, X::TensorScalarProduct)
    write(io,"(")
    show(io, X.T)
    write(io,"⋅")
    show(io, X.U)
    write(io,")")
end

"""
    dot(T::Tensor, U::Tensor)

Form the scalar product of the two tensors `T` and `U`, which need to
have the same rank.

# Examples

```jldoctest
julia> SphericalTensor(4)⋅SphericalTensor(4)
(𝐂̂⁽⁴⁾⋅𝐂̂⁽⁴⁾)
```
"""
LinearAlgebra.dot(T::Tensor, U::Tensor) =
    TensorScalarProduct(T, U)

system(::TensorScalarProduct{A,B}) where {A,B} = (system(A), system(B))

# * Linear combination of tensors

"""
    LinearCombinationTensor

Represents a linear combination of tensor components.
"""
const LinearCombinationTensor{T<:Tensor,N<:Number} = LinearCombination{<:TensorComponent{<:T},N}

@linearly_combinable TensorComponent

export Tensor, TensorComponent,
    TensorProduct, TensorScalarProduct,
    system


text:
using My-Julia-ML-Package
using Documenter

DocMeta.setdocmeta!(My-Julia-ML-Package, :DocTestSetup, :(using My-Julia-ML-Package); recursive=true)

makedocs(;
    modules=[My-Julia-ML-Package],
    authors="Sandesh Katakam",
    repo="https://github.com/sandeshkatakam/My-Julia-ML-Package.jl/blob/{commit}{path}#{line}",
    sitename="My-Julia-ML-Package.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://sandeshkatakam.github.io/My-Julia-ML-Package.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/sandeshkatakam/My-Julia-ML-Package.jl",
    devbranch="master",
)


text:
# 별칭 부여
IntOrBool = Union{Int, Bool}

function main()
    num::IntOrBool = 10 # Int값 할당
    println(num)

    num = true # Bool도 할당 가능
    println(num)
end

main()

max_stars_repo_name: myyrakle/Julia-Example

max_issues_repo_name: myyrakle/Julia-Example

max_forks_repo_name: myyrakle/Julia-Example


text:
using CSV
using DataFrames
using Plots
using JLD2
using DiscretisedFluidQueues
import Distributions
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/plot_utils.jl")
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/error_metrics.jl")

@halfwidth_plot_defaults()

plotlyjs()
function p()
ks_data = CSV.read((@__DIR__)*"/../data/ks_errors.csv",DataFrame)

plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot]
markerstyles_vec = [:cross,:diamond,:circle,:dot]
colours = [1;2;4]
for (c,col) in enumerate(names(ks_data))
    plot!(log10.(1:2:21),ks_data[:,col],
        label=col,linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2,color=colours[c])
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="KS error")
# plot!(legend=(0.7,0.7))
plot!()
@add_lines!(ks_data,("Unif","QBDRAP"),@__DIR__)

savefig((@__DIR__)*"/ks_error_formatted.pdf")
end
p()

function p()
l2_pdf_data = CSV.read((@__DIR__)*"/../data/l2_pdf_errors.csv",DataFrame)
colours = [1;2;4]

plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot]
markerstyles_vec = [:cross,:diamond,:circle,:cross]
for (c,col) in enumerate(names(l2_pdf_data))
    plot!(log10.(1:2:21),l2_pdf_data[:,col],
        label=col,linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2,color=colours[c])
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="L² error - PDF")
# plot!(legend=(0.7,0.8))
plot!()
@add_lines!(l2_pdf_data,("Unif","QBDRAP"),@__DIR__)

savefig((@__DIR__)*"/l2_pdf_error_formatted.pdf")
end
p()

function p()
l1_cdf_data = CSV.read((@__DIR__)*"/../data/l1_cdf_errors.csv",DataFrame)
@halfwidth_plot_defaults()
colours = [1;2;4]
plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot]
markerstyles_vec = [:cross,:diamond,:circle,:cross]
for (c,col) in enumerate(names(l1_cdf_data))
    plot!(log10.(1:2:21),l1_cdf_data[:,col],
        label=col,linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2,color=colours[c])
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="L¹ error - CDF")
# plot!(legend=(0.7,0.7))
plot!()
@add_lines!(l1_cdf_data,("Unif","QBDRAP"),@__DIR__)

savefig((@__DIR__)*"/l1_cdf_error_formatted.pdf")
end
p()

# @fullwidth_plot_defaults()
# dir = @__DIR__()
# os = 1:2:7
# plot(layout=(4,length(os)))
# linestyles_vec = [:solid,:solid,:solid,:solid]
# markerstyles_vec = [:cross,:diamond,:circle]
# xs = range(eps(),10-eps(),length=200)
# for (c,o) in enumerate(os)
#     model = BoundedFluidQueue(fill(0.0,2,2),[1.0; 0.0],10.0)

#     mesh = DGMesh(0.0:10.0,o)
#     dq = DiscretisedFluidQueue(model,mesh)
#     jldopen(dir*"/../../../coeffs_matrix.jld2") do f
#         coeffs = f["coeff_matrix"][fun_number,c,1]
#         dt = SFMDistribution(coeffs,dq)
#         plot!(xs,pdf(dt).(xs,1),subplot=c,label=(c==1 ? "DG" : false),
#             linestyle=linestyles_vec[1],ylims = (-0.1,1.1),  yticks=false,
#             # marker=markerstyles_vec[1],
#             linewidth=2,xticks=false)
#     end
#     plot!(subplot=c, yticks=(c==1 ? (0.1:0.4:1.2) : false),
#         ylabel=(c==1 ? "DG" : false), grid=false,
#         title = string("Dimension ",o))

#     mesh = DGMesh(0.0:(1.0./o):10.0,1)
#     dq = DiscretisedFluidQueue(model,mesh)
#     jldopen(dir*"/../../../coeffs_matrix.jld2") do f
#         coeffs = f["coeff_matrix"][fun_number,c,2]
#         dt = SFMDistribution(coeffs,dq)
#         plot!(xs,pdf(dt).(xs,1),subplot=c+4,label=(c==1 ? "unif." : false),
#             linestyle=linestyles_vec[2],ylims = (-0.1,1.1),  yticks=false,
#             # marker=markerstyles_vec[2],
#             linewidth=2,xticks=false)
#     end
#     (c==1)&&plot!(subplot=c+4, yticks=(c==1 ? (0.1:0.4:1.2) : false),
#         ylabel=(c==1 ? "unif." : false), grid=false)

#     mesh = DGMesh(0.0:10.0,o)
#     dq = DiscretisedFluidQueue(model,mesh)
#     jldopen(dir*"/../../../coeffs_matrix.jld2") do f
#         coeffs = f["coeff_matrix"][fun_number,c,3]
#         dt = SFMDistribution(coeffs,dq)
#         plot!(xs,pdf(dt).(xs,1),subplot=c+8,label=(c==1 ? "DG (limiter)" : false),
#             linestyle=linestyles_vec[3],ylims = (-0.1,1.1),  yticks=false,
#             # marker=markerstyles_vec[3],
#             linewidth=2,xticks=false)
#     end
#     (c==1)&&plot!(subplot=c+8, yticks=(c==1 ? (0.1:0.4:1.2) : false),
#         ylims = (-0.1,1.1),
#         ylabel=(c==1 ? "DG-limit" : false), grid=false)

#     mesh = FRAPMesh(0.0:10.0,o)
#     dq = DiscretisedFluidQueue(model,mesh)
#     jldopen(dir*"/../../../coeffs_matrix.jld2") do f
#         coeffs = f["coeff_matrix"][fun_number,c,4]
#         dt = SFMDistribution(coeffs,dq)
#         plot!(xs,pdf(dt).(xs,1),subplot=c+12,label=(c==1 ? "QBD-RAP" : false),
#             linestyle=linestyles_vec[3],ylims = (-0.1,1.1), yticks=false,
#             # marker=markerstyles_vec[3],
#             linewidth=2,xticks=0:2:10, grid=false)
#     end
#     (c==1)&&plot!(subplot=c+12, yticks=(c==1 ? (0.1:0.4:1.2) : false),
#         ylims = (-0.1,1.1),
#         ylabel=(c==1 ? "QBD-RAP" : false), grid=false)
#     # plot!(xs,x->2*x*(x<0.5) + Float64(x>=0.5),subplot=c,label=(c==1 ? "truth" : false),
#     #     linestyle=:dot,
#     #     # marker=markerstyles_vec[3],
#     #     color=:grey,
#     #     linewidth=2)

#     # plot!(subplot=c, yticks=(c==1 ? (0.1:0.3:1.2) : false),
#     #     ylims = (-0.1,1.1), xticks = 0.0:3:10,
#     #     ylabel=(c==1 ? "probability" : false), grid=false,
#     #     title = string("Dimension ",o),
#     #     legend=:outertopright)
#     plot!(legend=false)
# end
# plot!()
# savefig((@__DIR__)*"/pdf_formatted.pdf")

# @halfwidth_plot_defaults()
# L1_cell_probs_errors = DataFrame(DG = Float64[], Order_1 = Float64[], DG_limiter = Float64[], QBDRAP = Float64[])
# i=1; t=4
# F4(z) = exp(-1*t)*(exp(z)-1)*(z<=t) + (z>t)*(1-exp(-1*t))
# truth = F4.(1:10)-F4.(0:9)
# truth_pm = zeros(3)
# for (c_o,o) in enumerate(1:2:21)
#     row = zeros(4)
#     for (c_m,m) in enumerate(names(L1_cell_probs_errors))
#         jldopen(dir*"/../../../coeffs_matrix.jld2") do f
#             model = BoundedFluidQueue(fill(0.0,2,2),[1.0; 0.0],10.0)

#             mesh = (c_m<4) ? DGMesh(0.0:10.0,o) : FRAPMesh(0.0:10.0,o)
#             dq = DiscretisedFluidQueue(model,mesh)

#             coeffs = f["coeff_matrix"][fun_number,c_o,c_m]
#             dt = SFMDistribution(coeffs,dq)

#             x_vals = range(0.5,9.5;length=10)

#             approx = sum(cell_probs(dt).(x_vals,(1:2)'),dims=2)
#             approx_pm = [coeffs[1:N₋(model)];coeffs[end-N₊(model)+1:end]]

#             (o==3)&&display(approx)
#             row[c_m] = log10(sum(abs.(approx-truth)) + sum(abs.(approx_pm-truth_pm)))
#         end
#     end
#     push!(L1_cell_probs_errors,row)
# end

# file = (@__DIR__)*"/../data/meshs_l1_cell_"*"func_count_$(fun_number)"
# CSV.write(file*".csv",L1_cell_probs_errors)
# q = plot()
# linestyles_vec = [:solid,:dash,:dashdot, :dot]
# markerstyles_vec = [:cross,:diamond,:circle,:dot]
# for (c,reconstruction) in enumerate(names(L1_cell_probs_errors))
#     plot!(log10.(1:2:21),L1_cell_probs_errors[:,reconstruction],label=reconstruction,
#         linestyle=linestyles_vec[c],
#         marker=markerstyles_vec[c],
#         linewidth=2)
# end
# plot!(xlabel="Dimension", ylabel="log₁₀ Error",
#     title="Error between cell masses",
#     legend=:outerbottomright)
# savefig((@__DIR__)*"/L1_cell_probs.pdf")


text:
export iauEe00
"""
The equation of the equinoxes, compatible with IAU 2000 resolutions,
given the nutation in longitude and the mean obliquity.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  canonical model.

Given:
   date1,date2  double    TT as a 2-part Julian Date (Note 1)
   epsa         double    mean obliquity (Note 2)
   dpsi         double    nutation in longitude (Note 3)

Returned (function value):
               double    equation of the equinoxes (Note 4)

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. The obliquity, in radians, is mean of date.

   3. The result, which is in radians, operates in the following sense:

      Greenwich apparent ST = GMST + equation of the equinoxes

   4. The result is compatible with the IAU 2000 resolutions.  For
   further details, see IERS Conventions 2003 and Capitaine et al.
   (2002).

Called:
   iauEect00    equation of the equinoxes complementary terms

References:

   Capitaine, N., Wallace, P.T. and McCarthy, D.D., "Expressions to
   implement the IAU 2000 definition of UT1", Astronomy &
   Astrophysics, 406, 1135-1149 (2003)

   McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
   IERS Technical Note No. 32, BKG (2004)

This revision:  2008 May 16

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# double iauEe00(double date1, double date2, double epsa, double dpsi)

function iauEe00(date1::Real, date2::Real, epsa::Real, dpsi::Real)
   return ccall((:iauEe00, libsofa_c), Cdouble,
               (Cdouble, Cdouble, Cdouble, Cdouble),
               convert(Float64, date1), convert(Float64, date2),
               convert(Float64, epsa), convert(Float64, dpsi))
end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

# convert binary into HDF5 data
using HDF5

base_dir = "./"

if length(ARGS) > 0
  base_dir = ARGS[1]
end

datasets = [("train", ["data_batch_$i.bin" for i in 1:5]),
            ("test", ["test_batch.bin"])]

const width      = 32
const height     = 32
const channels   = 3
const batch_size = 10000

mean_model = zeros(Float32, width, height, channels, 1)

for (key, sources) in datasets
  h5open("$base_dir/$key.hdf5", "w") do h5
    dset_data = d_create(h5, "data", datatype(Float32),
        dataspace(width, height, channels, batch_size * length(sources)))
    dset_label = d_create(h5, "label", datatype(Float32),
        dataspace(1, batch_size * length(sources)))

    for n = 1:length(sources)
      open("$base_dir/cifar-10-batches-bin/$(sources[n])") do f
        println("Processing $(sources[n])...")
        mat = readbytes(f, (1 + width*height*channels) * batch_size)
        mat = reshape(mat, 1+width*height*channels, batch_size)

        # random shuffle within batch
        rp = randperm(batch_size)

        label = convert(Array{Float32},mat[1, rp])
        # If I divide by 256 as in the MNIST example, then
        # training on the giving DNN gives me random
        # performance: objective function not changing,
        # and test performance is always 10%...
        # The same results could be observed when
        # running Caffe, as our HDF5 dataset is
        # compatible with Caffe.
        img = convert(Array{Float32},mat[2:end, rp]) / 256.0
        img = reshape(img, width, height, channels, batch_size)

        if key == "train"
          # only accumulate mean from the training data
          global mean_model
          mean_model = (batch_size*mean_model + sum(img, 4)) / (n*batch_size)
        end

        index = (n-1)*batch_size+1:n*batch_size
        dset_data[:,:,:,index] = img
        dset_label[:,index] = label
      end
    end

    # but apply mean subtraction for both training and testing data
    println("Subtracting the mean...")
    for n = 1:length(sources)
      index = (n-1)*batch_size+1:n*batch_size
      dset_data[:,:,:,index] = dset_data[:,:,:,index] .- mean_model
    end
  end
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

if is_windows()
    const ERROR_ENVVAR_NOT_FOUND = UInt32(203)

    _getenvlen(var::Vector{UInt16}) = ccall(:GetEnvironmentVariableW,stdcall,UInt32,(Ptr{UInt16},Ptr{UInt16},UInt32),var,C_NULL,0)
    _hasenv(s::Vector{UInt16}) = _getenvlen(s) != 0 || Libc.GetLastError() != ERROR_ENVVAR_NOT_FOUND
    _hasenv(s::AbstractString) = _hasenv(cwstring(s))

    function access_env(onError::Function, str::AbstractString)
        var = cwstring(str)
        len = _getenvlen(var)
        if len == 0
            return Libc.GetLastError() != ERROR_ENVVAR_NOT_FOUND ? "" : onError(str)
        end
        val = zeros(UInt16,len)
        ret = ccall(:GetEnvironmentVariableW,stdcall,UInt32,(Ptr{UInt16},Ptr{UInt16},UInt32),var,val,len)
        if (ret == 0 && len != 1) || ret != len-1 || val[end] != 0
            error(string("getenv: ", str, ' ', len, "-1 != ", ret, ": ", Libc.FormatMessage()))
        end
        pop!(val) # NUL
        return transcode(String, val)
    end

    function _setenv(svar::AbstractString, sval::AbstractString, overwrite::Bool=true)
        var = cwstring(svar)
        val = cwstring(sval)
        if overwrite || !_hasenv(var)
            ret = ccall(:SetEnvironmentVariableW,stdcall,Int32,(Ptr{UInt16},Ptr{UInt16}),var,val)
            systemerror(:setenv, ret == 0)
        end
    end

    function _unsetenv(svar::AbstractString)
        var = cwstring(svar)
        ret = ccall(:SetEnvironmentVariableW,stdcall,Int32,(Ptr{UInt16},Ptr{UInt16}),var,C_NULL)
        systemerror(:setenv, ret == 0)
    end
else # !windows
    _getenv(var::AbstractString) = ccall(:getenv, Cstring, (Cstring,), var)
    _hasenv(s::AbstractString) = _getenv(s) != C_NULL

    function access_env(onError::Function, var::AbstractString)
        val = _getenv(var)
        val == C_NULL ? onError(var) : unsafe_string(val)
    end

    function _setenv(var::AbstractString, val::AbstractString, overwrite::Bool=true)
        ret = ccall(:setenv, Int32, (Cstring,Cstring,Int32), var, val, overwrite)
        systemerror(:setenv, ret != 0)
    end

    function _unsetenv(var::AbstractString)
        ret = ccall(:unsetenv, Int32, (Cstring,), var)
        systemerror(:unsetenv, ret != 0)
    end
end # os test

## ENV: hash interface ##

"""
    EnvHash() -> EnvHash

A singleton of this type provides a hash table interface to environment variables.
"""
mutable struct EnvHash <: Associative{String,String}; end

"""
    ENV

Reference to the singleton `EnvHash`, providing a dictionary interface to system environment
variables.
"""
const ENV = EnvHash()

similar(::EnvHash) = Dict{String,String}()

getindex(::EnvHash, k::AbstractString) = access_env(k->throw(KeyError(k)), k)
get(::EnvHash, k::AbstractString, def) = access_env(k->def, k)
in(k::AbstractString, ::KeyIterator{EnvHash}) = _hasenv(k)
pop!(::EnvHash, k::AbstractString) = (v = ENV[k]; _unsetenv(k); v)
pop!(::EnvHash, k::AbstractString, def) = haskey(ENV,k) ? pop!(ENV,k) : def
delete!(::EnvHash, k::AbstractString) = (_unsetenv(k); ENV)
setindex!(::EnvHash, v, k::AbstractString) = _setenv(k,string(v))
push!(::EnvHash, k::AbstractString, v) = setindex!(ENV, v, k)

if is_windows()
    start(hash::EnvHash) = (pos = ccall(:GetEnvironmentStringsW,stdcall,Ptr{UInt16},()); (pos,pos))
    function done(hash::EnvHash, block::Tuple{Ptr{UInt16},Ptr{UInt16}})
        if unsafe_load(block[1]) == 0
            ccall(:FreeEnvironmentStringsW, stdcall, Int32, (Ptr{UInt16},), block[2])
            return true
        end
        return false
    end
    function next(hash::EnvHash, block::Tuple{Ptr{UInt16},Ptr{UInt16}})
        pos = block[1]
        blk = block[2]
        len = ccall(:wcslen, UInt, (Ptr{UInt16},), pos)
        buf = Array{UInt16}(len)
        unsafe_copy!(pointer(buf), pos, len)
        env = transcode(String, buf)
        m = match(r"^(=?[^=]+)=(.*)$"s, env)
        if m === nothing
            error("malformed environment entry: $env")
        end
        return (Pair{String,String}(m.captures[1], m.captures[2]), (pos+(len+1)*2, blk))
    end
else # !windows
    start(::EnvHash) = 0
    done(::EnvHash, i) = (ccall(:jl_environ, Any, (Int32,), i) === nothing)

    function next(::EnvHash, i)
        env = ccall(:jl_environ, Any, (Int32,), i)
        if env === nothing
            throw(BoundsError())
        end
        env = env::String
        m = match(r"^(.*?)=(.*)$"s, env)
        if m === nothing
            error("malformed environment entry: $env")
        end
        return (Pair{String,String}(m.captures[1], m.captures[2]), i+1)
    end
end # os-test

#TODO: Make these more efficent
function length(::EnvHash)
    i = 0
    for (k,v) in ENV
        i += 1
    end
    return i
end

function show(io::IO, ::EnvHash)
    for (k,v) = ENV
        println(io, "$k=$v")
    end
end

"""
    withenv(f::Function, kv::Pair...)

Execute `f()` in an environment that is temporarily modified (not replaced as in `setenv`)
by zero or more `"var"=>val` arguments `kv`. `withenv` is generally used via the
`withenv(kv...) do ... end` syntax. A value of `nothing` can be used to temporarily unset an
environment variable (if it is set). When `withenv` returns, the original environment has
been restored.
"""
function withenv{T<:AbstractString}(f::Function, keyvals::Pair{T}...)
    old = Dict{T,Any}()
    for (key,val) in keyvals
        old[key] = get(ENV,key,nothing)
        val !== nothing ? (ENV[key]=val) : delete!(ENV, key)
    end
    try f()
    finally
        for (key,val) in old
            val !== nothing ? (ENV[key]=val) : delete!(ENV, key)
        end
    end
end
withenv(f::Function) = f() # handle empty keyvals case; see #10853


text:
#=
array3d_to_lar:
- Julia version:
- Author: miros
- Date: 2019-01-14
=#
include("../src/LarSurf.jl")
# include("../src/plasm.jl")
# include("../src/read.jl")
# include("../src/surface.jl")
# include("../src/representation.jl")
# include("../src/import3d.jl")
# LarSurf.

## Artifical data
data3d = zeros(Int16,5,4,3)
data3d[2:4,2:3,2] .= 10
# data3d

voxelsize_mm = [0.5, 0.9, 0.8]
threshold=0

# Get vertices and count the voxesl (*4 face number)
verts, nvoxels = LarSurf.vertices_and_count_voxels_as_the_square_face_preprocessing(data3d, voxelsize_mm)

# Create square faces

faces = LarSurf.create_square_faces(data3d, nvoxels, threshold)

new_verts, new_faces = LarSurf.keep_surface_faces(verts, faces)

verts = new_verts
faces = new_faces

# Triangulation

trifaces = LarSurf.triangulation(faces)

# Visualization
LarSurf.check_vf(verts, trifaces)

V, EV, FE = LarSurf.to_lar(verts, trifaces)
VV = [i for i=1:length(V)]

EVc = Lar.cop2lar(EV)
Vnew = convert(Array, transpose(V))
# LarSurf.createVVll()

Plasm.view((Vnew, [EVc]))

# Plasm.visualize_numbers
# this does not work yet
# LarSurf.visualize_numbers((V, EV, FE), 0.5)


text:
# reference implementation of a CPU-based array type

module JLArrays

using GPUArrays

export JLArray

#
# Host array
#

# the definition of a host array type, implementing different Base interfaces
# to make it function properly and behave like the Base Array type.

struct JLArray{T, N} <: AbstractGPUArray{T, N}
    data::Array{T, N}
    dims::Dims{N}

    function JLArray{T,N}(data::Array{T, N}, dims::Dims{N}) where {T,N}
        new(data, dims)
    end
end

## constructors

# type and dimensionality specified, accepting dims as tuples of Ints
JLArray{T,N}(::UndefInitializer, dims::Dims{N}) where {T,N} =
  JLArray{T,N}(Array{T, N}(undef, dims), dims)

# type and dimensionality specified, accepting dims as series of Ints
JLArray{T,N}(::UndefInitializer, dims::Integer...) where {T,N} = JLArray{T,N}(undef, dims)

# type but not dimensionality specified
JLArray{T}(::UndefInitializer, dims::Dims{N}) where {T,N} = JLArray{T,N}(undef, dims)
JLArray{T}(::UndefInitializer, dims::Integer...) where {T} =
  JLArray{T}(undef, convert(Tuple{Vararg{Int}}, dims))

# empty vector constructor
JLArray{T,1}() where {T} = JLArray{T,1}(undef, 0)

Base.similar(a::JLArray{T,N}) where {T,N} = JLArray{T,N}(undef, size(a))
Base.similar(a::JLArray{T}, dims::Base.Dims{N}) where {T,N} = JLArray{T,N}(undef, dims)
Base.similar(a::JLArray, ::Type{T}, dims::Base.Dims{N}) where {T,N} = JLArray{T,N}(undef, dims)

## array interface

Base.elsize(::Type{<:JLArray{T}}) where {T} = sizeof(T)

Base.size(x::JLArray) = x.dims
Base.sizeof(x::JLArray) = Base.elsize(x) * length(x)

## interop with other arrays

JLArray{T,N}(x::AbstractArray{S,N}) where {T,N,S} =
    JLArray{T,N}(convert(Array{T}, x), size(x))

# underspecified constructors
JLArray{T}(xs::AbstractArray{S,N}) where {T,N,S} = JLArray{T,N}(xs)
(::Type{JLArray{T,N} where T})(x::AbstractArray{S,N}) where {S,N} = JLArray{S,N}(x)
JLArray(A::AbstractArray{T,N}) where {T,N} = JLArray{T,N}(A)

# idempotency
JLArray{T,N}(xs::JLArray{T,N}) where {T,N} = xs

## conversions

Base.convert(::Type{T}, x::T) where T <: JLArray = x

## broadcast

using Base.Broadcast: BroadcastStyle, Broadcasted, ArrayStyle

BroadcastStyle(::Type{<:JLArray}) = ArrayStyle{JLArray}()

function Base.similar(bc::Broadcasted{ArrayStyle{JLArray}}, ::Type{T}) where T
    similar(JLArray{T}, axes(bc))
end

Base.similar(bc::Broadcasted{ArrayStyle{JLArray}}, ::Type{T}, dims...) where {T} = JLArray{T}(undef, dims...)

## memory operations

function Base.copyto!(dest::Array{T}, d_offset::Integer,
                      source::JLArray{T}, s_offset::Integer,
                      amount::Integer) where T
    @boundscheck checkbounds(dest, d_offset+amount-1)
    @boundscheck checkbounds(source, s_offset+amount-1)
    copyto!(dest, d_offset, source.data, s_offset, amount)
end

function Base.copyto!(dest::JLArray{T}, d_offset::Integer,
                      source::Array{T}, s_offset::Integer,
                      amount::Integer) where T
    @boundscheck checkbounds(dest, d_offset+amount-1)
    @boundscheck checkbounds(source, s_offset+amount-1)
    copyto!(dest.data, d_offset, source, s_offset, amount)
    dest
end

function Base.copyto!(dest::JLArray{T}, d_offset::Integer,
                      source::JLArray{T}, s_offset::Integer,
                      amount::Integer) where T
    @boundscheck checkbounds(dest, d_offset+amount-1)
    @boundscheck checkbounds(source, s_offset+amount-1)
    copyto!(dest.data, d_offset, source.data, s_offset, amount)
    dest
end

## fft

using AbstractFFTs

# defining our own plan type is the easiest way to pass around the plans in FFTW interface
# without ambiguities

struct FFTPlan{T}
    p::T
end

AbstractFFTs.plan_fft(A::JLArray; kw_args...) = FFTPlan(plan_fft(A.data; kw_args...))
AbstractFFTs.plan_fft!(A::JLArray; kw_args...) = FFTPlan(plan_fft!(A.data; kw_args...))
AbstractFFTs.plan_bfft!(A::JLArray; kw_args...) = FFTPlan(plan_bfft!(A.data; kw_args...))
AbstractFFTs.plan_bfft(A::JLArray; kw_args...) = FFTPlan(plan_bfft(A.data; kw_args...))
AbstractFFTs.plan_ifft!(A::JLArray; kw_args...) = FFTPlan(plan_ifft!(A.data; kw_args...))
AbstractFFTs.plan_ifft(A::JLArray; kw_args...) = FFTPlan(plan_ifft(A.data; kw_args...))

function Base.:(*)(plan::FFTPlan, A::JLArray)
    x = plan.p * A.data
    JLArray(x)
end

#
# AbstractGPUArray interface
#

# implementation of GPUArrays-specific interfaces

GPUArrays.unsafe_reinterpret(::Type{T}, A::JLArray, size::Tuple) where T =
    reshape(reinterpret(T, A.data), size)

## execution

struct JLBackend <: AbstractGPUBackend end

GPUArrays.backend(::Type{<:JLArray}) = JLBackend()

mutable struct JLState{N}
    blockdim::NTuple{N, Int}
    griddim::NTuple{N, Int}

    blockidx::NTuple{N, Int}
    threadidx::NTuple{N, Int}
    localmem_counter::Int
    localmems::Vector{Vector{Array}}
end

function JLState(threads::NTuple{N}, blockdim::NTuple{N}) where N
    idx = ntuple(i-> 1, Val(N))
    blockcount = prod(blockdim)
    lmems = [Vector{Array}() for i in 1:blockcount]
    JLState{N}(threads, blockdim, idx, idx, 0, lmems)
end

function JLState(state::JLState{N}, threadidx::NTuple{N}) where N
    JLState{N}(
        state.blockdim,
        state.griddim,
        state.blockidx,
        threadidx,
        0,
        state.localmems
    )
end

to_device(state, x::JLArray{T,N}) where {T,N} = JLDeviceArray{T,N}(x.data, x.dims)
to_device(state, x::Tuple) = to_device.(Ref(state), x)
to_device(state, x::Base.RefValue{<: JLArray}) = Base.RefValue(to_device(state, x[]))
to_device(state, x) = x

function GPUArrays._gpu_call(::JLBackend, f, A, args::Tuple, blocks_threads::Tuple{T, T}) where T <: NTuple{N, Integer} where N
    blocks, threads = blocks_threads
    idx = ntuple(i-> 1, length(blocks))
    blockdim = blocks
    state = JLState(threads, blockdim)
    device_args = to_device.(Ref(state), args)
    tasks = Array{Task}(undef, threads...)
    for blockidx in CartesianIndices(blockdim)
        state.blockidx = blockidx.I
        for threadidx in CartesianIndices(threads)
            thread_state = JLState(state, threadidx.I)
            tasks[threadidx] = @async @allowscalar f(thread_state, device_args...)
            # TODO: require 1.3 and use Base.Threads.@spawn for actual multithreading
            #       (this would require a different synchronization mechanism)
        end
        for t in tasks
            fetch(t)
        end
    end
    return
end

## device properties

struct JLDevice end

GPUArrays.device(x::JLArray) = JLDevice()

GPUArrays.threads(dev::JLDevice) = 256

## linear algebra

using LinearAlgebra

GPUArrays.blas_module(::JLArray) = LinearAlgebra.BLAS
GPUArrays.blasbuffer(A::JLArray) = A.data

#
# Device array
#

# definition of a minimal device array type that supports the subset of operations
# that are used in GPUArrays kernels

struct JLDeviceArray{T, N} <: AbstractDeviceArray{T, N}
    data::Array{T, N}
    dims::Dims{N}

    function JLDeviceArray{T,N}(data::Array{T, N}, dims::Dims{N}) where {T,N}
        new(data, dims)
    end
end

function GPUArrays.LocalMemory(state::JLState, ::Type{T}, ::Val{dims}, ::Val{id}) where {T, dims, id}
    state.localmem_counter += 1
    lmems = state.localmems[blockidx_x(state)]

    # first invocation in block
    data = if length(lmems) < state.localmem_counter
        lmem = fill(zero(T), dims)
        push!(lmems, lmem)
        lmem
    else
        lmems[state.localmem_counter]
    end

    N = length(dims)
    JLDeviceArray{T,N}(data, tuple(dims...))
end

## array interface

Base.size(x::JLDeviceArray) = x.dims

## indexing

@inline Base.getindex(A::JLDeviceArray, index::Integer) = getindex(A.data, index)
@inline Base.setindex!(A::JLDeviceArray, x, index::Integer) = setindex!(A.data, x, index)

for (i, sym) in enumerate((:x, :y, :z))
    for f in (:blockidx, :blockdim, :threadidx, :griddim)
        fname = Symbol(string(f, '_', sym))
        @eval GPUArrays.$fname(state::JLState) = Int(state.$f[$i])
    end
end

## synchronization

@inline function GPUArrays.synchronize_threads(::JLState)
    # All threads are getting started asynchronously, so a yield will yield to the next
    # execution of the same function, which should call yield at the exact same point in the
    # program, leading to a chain of yields effectively syncing the tasks (threads).
    yield()
    return
end

end


text:
using Mimi
using Test
using DataFrames
using VegaLite
using Electron

import Mimi:
    dataframe_or_scalar, _spec_for_item, menu_item_list, getdataframe, dim_names

# Helper function returns true if VegaLite is verison 3 or above, and false otherwise
function _is_VegaLite_v3()
    return isdefined(VegaLite, :vlplot) ? true : false
end

@defcomp MyComp begin
    a = Parameter(index=[time, regions])
    b = Parameter(index=[time])
    c = Parameter(index=[regions])
    d = Parameter()
    e = Parameter(index=[four])
    f = Parameter{Array{Float64, 2}}()

    x = Variable(index=[time, regions])

    function run_timestep(p, v, d, t)
        for r in d.regions
            v.x[t, r] = 0
        end
    end
end

m = Model()
set_dimension!(m, :time, 2000:2100)
set_dimension!(m, :regions, 3)
set_dimension!(m, :four, 4)

add_comp!(m, MyComp)
update_param!(m, :MyComp, :a, ones(101,3))
update_param!(m, :MyComp, :b, 1:101)
update_param!(m, :MyComp, :c, [4,5,6])
update_param!(m, :MyComp, :d, .5)
update_param!(m, :MyComp, :e, [1,2,3,4])
update_param!(m, :MyComp, :f, [1.0 2.0; 3.0 4.0])

run(m)

# 1.  dataframe helper functions
@test typeof(dataframe_or_scalar(m, :MyComp, :a)) == DataFrame
@test typeof(dataframe_or_scalar(m, :MyComp, :d)) == Float64

# 2.  Specs and menu
items = [:a, :b, :c, :d, :e, :f, :x]
for item in items
    static_spec = _spec_for_item(m, :MyComp, item; interactive = false)
    interactive_spec = _spec_for_item(m, :MyComp, item)
    if length(dim_names(m, :MyComp, item)) == 0
        name =  string(:MyComp, " : ", item, " = ", m[:MyComp, item])
    else
        name = string(:MyComp, " : ", item)
    end
    @test static_spec["name"] == interactive_spec["name"] == name
end

s = menu_item_list(m)
@test collect(keys(s)) == ["pars", "vars"]
@test length(collect(keys(s["pars"]))) == 6
@test length(collect(keys(s["vars"]))) == 1

# 3.  explore(m::Model)
w = explore(m)
@test typeof(w) == Electron.Window
close(w)

# 4.  Mimi.plot(m::Model, comp_name::Symbol, datum_name::Symbol;
#       dim_name::Union{Nothing, Symbol} = nothing)
items = [:a, :b, :c, :d, :e, :f, :x]
for item in items
    p_type = _is_VegaLite_v3() ? VegaLite.VLSpec : VegaLite.VLSpec{:plot}
    @test typeof(Mimi.plot(m, :MyComp, item)) == p_type
end

# 5.  errors and warnings
@defcomp MyComp2 begin

    a = Parameter(index = [time, regions, four])
    x = Variable(index=[time, regions])

    function run_timestep(p, v, d, t)
        for r in d.regions
            v.x[t, r] = rand(10)[1]
        end
    end
end

m2 = Model()
set_dimension!(m2, :time, 2000:2100)
set_dimension!(m2, :regions, 3)
set_dimension!(m2, :four, 4)

add_comp!(m2, MyComp2)
update_param!(m2, :MyComp2, :a, ones(101, 3, 4))

run(m2)

# spec creation for MyComp.a should warn because over 2 indexed dimensions
# @test_logs (:warn, "MyComp2.a has > 2 indexed dimensions, not yet implemented in explorer") explore(m2) # URI PARSER WARNING IN CI
# @test_logs (:warn, "MyComp2.a has > 2 indexed dimensions, not yet implemented in explorer") _spec_for_item(m2, :MyComp2, :a) # URI PARSER WARNING IN CI

#7. Test TimestepArrays with time not as the first dimension

@defcomp gdp begin
    growth = Parameter(index=[regions, foo, time, 2])   # test that time is not first but not last
    gdp = Variable(index=[regions, foo, time, 2])
    gdp0 = Parameter(index=[regions, foo, 2])

    pgrowth = Parameter(index=[regions, 3, time])       # test time as last
    pop = Variable(index=[regions, 3, time])

    mat = Parameter(index=[regions, time])              # test time as last for a matrix
    mat2 = Variable(index=[regions, time])

    function run_timestep(p, v, d, ts)
        if is_first(ts)
            v.gdp[:, :, ts, :] = (1 .+ p.growth[:, :, ts, :]) .* p.gdp0
            v.pop[:, :, ts] = zeros(2, 3)
        else
            v.gdp[:, :, ts, :] = (1 .+ p.growth[:, :, ts, :]) .* v.gdp[:, :, ts-1, :]
            v.pop[:, :, ts] = v.pop[:, :, ts-1] .+ p.pgrowth[:, :, ts]
        end
        v.mat2[:, ts] = p.mat[:, ts]
    end
end

time_index = 2000:2100
regions = ["OECD","non-OECD"]
nsteps=length(time_index)

m = Model()
set_dimension!(m, :time, time_index)
set_dimension!(m, :regions, regions)
set_dimension!(m, :foo, 3)
add_comp!(m, gdp)
update_param!(m, :gdp, :gdp0, [3; 7] .* ones(length(regions), 3, 2))
update_param!(m, :gdp, :growth, [0.02; 0.03] .* ones(length(regions), 3, nsteps, 2))
set_leftover_params!(m, Dict{String, Any}([
    "pgrowth" => ones(length(regions), 3, nsteps),
    "mat" => rand(length(regions), nsteps)
]))
run(m)
w = explore(m)
close(w)

# 8. More plot types

@defcomp example begin
    p0 = Parameter(index = [time]) # line plot
    p1 = Parameter(index = [foo]) # scatter plot
    p2 = Parameter(index = [baz]) # bar plot
    p3 = Parameter(index = [foo, bar]) # multi-scatter
    p4 = Parameter(index = [baz, foo]) # multi-bar
    p5 = Parameter(index = [time, bar]) # multi-line
    p6 = Parameter(index = [time, baz]) # multi-line

    x = Variable(index = [time])
    function run_timestep(p,v,d,t)
        v.x[t] = 0
    end
end

m = Model()

set_dimension!(m, :time, 2000:2009)
set_dimension!(m, :foo, 1:5)
set_dimension!(m, :bar, 1:3)
set_dimension!(m, :baz, [:A, :B, :C])

add_comp!(m, example)

update_param!(m, :example, :p0, 1:10)
update_param!(m, :example, :p1, 6:10)
update_param!(m, :example, :p2, 4:6)

update_param!(m, :example, :p3, reshape(1:15, 5, 3))
update_param!(m, :example, :p4, reshape(1:15, 3, 5))
update_param!(m, :example, :p5, reshape(1:30, 10, 3))
update_param!(m, :example, :p6, reshape(1:30, 10, 3))

run(m)
explore(m)

items = [:p0, :p1, :p2, :p3, :p4, :p5, :p6]
for item in items
    static_spec = _spec_for_item(m, :example, item; interactive = false)
    interactive_spec = _spec_for_item(m, :example, item)
    if length(dim_names(m, :example, item)) == 0
        name =  string(:example, " : ", item, " = ", m[:example, item])
    else
        name = string(:example, " : ", item)
    end
    @test static_spec["name"] == interactive_spec["name"] == name
end


text:
#=Given a posititve number n find all the numbers from 2 to n which
are prime and print them. Sieve of Eratosthenes is an optimal approach
to find the prime numbers.We get an array of length n and if the value
of array at index i is zero the number i is prime.=#

## Function

function primeSieve(n)
    a = zeros(Int64, n + 5)
    for i = 2:n
        if (a[i] == 0)
            for j = (i*i):i:n
                a[j] = 1
            end
        end
    end
    for i = 2:n
        if (a[i] == 0)
            print("$(i) ")
        end
    end
end

##Input

n = readline()
n = parse(Int64, n)

#Calling the function

primeSieve(n)

#=
Sample test case:
Input:
    n = 20
Output:
    2 3 5 7 11 13 17 19

Time complexity: O( N * log(log(N)) )
=#


text:
"""
    find_ccalls(f, tt)

Returns a `Dict` mapping function addresses to symbol names for all `ccall`s and
`cglobal`s called from the method. This descends into other invocations
within the method.
"""
find_ccalls(@nospecialize(f), @nospecialize(tt)) = find_ccalls(reflect(f, tt))

function find_ccalls(ref::Reflection)
    result = Dict{Ptr{Nothing}, Symbol}()
    idx = VERSION > v"1.2" ? 5 : 4
    foreigncalls = TypedCodeUtils.filter((c) -> lookthrough((c) -> c.head === :foreigncall && !(c.args[idx] isa QuoteNode && c.args[idx].value == :llvmcall), c), ref.CI.code)
    # foreigncalls = TypedCodeUtils.filter((c) -> lookthrough((c) -> c.head === :foreigncall, c), ref.CI.code)
    for fc in foreigncalls
        sym = getsym(fc[2].args[1])
        address = eval(:(cglobal($(sym))))
        result[address] = Symbol(sym isa Tuple ? sym[1] : sym.value)
    end
    cglobals = TypedCodeUtils.filter((c) -> lookthrough(c -> c.head === :call && iscglobal(c.args[1]), c), ref.CI.code)
    for fc in cglobals
        sym = getsym(fc[2].args[2])
        address = eval(:(cglobal($(sym))))
        result[address] = Symbol(sym isa Tuple ? sym[1] : sym.value)
    end
    invokes = TypedCodeUtils.filter((c) -> lookthrough(identify_invoke, c), ref.CI.code)
    invokes = map((arg) -> process_invoke(DefaultConsumer(), ref, arg...), invokes)
    for fi in invokes
        canreflect(fi) || continue
        merge!(result, find_ccalls(reflect(fi)))
    end
    return result
end

getsym(x) = x
getsym(x::String) = QuoteNode(Symbol(x))
getsym(x::QuoteNode) = x
getsym(x::Expr) = eval.((x.args[2], x.args[3]))

iscglobal(x) = x == cglobal || x isa GlobalRef && x.name == :cglobal

"""
    fix_ccalls!(mod::LLVM.Module, d)

Replace function addresses with symbol names in `mod`. The symbol names are
meant to be linked to `libjulia` or other libraries.
`d` is a `Dict` mapping a function address to symbol name for `ccall`s.
"""
function fix_ccalls!(mod::LLVM.Module, d)
    for fun in functions(mod), blk in blocks(fun), instr in instructions(blk)
        if instr isa LLVM.CallInst
            dest = called_value(instr)
            if dest isa ConstantExpr && occursin("inttoptr", string(dest))
                # @show instr
                # @show dest
                argtypes = [llvmtype(op) for op in operands(instr)]
                nargs = length(parameters(eltype(argtypes[end])))
                # num_extra_args = 1 + length(collect(eachmatch(r"jl_roots", string(instr))))
                ptr = Ptr{Cvoid}(convert(Int, first(operands(dest))))
                if haskey(d, ptr)
                    sym = d[ptr]
                    newdest = LLVM.Function(mod, string(sym), LLVM.FunctionType(llvmtype(instr), argtypes[1:nargs]))
                    LLVM.linkage!(newdest, LLVM.API.LLVMExternalLinkage)
                    replace_uses!(dest, newdest)
                end
            end
        elseif instr isa LLVM.LoadInst && occursin("inttoptr", string(instr))
            # dest = called_value(instr)
            for op in operands(instr)
                lastop = op
                if occursin("inttoptr", string(op))
                    # @show instr
                    if occursin("addrspacecast", string(op)) || occursin("getelementptr", string(op))
                        op = first(operands(op))
                    end
                    first(operands(op)) isa LLVM.ConstantInt || continue
                    ptr = Ptr{Cvoid}(convert(Int, first(operands(op))))
                    if haskey(d, ptr)
                        obj = d[ptr]
                        newdest = GlobalVariable(mod, llvmtype(instr), string(d[ptr]))
                        LLVM.linkage!(newdest, LLVM.API.LLVMExternalLinkage)
                        replace_uses!(op, newdest)
                    end
                end
            end
        end
    end
end


text:
module BitsMNIST

	using Reexport
	@reexport using Flux
	@reexport using TinyML

	include("./Modules/Misc.jl")
	using .Misc
	include("./Modules/Datasets.jl")
	using .Datasets
	include("./Modules/IO.jl")
	using .IO
	include("./Modules/Statistics.jl")
	using .Statistics
	include("./Modules/ZeroOne.jl")
	using .ZeroOne

end


text:
using DelimitedFiles
using Printf

# Tests all passports in file with requested keys
# Call function with valcheck=true to test values
function test_passports(file, req_keys; valcheck=false)
    # Initialize variables
    total = 0
    inval = 0

    # Open the files and read through lines
    open(file) do f
        # Empty dict
        dict = Dict()

        # While not at the end of file
        while !eof(f)
            line = readline(f) # read line

            if line != "" # if line is not empty
                # Split by space and then by key:val pairs
                pairs = split(line, " ")

                for p in pairs
                    a,b = split(p, ":")
                    dict[a] = b # add value to dictionary
                end

            # If line is empty, assume passport entry is over.
            else
                total += 1 # increment total passports

                # Check if any keys are missing
                if check_keys(dict, req_keys) == false
                    inval += 1

                # Check if values are valid (if valcheck == true)
                elseif valcheck && check_vals(dict) == false
                    inval += 1
                end

                # Empty dict
                dict = Dict()
            end
        end

        # Do a check one more time for last passport entry
        # (last line doesn't catch the [if line != ""] statement)
        if check_keys(dict, req_keys) == false
            inval += 1
        elseif valcheck && check_vals(dict) == false
            inval += 1
        end
        total +=1
    end

    # Print output
    println(file)
    println("Total\t",total)
    println("Invalid\t",inval)
    println("Valid\t",total-inval)
end

# Function to loop through all req_keys and check for existence
function check_keys(dict, req_keys)
    for key in req_keys
        if haskey(dict, key) == false
            return false
        end
    end
    return true
end

# Function to check values for correctness
function check_vals(dict)

    # Regex matching for each key
    tests = [
        ["byr", r"^(19[2-9]\d|20[0-1]\d|2020)$"],
        ["iyr", r"^(20[1]\d|2020)$"],
        ["eyr", r"^(202\d|2030)$"],
        ["hgt", r"^(1[5-8]\d|19[0-3])cm|(59|6\d|7[0-6])in$"],
        ["hcl", r"^#[a-f0-9]{6}$"],
        ["ecl", r"^(amb|blu|brn|gry|grn|hzl|oth)$"],
        ["pid", r"^[0-9]{9}$"],
    ]

    # Check each value and return false if not valid
    for t in tests
        key = t[1]
        test = t[2]
        if occursin(test, dict[key]) == false
            return false
        end
    end
    return true
end

req_keys = [
    "byr",
    "iyr",
    "eyr",
    "hgt",
    "hcl",
    "ecl",
    "pid",
]

test_passports("./input.dat", req_keys)
println("============")
test_passports("./input.dat", req_keys, valcheck=true)


text:
push!(LOAD_PATH, @__DIR__)

using NodeJS, FF
import JSON
import BitIntegers: UInt256, @uint256_str

include("exec.jl")

path = nodejs_cmd()

function parse(src::String)
  astRep = read(Cmd([nodejs_cmd().exec[1], "parser/parser.js", src]), String);
  return JSON.parse(astRep)
end

function compile(srcFile)
  ast = parse(srcFile)
  @assert ast["type"]=="BLOCK" "ast.type must be BLOCK"
  ctx = Ctx()
  exec!(ctx, ast)
end

compile("localtests/circuit.circom");


text:
module PerformanceProfilingHttpEndpoints

import HTTP
import Profile
import PProf

using Serialization: serialize

#----------------------------------------------------------
#  Description
#
# It would be nice if visiting the `/profile` page from a web browser gives you a form where you can:
# - configure the Profiling parameters:
#    - Set julia Profile configuration, either by:
#       - Manually setting `n` and `delay`, for `Profile.init(n, delay)`, or
#       - provide a `duration` and a `delay`, and we will compute suggested `n` and `delay` for you. (`n` is the number of stack _frames_  the profiling will hold, so it also depends on the expected stack _depth_ for each sample, and the number of threads currently running. It's good to leave a wide buffer.)
# - Export the profiling data (either with or without C-lang functions - default is `C=true`)
#   - Download the raw Julia profile data (`Profile.retrieve()`)
#   - Convert the profile data to PProf, and download it (`PProf.pprof()`)
#   - Upload the resulting profile to DataDog (as PProf profile) -- see #4937.
#
#----------------------------------------------------------

default_n() = "1e8"
default_delay() = "0.01"
default_duration() = "10.0"
default_pprof() = "true"
default_alloc_sample_rate() = "0.0001"

cpu_profile_error_message() = """Need to provide query params:
    - duration=$(default_duration())
    - delay=$(default_delay())
    - n=$(default_n())
    - pprof=$(default_pprof())

Hint: A good goal is to shoot for around 1,000 to 10,000 samples. So if you know what
duration you want to profile for, you can pick a delay via e.g. `delay = duration / 1,000`.

For example, for a 30 second profile:
30 / 1_000 = 0.03
for `duration=30&delay=0.03`

But note that if your delay gets too small, you can really slow down the program you're
profiling, and thus end up with an inaccurate profile.

If you wanted to get an extremely large number of samples (or you have very very deep stack
traces), you'll have to start worrying about filling up your julia profiling buffer,
controlled by `n=`. If you assume an average stack depth of 100, and you were aiming to get
50,000 samples, you'd need a buffer at least 50,000 * 100 big, or at least 5e6.

The default `n` is 1e8, which should be big enough for most profiles.
"""
allocs_profile_error_message() = """Need to provide query params:
    - duration=$(default_duration())
    - sample_rate=$(default_alloc_sample_rate())

Hint: A good goal is to shoot for around 1,000 to 10,000 samples. So if you know what
duration you want to profile for, and you *already have an expectation for how much your
program will allocate,* you can pick a sample_rate via `sample_rate = 1,000 / expected_allocations`.

For example, if you expect your program will actually perform 1 million allocations:
1_000 / 1_000_000 = 0.001
for `duration=30&sample_rate=0.001`

Note that if your sample_rate gets too large, you can really slow down the program you're
profiling, and thus end up with an inaccurate profile.

Finally, if you think your program only allocates a small amount, you can capture *all*
allocations by passing sample_rate=1.
"""

function cpu_profile_endpoint(req::HTTP.Request)
    uri = HTTP.URI(HTTP.Messages.uri(req))
    qp = HTTP.queryparams(uri)
    if isempty(qp)
        @info "TODO: interactive HTML input page"
        return HTTP.Response(400, cpu_profile_error_message())
    end

    # Run the profile
    n = convert(Int, parse(Float64, get(qp, "n", default_n())))
    delay = parse(Float64, get(qp, "delay", default_delay()))
    duration = parse(Float64, get(qp, "duration", default_duration()))
    with_pprof = parse(Bool, get(qp, "pprof", default_pprof()))

    return _do_cpu_profile(n, delay, duration, with_pprof)
end

function _do_cpu_profile(n, delay, duration, with_pprof)
    @info "Starting CPU Profiling from PerformanceProfilingHttpEndpoints with configuration:" n delay duration

    Profile.clear()

    Profile.init(n, delay)

    Profile.@profile sleep(duration)

    data = Profile.retrieve()
    if with_pprof
        prof_name = tempname()
        PProf.pprof(out=prof_name, web=false)
        prof_name = "$prof_name.pb.gz"
        return _http_response(read(prof_name),
            "cpu_profile-duration=$duration&delay=$delay&n=$n.pb.gz")
    else
        iobuf = IOBuffer()
        serialize(iobuf, data)
        return _http_response(take!(iobuf),
            "cpu_profile&duration=$duration&delay=$delay&n=$n.prof.bin")
    end
end

function _http_response(binary_data, filename)
    return HTTP.Response(200, [
        "Content-Type" => "application/octet-stream"
        "Content-Disposition" => "attachment; filename=$(repr(filename))"
    ], body = binary_data)
end

function heap_snapshot_endpoint(req::HTTP.Request)
    # TODO: implement this once https://github.com/JuliaLang/julia/pull/42286 is merged
end

@static if !(isdefined(Profile, :Allocs) && isdefined(PProf, :Allocs))

function allocations_profile_endpoint(::HTTP.Request)
    return HTTP.Response(501, "You must use a build of Julia (1.8+) and PProf that support Allocations profiling.")
end

else

function allocations_profile_endpoint(req::HTTP.Request)

    uri = HTTP.URI(HTTP.Messages.uri(req))
    qp = HTTP.queryparams(uri)
    if isempty(qp)
        @info "TODO: interactive HTML input page"
        return HTTP.Response(400, allocs_profile_error_message())
    end

    sample_rate = convert(Float64, parse(Float64, get(qp, "sample_rate", default_alloc_sample_rate())))
    duration = parse(Float64, get(qp, "duration", default_duration()))

    return _do_alloc_profile(duration, sample_rate)
end

function _do_alloc_profile(duration, sample_rate)
    @info "Starting allocation Profiling from PerformanceProfilingHttpEndpoints with configuration:" duration sample_rate

    Profile.Allocs.clear()

    Profile.Allocs.@profile sample_rate=sample_rate sleep(duration)

    prof_name = tempname()
    PProf.Allocs.pprof(out=prof_name, web=false)
    prof_name = "$prof_name.pb.gz"
    return _http_response(read(prof_name),
            "allocs_profile-duration=$duration&sample_rate=$sample_rate.pb.gz")
end

end  # if isdefined

function serve_profiling_server(;addr="127.0.0.1", port=16825)
    @info "Starting HTTP profiling server on port $port"
    HTTP.serve(addr, port) do req
        # Invoke latest for easier development with Revise.jl :)
        Base.invokelatest(_server_handler, req)
    end
end

function _server_handler(req::HTTP.Request)
    @info "DEBUG REQUEST: $(HTTP.Messages.uri(req))"

    uri = HTTP.URI(HTTP.Messages.uri(req))
    segments = HTTP.URIs.splitpath(uri)
    @assert length(segments) >= 1
    path = segments[1]

    if path == "profile"
        return cpu_profile_endpoint(req)
    elseif path == "allocs_profile"
        return allocations_profile_endpoint(req)
    end

    @info "Unsupported Path: $path"

    return HTTP.Response(404)
end

# Precompile the endpoints as much as possible, so that your /profile attempt doesn't end
# up profiling compilation!
function __init__()
    precompile(serve_profiling_server, ()) || error("precompilation of package functions is not supposed to fail")
    precompile(cpu_profile_endpoint, (HTTP.Request,)) || error("precompilation of package functions is not supposed to fail")
    precompile(_do_cpu_profile, (Int,Float64,Float64,Bool)) || error("precompilation of package functions is not supposed to fail")
    precompile(allocations_profile_endpoint, (HTTP.Request,)) || error("precompilation of package functions is not supposed to fail")
    if isdefined(Profile, :Allocs) && isdefined(PProf, :Allocs)
        precompile(_do_alloc_profile, (Float64,Float64,)) || error("precompilation of package functions is not supposed to fail")
    end
end

end # module PerformanceProfilingHttpEndpoints

max_stars_repo_name: NHDaly/ServerIntrospectionEndpoints

max_issues_repo_name: NHDaly/ServerIntrospectionEndpoints

max_forks_repo_name: NHDaly/ServerIntrospectionEndpoints


text:
using Snowball
using Test

@testset "Snowball.jl" begin

    algs = stemmer_types()
    @test !isempty(algs)
    @test length(algs) >= 26
    @test "english" in algs

    for alg in algs
        stmr = Stemmer(alg)
        @test stmr.cptr != C_NULL
        Snowball.release(stmr)
    end

    test_cases = Dict{String, Any}(
        "english" => Dict{AbstractString, AbstractString}(
            "working" => "work",
            "worker" => "worker",
            "aβc" => "aβc",
            "a∀c" => "a∀c",
            "fishes" => "fish",
            "fishing" => "fish",
            "fish" => "fish",
            "fisher" => "fisher",
        )
    )

    for (alg, test_words) in test_cases
        stmr = Stemmer(alg)
        for (n,v) in test_words
            @test v == stem(stmr, n)
        end
    end
    stmr = Stemmer("english")
    @test  stem_all(stmr, "The fisher is fishing the fishes") == "The fisher is fish the fish"

    @test stem(stmr, [ "The", "fishers", "are", "fishing", "the", "fishes"]) ==
        [ "The", "fisher", "are", "fish", "the", "fish"]
end

max_stars_repo_name: JuliaText/Snowball


text:
function decompose_RX_on_q3()

    println(">>>>> RX Gate on third qubit using U3Gate <<<<<")

    return Dict{String, Any}(

    "num_qubits" => 3,
    "maximum_depth" => 3,
    "elementary_gates" => ["U3_3", "Identity"],
    "target_gate" => QCOpt.kron_single_qubit_gate(3, QCOpt.RXGate(π/4), "q3"),
    "objective" => "minimize_depth",
    "decomposition_type" => "exact",

    "U3_θ_discretization" => [0, π/4],
    "U3_ϕ_discretization" => [0, -π/2],
    "U3_λ_discretization" => [0, π/2],
    )

end

function decompose_toffoli()

    function toffoli_circuit()
        # [(depth, gate)]
        return [(1, "T_1"),
                (2, "T_2"),
                (3, "H_3"),
                (4, "CNot_2_3"),
                (5, "Tdagger_3"),
                (6, "CNot_1_3"),
                (7, "T_3"),
                (8, "CNot_2_3"),
                (9, "Tdagger_3"),
                (10, "CNot_1_3"),
                (11, "T_3"),
                (12, "H_3"),
                (13, "CNot_1_2"),
                (14, "Tdagger_2"),
                (15, "CNot_1_2")
                ]
    end

    println(">>>>> Toffoli gate <<<<<")

    return Dict{String, Any}(

    "num_qubits" => 3,
    "maximum_depth" => 15,
    "elementary_gates" => ["T_1", "T_2", "T_3", "H_3", "CNot_1_2", "CNot_1_3", "CNot_2_3", "Tdagger_1", "Tdagger_2", "Tdagger_3", "Identity"],
    "target_gate" => QCOpt.ToffoliGate(),
    "objective" => "minimize_depth",
    "decomposition_type" => "exact",

    "input_circuit" => toffoli_circuit(),
    "set_cnot_lower_bound" => 6,
    "set_cnot_upper_bound" => 6,
    )

end

function decompose_toffoli_using_kronecker()

    println(">>>>> Toffoli gate using Kronecker <<<<<")

    return Dict{String, Any}(

    "num_qubits" => 3,
    "maximum_depth" => 12,
    "elementary_gates" => ["T_3", "H_3", "CNot_1_2", "CNot_1_3", "CNot_2_3", "Tdagger_3", "I_1xT_2xT_3", "CNot_1_2xH_3", "T_1xTdagger_2xI_3", "Identity"],
    "target_gate" => QCOpt.ToffoliGate(),
    "objective" => "minimize_depth",
    "decomposition_type" => "exact",

    "set_cnot_lower_bound" => 6,
    )

end

function decompose_toffoli_with_controlled_gates()

    # Reference: https://doi.org/10.1109/TCAD.2005.858352
    println(">>>>> Toffoli gate with controlled gates <<<<<")

    return Dict{String, Any}(

    "num_qubits" => 3,
    "maximum_depth" => 5,
    "elementary_gates" => ["CV_1_3", "CV_2_3", "CV_1_2", "CVdagger_1_3", "CVdagger_2_3", "CVdagger_1_2", "CNot_2_1", "CNot_1_2", "Identity"],
    # "elementary_gates" => ["CV_1_3", "CV_2_3", "CVdagger_1_3", "CNot_1_2", "CNot_2_1", "Identity"],
    "target_gate" => QCOpt.ToffoliGate(),
    "objective" => "minimize_depth",
    "decomposition_type" => "exact",
    )

end

function decompose_CNot_1_3()

    return Dict{String, Any}(

    "num_qubits" => 3,
    "maximum_depth" => 8,
    # "elementary_gates" => ["CNot_1_2", "CNot_2_3", "Identity"],
    "elementary_gates" => ["H_1", "H_3", "H_2", "CNot_2_1", "CNot_3_2", "Identity"],
    "target_gate" => QCOpt.get_full_sized_gate("CNot_1_3", 3),
    "objective" => "minimize_depth",
    )

end

function decompose_FredkinGate()

    return Dict{String, Any}(
    "num_qubits" => 3,
    "maximum_depth" => 7,
    # Reference: https://doi.org/10.1103/PhysRevA.53.2855
    "elementary_gates" => ["CV_1_2", "CV_2_3", "CV_1_3", "CVdagger_1_2", "CVdagger_2_3", "CVdagger_1_3", "CNot_1_2", "CNot_3_2", "CNot_2_3", "CNot_1_3", "Identity"],
    "target_gate" => QCOpt.CSwapGate(), #also Fredkin
    "objective" => "minimize_depth"
    )

end

function decompose_toffoli_left()

    # Reference: https://arxiv.org/pdf/0803.2316.pdf

    function target_gate()
        H_3 = QCOpt.get_full_sized_gate("H_3", 3)
        Tdagger_3 = QCOpt.get_full_sized_gate("Tdagger_3", 3)
        T_3 = QCOpt.get_full_sized_gate("T_3", 3)
        cnot_23 = QCOpt.get_full_sized_gate("CNot_2_3", 3)
        CNot_1_3 = QCOpt.get_full_sized_gate("CNot_1_3", 3)

        return H_3 * cnot_23 * Tdagger_3 * CNot_1_3 * T_3 * cnot_23 * Tdagger_3
    end

    return Dict{String, Any}(

    "num_qubits" => 3,
    "maximum_depth" => 7,
    "elementary_gates" => ["H_3", "T_1", "T_2", "T_3", "Tdagger_1", "Tdagger_2", "Tdagger_3", "CNot_1_2", "CNot_2_3", "CNot_1_3", "Identity"],
    "target_gate" => target_gate(),
    "objective" => "minimize_depth"
    )

end

function decompose_toffoli_right()

    # Reference: https://arxiv.org/pdf/0803.2316.pdf

    function target_gate()
        H_3 = QCOpt.get_full_sized_gate("H_3", 3)
        Tdagger_2 = QCOpt.get_full_sized_gate("Tdagger_2", 3)
        T_1 = QCOpt.get_full_sized_gate("T_1", 3)
        T_2 = QCOpt.get_full_sized_gate("T_2", 3)
        T_3 = QCOpt.get_full_sized_gate("T_3", 3)
        CNot_1_2 = QCOpt.get_full_sized_gate("CNot_1_2", 3)
        CNot_1_3 = QCOpt.get_full_sized_gate("CNot_1_3", 3)

        return CNot_1_3 * T_2 * T_3 * CNot_1_2 * H_3 * T_1 * Tdagger_2 * CNot_1_2
    end

    return Dict{String, Any}(

    "num_qubits" => 3,
    "maximum_depth" => 8,
    "elementary_gates" => ["H_3", "T_1", "T_2", "T_3", "Tdagger_1", "Tdagger_2", "Tdagger_3", "CNot_1_2", "CNot_2_3", "CNot_1_3", "Identity"],
    "target_gate" => target_gate(),
    "objective" => "minimize_depth"
    )

end

function decompose_miller()
    # Reference: https://doi.org/10.1109/TCAD.2005.858352

    function target_gate()
        CV_1_3 = QCOpt.get_full_sized_gate("CV_1_3", 3)
        CV_2_3 = QCOpt.get_full_sized_gate("CV_2_3", 3)
        CVdagger_2_3 = QCOpt.get_full_sized_gate("CVdagger_2_3", 3)
        CNot_1_2 = QCOpt.get_full_sized_gate("CNot_1_2", 3)
        CNot_3_1 = QCOpt.get_full_sized_gate("CNot_3_1", 3)
        CNot_3_2 = QCOpt.get_full_sized_gate("CNot_3_2", 3)

        return CNot_3_1 * CNot_3_2 * CV_2_3 * CNot_1_2 * CVdagger_2_3 * CV_1_3 * CNot_3_1 * CNot_1_2
    end

    return Dict{String, Any}(

    "num_qubits" => 3,
    "maximum_depth" => 8,
    "elementary_gates" => ["CV_1_3", "CV_2_3", "CVdagger_2_3", "CNot_1_2", "CNot_3_1", "CNot_3_2", "Identity"],
    "target_gate" => target_gate(),
    "objective" => "minimize_depth"
    )

end

function decompose_relative_toffoli()
    #Reference: https://arxiv.org/pdf/1508.03273.pdf

    return Dict{String, Any}(

        "num_qubits" => 3,
        "maximum_depth" => 9,
        "elementary_gates" => ["H_3", "T_3", "Tdagger_3", "CNot_1_2", "CNot_2_3", "CNot_1_3", "Identity"],
        "target_gate" => QCOpt.RCCXGate(),
        "objective" => "minimize_depth",

        "set_cnot_upper_bound" => 3
        )

end

function decompose_margolus()
    #Reference: https://arxiv.org/pdf/quant-ph/0312225.pdf

    return Dict{String, Any}(

        "num_qubits" => 3,
        "maximum_depth" => 7,
        "elementary_gates" => ["RY_3", "CNot_1_2", "CNot_2_3", "CNot_1_3", "Identity"],
        "target_gate" => QCOpt.MargolusGate(),
        "objective" => "minimize_depth",

        "RY_discretization" => -π:π/4:π,
        # "set_cnot_lower_bound" => 3,
        # "set_cnot_upper_bound" => 3
        )

end

function decompose_CiSwap()
    #Reference: https://doi.org/10.1103/PhysRevResearch.2.033097

    return Dict{String, Any}(

        "num_qubits" => 3,
        "maximum_depth" => 12,
        "elementary_gates" => ["H_3", "T_3", "Tdagger_3", "CNot_3_2", "CNot_2_3", "CNot_1_3", "Identity"],
        "target_gate" => QCOpt.CiSwapGate(),
        "objective" => "minimize_depth",
        )

end


text:
# ---
# title: 1489. Find Critical and Pseudo-Critical Edges in Minimum Spanning Tree
# id: problem1489
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Depth-first Search, Union Find
# link: <https://leetcode.com/problems/find-critical-and-pseudo-critical-edges-in-minimum-spanning-tree/description/>
# hidden: true
# ---
#
# Given a weighted undirected connected graph with `n` vertices numbered from
# `0` to `n - 1`, and an array `edges` where `edges[i] = [ai, bi, weighti]`
# represents a bidirectional and weighted edge between nodes `ai` and `bi`. A
# minimum spanning tree (MST) is a subset of the graph's edges that connects all
# vertices without cycles and with the minimum possible total edge weight.
#
# Find _all the critical and pseudo-critical edges in the given graph 's minimum
# spanning tree (MST)_. An MST edge whose deletion from the graph would cause
# the MST weight to increase is called a  _critical edge_. On the other hand, a
# pseudo-critical edge is that which can appear in some MSTs but not all.
#
# Note that you can return the indices of the edges in any order.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/06/04/ex1.png)
#
#
#
#     Input: n = 5, edges = [[0,1,1],[1,2,1],[2,3,2],[0,3,2],[0,4,3],[3,4,3],[1,4,6]]
#     Output: [[0,1],[2,3,4,5]]
#     Explanation: The figure above describes the graph.
#     The following figure shows all the possible MSTs:
#     ![](https://assets.leetcode.com/uploads/2020/06/04/msts.png)
#     Notice that the two edges 0 and 1 appear in all MSTs, therefore they are critical edges, so we return them in the first list of the output.
#     The edges 2, 3, 4, and 5 are only part of some MSTs, therefore they are considered pseudo-critical edges. We add them to the second list of the output.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/06/04/ex2.png)
#
#
#
#     Input: n = 4, edges = [[0,1,1],[1,2,1],[2,3,1],[0,3,1]]
#     Output: [[],[0,1,2,3]]
#     Explanation: We can observe that since all 4 edges have equal weight, choosing any 3 edges from the given 4 will yield an MST. Therefore all 4 edges are pseudo-critical.
#
#
#
#
# **Constraints:**
#
#   * `2 <= n <= 100`
#   * `1 <= edges.length <= min(200, n * (n - 1) / 2)`
#   * `edges[i].length == 3`
#   * `0 <= ai < bi < n`
#   * `1 <= weighti <= 1000`
#   * All pairs `(ai, bi)` are **distinct**.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# OSSB (optimal sampling for structured bandits).
# Based on https://arxiv.org/abs/1711.00400.

abstract type OSSBOptimisationAlgorithm end

function optimise_graves_lai(::CombinatorialInstance{T}, algo::OSSBOptimisationAlgorithm, ::Dict{T, Float64}; kwargs...) where T
  error("OSSB optimiser $algo has not yet been implemented.")
end

struct OSSB <: Policy
  algo::OSSBOptimisationAlgorithm
  γ::Float64
  ε::Float64
end

mutable struct GravesLaiResults{T}
  objective::Float64
  points::Vector{Vector{T}} # Vector of solutions.
  weights::Vector{Float64} # Vector of weights, one per solution.
  delayed::Union{Function, Nothing}
end

function GravesLaiResults(objective::Float64, points::Vector{Vector{T}}, weights::Vector{Float64}) where T
  return GravesLaiResults{T}(objective, points, weights, nothing)
end

function GravesLaiResults(instance::CombinatorialInstance{T}, t_val::Dict{T, Float64}, objective::Float64) where T
  return GravesLaiResults{T}(objective, Vector{Vector{T}}(undef, 0), Vector{Float64}(undef, 0),
                             () -> _ossb_exact_convex_combination(copy(instance), t_val))
end

function Base.getproperty(obj::GravesLaiResults, sym::Symbol)
  # If one of the delayed field is requested, start the computations.
  if getfield(obj, :delayed) !== nothing && (sym === :points || sym === :weights)
    obj.points, obj.weights = obj.delayed()
    obj.delayed = nothing
  end

  # Return the requested value.
  return getfield(obj, sym)
end

mutable struct GravesLaiDetails <: PolicyDetails
  n_iterations::Int
  time_per_iteration::Vector{Float64} # ms
  objectives::Vector{Float64}
  gl_bounds::Vector{Float64} # Usually equal to the objective, but not always: the objective may include a penalisation.
end

GravesLaiDetails(time::Float64, objective::Float64) = GravesLaiDetails(1, [time], [objective], [objective])
GravesLaiDetails() = GravesLaiDetails(0, Float64[0.0], Float64[0.0], Float64[0.0])

mutable struct OSSBDetails <: PolicyDetails
  solver_time::Float64 # ms
  graves_lai::GravesLaiDetails
end

OSSBDetails() = OSSBDetails(0.0, GravesLaiDetails())

function choose_action(instance::CombinatorialInstance{T}, policy::OSSB, state::State{T}; with_trace::Bool=false) where T
  # Ensure that the state contains a counter for the number of times each solution is played (and not every arm).
  # This requirement is very OSSB-specific.
  # Do not pre-populate this data structure: in the worst case, it contains one entry per combinatorial solution.
  # OSSB also requires a counter of rounds where there has been no exploitation.
  if ! (:ossb in state.policy_extension)
    state.policy_extension[:ossb] = (0, Dict{Vector{T}, Int}())
  end

  # For initialisation, if some arms have never been tried, force them to be tried.
  # This should not be required for OSSB, but it improves running times for the first few iterations.
  # TODO: propose a helper method for this part (including solve_linear and trace), it's common to almost all bandit algorithms.
  if any(v == 0 for v in values(state.arm_counts))
    weights = Dict(arm => iszero(state.arm_counts[arm]) ? 1.0 : 0.0 for arm in keys(state.arm_counts))

    t0 = time_ns()
    sol = solve_linear(instance, weights)
    t1 = time_ns()

    if with_trace
      return sol, OSSBDetails((t1 - t0) / 1_000_000, GravesLaiDetails())
    else
      return sol
    end
  end

  # Solve OSSB's problem to get a distribution over actions to play (not really a probability distribution, however).
  t0 = time_ns()
  if with_trace
    gl, gl_trace = optimise_graves_lai(instance, policy.algo, state.arm_average_reward, with_trace=true)
  else
    gl = optimise_graves_lai(instance, policy.algo, state.arm_average_reward, with_trace=false)
  end
  t1 = time_ns()

  # Sometimes, the (0, 0...) vertex is found. For bandits, this is very bad: no arm is played, no reward can be
  # obtained, no exploration is performed.
  if any(iszero.(gl.points))
    idx = findall(iszero.gl.(points))
    deleteat!(gl.points, idx)
    deleteat!(gl.weights, idx)
  end

  # If the solution is zero, there is not much to do. This probably indicates a bug in the underlying solver.
  if length(gl.points) == 0
    println("Uh ho...")
    println(gl.points)
    println(gl.weights)

    t0 = time_ns()
    sol = solve_linear(instance, state.arm_average_reward)
    t1 = time_ns()

    if with_trace
      return sol, OSSBDetails((t1 - t0) / 1_000_000, gl_trace)
    else
      return sol
    end
  end

  # Exploit phase? Check if all solutions have already been played enough according to the new solution.
  enters_exploit = true
  unplayed_solution = nothing
  for (idx, point) in enumerate(gl.points)
    # Has this solution ever been played?
    if !(point in keys(state.policy_extension[:ossb][1]))
      enters_exploit = false
      unplayed_solution = point
      break
    end

    # Has this solution been played enough?
    if gl.weights[idx] * (1 + policy.γ) * log(state.round) <= state.policy_extension[:ossb][1][point]
      enters_exploit = false
      break
    end

    # Ignore the other solutions, as they have a zero weight: they already have been played enough, as the current
    # solution recommends to play them at least zero times.
  end

  if enters_exploit
    t0 = time_ns()
    sol = solve_linear(instance, state.arm_average_reward)
    t1 = time_ns()

    if with_trace
      return sol, OSSBDetails((t1 - t0) / 1_000_000, gl_trace)
    else
      return sol
    end
  end

  # Explore or sample phase?
  state.policy_extension[:ossb][2] += 1 # One more round without exploitation.

  if unplayed_solution !== nothing
    # At least one solution has never been played: 0 <= ε S, so play it right now.
    # Special case of the explore phase.
    if with_trace
      return unplayed_solution, OSSBDetails(0.0, gl_trace)
    else
      return unplayed_solution
    end
  end

  least_played = reduce((x, y) -> state.policy_extension[:ossb][1][x] ≤ state.policy_extension[:ossb][1][y] ? x : y, keys(state.policy_extension[:ossb][1]))

  # Explore.
  # This search can be limited to the points in the returned convex combination: if another solution was selected as
  # least played, the comparison would reduce to 0 <= something.
  if gl.weights[least_played] * (1 + policy.γ) * log(state.round) <= state.policy_extension[:ossb][1][least_played]
    if with_trace
      return least_played, OSSBDetails(0.0, gl_trace)
    else
      return least_played
    end
  end

  # Sample.
  # Again, this search can be limited to the points in the returned convex combination: the other solutions would
  # be associated to either 0/0 or x/0 with x > 0, and these values are either undefined or infinite.
  furthest_away = nothing
  furthest_away_obj = Inf

  t0 = time_ns()
  for (idx, point) in enumerate(gl.points)
    obj = state.policy_extension[:ossb][1][point] / gl.weights[idx]
    if obj <= furthest_away_obj
      furthest_away = point
      furthest_away_obj = obj
    end
  end
  t1 = time_ns()

  if with_trace
    return least_played, OSSBDetails((t1 - t0) / 1_000_000, gl_trace)
  else
    return least_played
  end
end


text:
#ДАНО: Робот - в произвольной клетке ограниченного прямоугольной рамкой поля без внутренних перегородок и маркеров.
#РЕЗУЛЬТАТ: Робот - в исходном положении в центре косого креста (в форме X) из маркеров.

function reverse(side)
    for i=0:3
        if side == HorizonSide(i)
            return HorizonSide(mod(i+2, 4))
        end
    end
end

function next_side(side)
    for i=0:3
        if side == HorizonSide(i)
            return HorizonSide(mod(i+3,4))
        end
    end
end

function krest_13(r)
    side = Nord
    for i in 1:4
        while (isborder(r,side)||isborder(r,next_side(side))) == false
            move!(r,side)
            move!(r,next_side(side))
            putmarker!(r)
        end
        side = reverse(side)
        while ismarker(r)
            move!(r,side)
            move!(r,next_side(side))
        end
        side = reverse(next_side(side))
    end
    putmarker!(r)
end

max_stars_repo_name: gardend/mirea-progs

max_issues_repo_name: gardend/mirea-progs

max_forks_repo_name: gardend/mirea-progs


text:
using Pkg; Pkg.activate("..")

using OceanTurb, Printf, StaticArrays, BenchmarkTools

Ns = (100,)
steppers = (:ForwardEuler,)
nts = (100,)

function manyrhs!(m, n)
    for i = 1:n
        OceanTurb.calc_explicit_rhs!(m.timestepper.rhs, m.timestepper.eqn, m.solution, m)
    end
    return nothing
end

function bad_manyrhs!(m, n)
    s1 = m.solution
    s2 = m.timestepper.rhs
    for i = 1:n
        for j in eachindex(s1)
            ϕ1 = s1[j]
            ϕ2 = s2[j]
            for i in eachindex(ϕ1)
                @inbounds ϕ1[i] = 2*ϕ2[i]
            end
        end
    end
    return nothing
end

function good_functionbarrier_manyrhs!(m, n)
    s1 = m.solution
    s2 = m.timestepper.rhs
    for i = 1:n
        fakerhs!(s1, s2)
    end
    return nothing
end

function fakerhs!(s1, s2)
    for j in eachindex(s1)
        ϕ1 = s1[j]
        ϕ2 = s2[j]
        for i in eachindex(ϕ1)
            @inbounds ϕ1[i] = 2*ϕ2[i]
        end
    end
    return nothing
end

@printf "\nTesting OceanTurb calc rhs...\n"
for N in Ns
    for stepper in (:ForwardEuler,)
        model = Diffusion.Model(N=N, stepper=stepper)
        for nt in nts
            @printf "stepper: % 16s, N: % 6d, nt: % 6d" stepper N nt
            manyrhs!(model, nt)
            @btime manyrhs!($model, $nt)
        end
    end
end

@printf "\nTesting a bad method to loop over rhs...\n"
for N in Ns
    for stepper in (:ForwardEuler,)
        model = Diffusion.Model(N=N, stepper=stepper)
        for nt in nts
            @printf "stepper: % 16s, N: % 6d, nt: % 6d" stepper N nt
            bad_manyrhs!(model, nt)
            @btime bad_manyrhs!($model, $nt)
        end
    end
end

@printf "\nTesting a good function barrier method to loop over rhs...\n"
for N in Ns
    for stepper in (:ForwardEuler,)
        model = Diffusion.Model(N=N, stepper=stepper)
        for nt in nts
            @printf "stepper: % 16s, N: % 6d, nt: % 6d" stepper N nt
            good_functionbarrier_manyrhs!(model, nt)
            @btime good_functionbarrier_manyrhs!($model, $nt)
        end
    end
end